## Code

```
# Load the required packages
library(ggplot2)
library(scales) #use the percent scale later
# Create a function for running of the bets.
<- function(p, n, t, start=100, gain, loss, ergodic=FALSE, absorbing=FALSE){
bet
#p is probability of a gain
#n is how many people in the simulation
#t is the number of coin flips simulated for each person
#start is the number of dollars each person starts with
#if ergodic=FALSE, gain and loss are the multipliers
#if ergodic=TRUE, gain and loss are the dollar amounts
#if absorbing=TRUE, zero wealth ends the series of flips for that person
<- as.data.frame(c(p, n, t, start, gain, loss, ergodic, absorbing))
params rownames(params) <- c("p", "n", "t", "start", "gain", "loss", "ergodic", "absorbing")
colnames(params) <- "value"
<- matrix(data = NA, nrow = t, ncol = n)
sim
if(ergodic==FALSE){
for (j in 1:n) {
<- start
x for (i in 1:t) {
<- rbinom(n=1, size=1, prob=p)
outcome ifelse(outcome==0, x <- x*loss, x <- x*gain)
<- x
sim[i,j]
}
}
}
if(ergodic==TRUE){
for (j in 1:n) {
<- start
x for (i in 1:t) {
<- rbinom(n=1, size=1, prob=p)
outcome ifelse(outcome==0, x <- x-loss, x <- x+gain)
<- x
sim[i,j] if(absorbing==TRUE){
if(x<0){
:t,j] <- 0
sim[ibreak
}
}
}
}
}
<- rbind(rep(start,n), sim) #placing the starting sum in the first row
sim <- cbind(seq(0,t), sim) #number each period
sim <- data.frame(sim)
sim colnames(sim) <- c("period", paste0("p", 1:n))
<- list(params=params, sim=sim)
sim
sim
}
# Simulate 10,000 people who accept a series of 1000 50:50 bets to win \$50 or lose \$40 from a starting wealth of \$100.
set.seed(20191215)
<- bet(p=0.5, n=10000, t=1000, gain=1.5, loss=0.6, ergodic=FALSE)
nonErgodic
# Create a function to generate summary statistics.
<- function(sim, t=100){
summaryStats
<- mean(as.matrix(sim$sim[(t+1),2:(sim$params[2,]+1)])) # mean wealth
meanW <- median(as.matrix(sim$sim[(t+1),2:(sim$params[2,]+1)])) # median wealth
medianW <- sum(sim$sim[(t+1),2:(sim$params[2,]+1)]<(sim$params[4,]/100)) #number who lost more than 99% of their wealth
num99 <- num99/sim$params[2,]*100 #percentage who lost more than 99% of their wealth
per99 <- sum(sim$sim[(t+1),2:(sim$params[2,]+1)]>sim$params[4,]) #number who gain
numGain <- numGain/sim$params[2,]*100 #percentage who gain
perGain <- sum(sim$sim[(t+1),2:(sim$params[2,]+1)]>(sim$params[4,]*100)) #number who increase their wealth more than 100-fold
num100 <- max(sim$sim[(t+1),2:(sim$params[2,]+1)]) #wealth of wealthiest person
winner <- winner / sum(sim$sim[(t+1),2:(sim$params[2,]+1)]) #wealth share of wealthiest person
winnerShare
<- data.frame(meanW = meanW, medianW = medianW, num99 = num99, per99 = per99, numGain = numGain, perGain = perGain, num100 = num100, winner = winner, winnerShare = winnerShare)
stats
}
<- summaryStats(nonErgodic, 100)
nonErgodicStats
# Create a function for plotting the average wealth of the population over a set number of periods.
<- function(sim, t=100){
averagePlot
<- ggplot(sim$sim[c(1:(t+1)),], aes(x=period)) +
basePlot labs(y = "Average Wealth ($)")
<- basePlot +
averagePlot geom_line(aes(y = rowMeans(sim$sim[c(1:(t+1)),2:(sim$params[2,]+1)])), color = 1, linewidth=1)
averagePlot
}
averagePlot(nonErgodic, 100)
```