Gamblers Ruin with multiple coins tossed at the same time

Question

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 
    added to equity.
    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)

Answer 1

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)

Pxk[11,1] = 1 ### start with x=10 money in step k=1

### 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)

Answer 2

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).