# Gamblers Ruin with multiple coins tossed at the same time

I'm quite ashamed to be stuck with a Gambler's Ruin problem, I guess I'm missing some basic statistical intuition here:

Three fair coins tossed. Heads gets +1, tails -1, pay-offs are added and net pay-off added to equity. The 3 tosses are repeated 1000 times. Initial equity is 10 \$ . What is the probability of total ruin (within +/- 0.05 error)?

I simulated the problem as 3 iid coin tosses in one round which is then repeated, the same as it would be with a repeated one coin toss. My simulated probability of ruin converges to ca. 83%, while 100% would be the correct answer. The only hint I have is 'Flipping a coin in succession is different from flipping three concurrently from markov lens'. Could someone help me and explain?

Thanks!!

Tobi

import numpy as np

class GamblersRuin(object):
"""
Three fair coins tossed. Heads gets +1, tails -1, pay-offs are added and net pay-off
The 3 tosses are repeated 1000 times. Initial equity is 10 dollars
p: probability that gambler is successful/ wins at each round.
i: gambler's initial amount of money/reserves
"""

def __init__(self, p, init_bal):
self.p = p
self.init_bal = init_bal
self.bal = init_bal
self.q = 1 - self.p
self.realizations = np.array(self.init_bal)
self.simulation_results = []

def coin_toss(self):
"""
One coin flip with payoff (1, -1) with probability (p,q)
"""
outcome = np.random.uniform(0, 1)

if outcome < self.p:
result = 1
else:
result = -1

return result

def play_one_round(self):
"""
Three coin tosses in one round round
"""
result_round = 0
for i in range(0,3):
result_round += self.coin_toss()
return result_round

def gamble(self, no_rounds):
"""
One round is played until ruin or no_rounds times
"""
self.realizations = np.array(self.init_bal)
self.bal = self.init_bal

round = 1
while round < no_rounds:
round_result = self.play_one_round()
if (self.bal + round_result) >= 0:
self.bal += round_result
else:
break
self.realizations = np.append(self.realizations, self.bal)
round += 1

def simulate(self, no_simulations, no_rounds):
# Gamble multiple times and store realization paths
self.simulation_results = []

for game in range(1,no_simulations+1):
self.gamble(no_rounds=no_rounds)
self.simulation_results.append(self.realizations)


### Monte Carlo method

DyedPurple already showed that your simulation is not wrong and you should get a probability of ~0.84 for a run length of 1000. It is only when the run length goes towards infinity that you are almost certain to get gambler's ruin (If you have a stopping rule for some upper boundary, as in this question, then you can escape the gambler's ruin with some non-zero probability).

In this answer, I show how you can compute it exactly instead of simulating it with a Monte Carlo method. (and there is also an approximate analytic solution by comparing the situation with Brownian motion).

### Computation as a Markov chain

The problem is similar to this question Amoeba Interview Question or this question The Frog Problem (puzzle in YouTube video)

The probabilities, $$P_k(x)$$, to have $$x$$ money after $$k$$ tosses can be expressed in terms of the probabilities for earlier tosses:

$$P_k(x) = \frac{1}{8} P_{k-1}(x-3) + \frac{3}{8} P_{k-1}(x-1) + \frac{3}{8} P_{k-1}(x+1) + \frac{1}{8} P_{k-1}(x+3)$$

With this formula, you can already compute the result for 1000 steps (see the R-code and the image below).

### Comparison with a diffusion process

You can also model the amount of money as approximately a one-dimensional diffusion process or a Brownian motion (the solution is given in 1916 by Smoluchowski, more on that in the answer here https://stats.stackexchange.com/a/401539).

The amount of money $$M_k$$ in step $$k$$ changes relatively to the amount in the previous step $$M_{k-1}$$ by the addition of a random variable

$$M_k = M_{k-1} + \epsilon_k$$

In this case the random variable $$\epsilon_k$$ is a scaled and shifted binomial distributed variable that takes values $$-3$$, $$-1$$, $$1$$, $$3$$, with probabilities $$1/8$$, $$3/8$$, $$3/8$$, $$1/8$$. This variable has a variance equal to 3.

We can relate this to a diffusion process or Brownian motion where the diffusivity is equal to the variance of the variable $$\epsilon$$.

The time to reach a certain point, the first hitting time, follows an inverse Gaussian distribution. Or since there is no drift it is a Levy Distribution. Then the hitting time is distributed according to a Levy distribution with parameters $$m=0$$ and $$s = (10/\sqrt{3})^2$$. We can use the cumulative distribution function to model the fraction of cases that have hit the point of zero money after 1000 steps.

### Example

The graph and code below demonstrate the computation with the Markov chain and the estimation with the Levy distribution.

kmax <- 3000

### a kmax times 3kmax matrix for the
### probability to be with profit x in step k
###
### note: in R code the index starts with 1, and this relates to 0 money
###
Pxk <- matrix(rep(0,3*kmax^2),3*kmax)

### compute each coin toss
for (i in 2:kmax) {
### compute the cases when money is 5 or larger
for(j in 4:(3*kmax-5)) {
Pxk[j,i] <- (1/8)*Pxk[j-3,i-1] + (3/8)*Pxk[j-1,i-1] + (3/8)*Pxk[j+1,i-1] + (1/8)*Pxk[j+3,i-1]
}
### compute the special cases when money is 0,1,2 or 3 or smaller
Pxk[1,i] <- Pxk[1,i-1] +  (4/8)*Pxk[2,i-1] + (1/8)*Pxk[3,i-1] + (1/8)*Pxk[4,i-1]
Pxk[2,i] <- (3/8)*Pxk[3,i-1] + (1/8)*Pxk[5,i-1]
Pxk[3,i] <- (3/8)*Pxk[2,i-1] + (3/8)*Pxk[4,i-1] + (1/8)*Pxk[6,i-1]
Pxk[4,i] <- (3/8)*Pxk[3,i-1] + (3/8)*Pxk[5,i-1] + (1/8)*Pxk[7,i-1]
}

### plot the simulation
plot(Pxk[1,], type = "l",
ylab = expression(P[ruin]), xlab = "number of tosses", ylim = c(0,1))

### add a curve based on the Levy distribution
n <- c(1:3000)
dist <- 10
sigma <- 2 * 1/8 * 3^2 + 2 * 3/8 * 1^2 ### variance of steps relates to diffusion rate
lines(n,rmutil::plevy(n, m = 0, s = dist^2/sigma), col= 2)

### highlight the point for 1000 tosses
points(1000,Pxk[1,1000], pch = 21, col = 1, bg = 0)
text(1000,Pxk[1,1000], expression(P %~~% 0.85), pos = 1, cex = 0.7)

legend(1000,0.4, c("exact computation","Levy distribution (diffusion model)"),
col = c(1,2), lty = 1, cex = 0.7)


I think you are correct. I wrote the following simulation (Python 3) and got the same result as you (i.e. that the probability of ruin is ~0.84).

import random

def flip_3_coins():
return sum(random.choice([1,-1]) for _ in range(3))

num_ruined = 0
num_trials = 1000

for trial in range(num_trials):
equity = 10
for flip in range(1000):
equity += flip_3_coins()
if equity <= 0:
num_ruined += 1
break

print(num_ruined/num_trials)


The probability of ruin converges to 1 if you increase the number of flips (e.g. if you change this from 1000 to 10000 then the probability of ruin becomes ~0.95).