Gambler’s Ruin and First Passage Time

Jan Vrbik

We investigate the classical problem of a gambler repeatedly betting $1 on the flip of a potentially biased coin until he either loses all his money or wins the money of his opponent. This is then extended to the case of his adversary (a casino) having practically unlimited resources and used to derive the inverse Gaussian distribution of the first passage time.

Probability of Winning the Game

We assume that the probability of winning a single round of this game is p, which is usually close, but not necessarily equal, to (1/2). Thus, for example, when betting $1 on black in the game of roulette, with 18 black, 18 red, and 1 green, all with equally likely outcomes, p is equal to (18/37)=0.4865. We also assume that our player and his adversary start with a dollars and b dollars, respectively, and that they are both determined to continue playing until one of them goes broke.

Our first objective is to find the probability that our player wins the game, ending up with a+b dollars. To do this, we have to imagine that the game has been going on for some time, and the player has reached the point of having exactly i dollars in his pocket, so that his opponent has a+b-i. Given that, we denote our player’s probability of winning the game by (Subscript[ω, i]). If one more round is played, one can see that

(Subscript[ω, i])=p (Subscript[ω, i+1])+(1-p)(Subscript[ω, i-1]),

depending on whether the player wins (the first term) or loses (the second one) the next round []. The end conditions are (Subscript[ω, a+b])=1 (the game is already won) and (Subscript[ω, 0])=0 (the game is lost). We thus have a+b-1 linear equations for the corresponding (Subscript[ω, i]) probabilities, where 0<i<a+b. The equations are of a rather special type, as each of them involves only three consecutive values of (Subscript[ω, i]), so rather than using the usual LinearSolve procedure, it is more appropriate to switch to RSolve.

sol=RSolve[{ω[i]==p ω[i+1]+(1-p)ω[i-1],ω[0]==0,ω[a+b]==1},ω[i],i][[1]]

When p=(1/2), the solution is indeterminate, but we get the correct answer in the limit.

Limit[ω[i]/.sol,p->(1/2)]

All we have to do now is to evaluate (Subscript[ω, i]) using i=a. To this end, we build the following function.

Pw[a_,b_,p_]:=(1-((((1-p)/p))^a))/(1-((1-p)/p)^(a+b))

Pw[a_,b_,1/2]:=a/(a+b)

We apply the function to the cases discussed so far.

Pw[10,10,1/2]//N
Pw[10,10,18/37]//N
Pw[100,100,18/37]//N

From these examples we can see that when both players start with the same amount of money and the coin is fair, each of them can win with the same probability of 50%. But as soon as the coin is slightly biased, the probability of the disadvantaged player winning decreases; this disparity becomes more pronounced as the initial capital of both players increases.

In practice, things are even worse than that: the casino starts with practically unlimited resources, which implies that a disadvantaged player does not stand any chance of winning whenever p≤(1/2). Should p be slightly higher than (1/2) though, his chances to continue winning indefinitely are computed by taking the limit of (Subscript[ω, a]) as b->∞, resulting in

1-(((((1-p)/p)))^a).

Distribution of the Game’s Duration

To find the distribution of the random number of rounds it takes to complete a game with two players, each with finite resources, given that currently the first player has i dollars in his pocket, we introduce the corresponding probability generating function (Subscript[H, i])(z). Similarly to how we solved for (Subscript[ω, i]), we can now set up the equations

(Subscript[H, i])(z)=z (p (Subscript[H, i+1])(z)+(1-p)(Subscript[H, i-1])(z)),

where the multiplication by z is necessary to account for the extra round played. We also know that (Subscript[H, 0])(z)=(Subscript[H, a+b])(z)=1 (a probability generating function of 1 indicates that there are no more rounds to be played).

We now solve these equations.

RSolve[{H[i]==z p H[i+1]+z(1-p)H[i-1],H[0]==1,H[a+b]==1},H[i],i][[1]]//FullSimplify

Substituting a for i yields the corresponding result.

H[a_,b_,p_]:=Module[{r=((1+(Sqrt[1-4p(1-p)(z^2)]))/(2p z)),s=((1-(Sqrt[1-4p(1-p)(z^2)]))/(2p z))},(((r^a)(1-(s^(a+b)))-(s^a)(1-(r^(a+b))))/((r^(a+b))-(s^(a+b))))]

Using this function, we can now easily find and display the corresponding distribution of the game’s duration.

aux=Series[N[H[10,10,18/37],36],{z,0,400}]; 
ListPlot[Table[{i,Coefficient[aux,z,i]},{i,10,400,2}],AspectRatio->.3,ImageSize->700]

One can see that these games can easily last hundreds of rounds (things get worse when p approaches (1/2)).

Sometimes it is sufficient to know only the mean and variance of this distribution.

μ[a_,b_,p_]:=((a-(a+b)Pw[a,b,p])/(1-2p))

Var[a_,b_,p_]:=((a+b)/(1-2p))(((4(((((1-p)/p)))^a)μ[a,b,p])/((((((1-p)/p)))^(a+b))-1))-((3(a+b)Pw[a,b,p](1-Pw[a,b,p]))/(1-2p)))+((4p(1-p)μ[a,b,p])/(((1-2p))^2))

The correctness of these formulas can be verified by recalling that the mean μ of a distribution is given by H'(z) evaluated at z=1, and its variance similarly by H''(z=1)+μ-(μ^2).

(D[H[a,b,p],z]/.z->1/.(((1-4(1-p)p))^(q_)):>(((1-2p))^(2q)))-μ[a,b,p]//Simplify

(D[H[a,b,p],{z,2}]/.z->1/.(((1-4(1-p)p))^(q_)):>(((1-2p))^(2q)))+μ[a,b,p]-((μ[a,b,p])^2)-Var[a,b,p]//Simplify

The two formulas can be extended to the case of p=(1/2) by taking the corresponding limit.

Limit[μ[a,b,p],p->1/2]

Limit[Var[a,b,p],p->1/2]

These impose the following.

μ[a_,b_,1/2]:=a b

Var[a_,b_,1/2]:=((a b((a^2)+(b^2)-2))/3)

The Case of an Infinitely Rich Adversary

By analyzing the local variables r and s in the definition of H, we can see that r>1 and 0<s<1, for any 0≤z<1. As b->∞, we thus get

(Subscript[H, a])(z)=(s^2)=(((((1-(Sqrt[1-4p(1-p)(z^2)]))/(2p z))))^a).

When p>(1/2), the probability that the game finishes in finite time is (((((1-p)/p)))^a), according to (). The conditional probability generating function of the total number of rounds, given that the game does not continue indefinitely, is thus

(((((1-(Sqrt[1-4p(1-p)(z^2)]))/(2(1-p) z))))^a).

The more interesting case is when p≤(1/2), which implies that the game cannot continue indefinitely. Expanding () in powers of z would enable us to plot the corresponding distribution, as was done at the beginning of this section; we leave this as an exercise. Here is the mean and variance of the total number of rounds.

aux=(((((1-(Sqrt[1-4p(1-p)(z^2)]))/(2p z))))^a);

D[aux,z]/.z->1/.(((1-4(1-p)p))^(q_)):>(((1-2p))^(2q))//Simplify

(D[aux,{z,2}]/.z->1/.(((1-4(1-p)p))^(q_)):>(((1-2p))^(2q)))+%-((%)^2)//Simplify

Brownian Motion

In this last section, we explore yet another interesting limit of the gambler’s ruin problem. First we assume that the game is happening in “real time” t and that n rounds are played during each hour (or any other unit of time). We also assume that

p=(1/2)+(d/(2(Sqrt[c·n]))),

where d and c>0 are two constants, and that (Sqrt[(c/n)]) dollars are bet in each round instead of the original $1. In the limit as n->∞, this results in a new process called Brownian motion that runs in real (i.e. continuous) time and whose values constitute a continuous, even though nowhere differentiable, function of t [].

To visualize its possible course, we take d=-0.35, c=2, and n=20 and display a random realization of the process; that is, we display its current values—the amount of money in the gambler’s pocket—as a function of t during the next 50 hours, taking the initial value to be $100.

Module[{n=20,c=2,d=-0.35,t=50,A=100,aux},
aux=(2RandomVariate[BinomialDistribution[1,(1/2)+(d/(2(Sqrt[c n])))],t n]-1)(Sqrt[(c/n)]);
ListPlot[Table[{i,A+Apply[Plus,Take[aux,i]]},{i,t n}],AspectRatio->.3,ImageSize->700]]

Assuming the gambler can borrow money and continue playing even when the value of the process (his current worth) is negative, it is easy to see that this value at time t is a normally distributed random variable with mean A+d t and variance c t, where A is the initial value, and d and c (introduced earlier) are the so-called drift and diffusion coefficients, respectively. This assertion can be proved by realizing that the moment-generating function of the net win (equal to 1 with probability p and to -1 with probability 1-p) in one round of the original game is M(u)=p (e^u)+(1-p)(e^(-u)). To adjust this to the new game (winning or losing (Sqrt[(c/n)]) dollars, instead of $1), replace each u of the last expression by u(Sqrt[(c/n)]). To get the moment-generating function of the total win (or loss, when negative) accumulated during a time interval of length t, we have to add the results of n t independent rounds, which corresponds to raising M(u) to the power of n t. Finally, we need to take the limit as n->∞.

Limit[(((((1/2)+(d/(2(Sqrt[c n]))))Exp[u(Sqrt[(c/n)])]+((1/2)-(d/(2(Sqrt[c n]))))Exp[-u(Sqrt[(c/n)])]))^(n t)),n->∞]

This is the moment-generating function of a normal distribution with expected value d t (the exponent’s coefficient of u) and variance of c t (the coefficient of ((u^2)/2)). Since it represents only the total net winnings up to time t, we have to add the initial value of A to get the distribution of the gambler’s net worth at time t; this will increase the expected value to A+d t.

Inverse Gaussian Distribution

We now establish the distribution of T, the first passage time through the value of zero, which is the time at which the gambler has to stop playing, since he has lost all his money. We assume that he starts with the amount A and that d<0; reaching the value of zero in finite time is thus guaranteed.

We proceed similarly to deriving the distribution of his net worth: we first find the moment-generating function of T in terms of the old game, with n rounds played every hour, and then take the limit as n->∞.

To get the moment-generating function of the number of rounds of the old game needed to go broke, we simply replace z in () by (e^u), and a by A(Sqrt[(n/c)]), the initial amount expressed in the new monetary units. To convert the result into hours, each consisting of n rounds, we need to replace u by (u/n).

Limit[(((((1-(Sqrt[1-(1-((d^2)/(c n)))Exp[2u/n]]))/((1+(d/(Sqrt[c n]))) Exp[u/n]))))^(A(Sqrt[(n/c)]))),n->∞]

The resulting moment-generating function can be rewritten in the following simple form, using only two parameters: Λ=-((d A)/c) and (Subscript[t, 0])=-(A/d) (since d<0, both are positive):

(e^(Λ(1-(Sqrt[1-2((Subscript[t, 0])/Λ)u])))).

The corresponding probability density function is

f(t)=(Sqrt[((Λ·(Subscript[t, 0]))/(2π (t^3)))])exp(-(((Λ(t-(Subscript[t, 0])))^2)/(2t (Subscript[t, 0]))))

for t>0 and f(t)=0 otherwise. To prove this, we compute the moment-generating function of (), which agrees with ().

Integrate[(Sqrt[((Λ t0)/(2π (t^3)))])Exp[-((Λ(((t-t0))^2))/(2t t0))+t u],{t,0,∞},Assumptions->{Re[t0 Λ]>0,2Re[u]<Re[Λ/t0]}]//PowerExpand

The corresponding distribution is called the inverse Gaussian; it has mean (Subscript[t, 0]) and variance (((Subscript[t, 0])^2)/Λ).

Module[{M=Exp[Λ(1-(Sqrt[1-2t0 u/Λ]))],m},m=D[M,u]/.u->0;{m,(D[M,{u,2}]/.u->0)-(m^2)}]

Its distribution function of cumulative probabilities can be expressed in terms of the distribution function of the standard normal distribution (denoted Φ) as:

Φ((Sqrt[(Λ/(t·(Subscript[t, 0])))])·(t-(Subscript[t, 0])))+(e^(2Λ))·Φ(-(Sqrt[(Λ/(t·(Subscript[t, 0])))])·(t+(Subscript[t, 0]))).

We verify this, getting an expression that agrees with ().

D[CDF[NormalDistribution[0,1],(Sqrt[(Λ/(t t0))])(t-t0)]
+Exp[2Λ]CDF[NormalDistribution[0,1],-(Sqrt[(Λ/(t t0))])(t+t0)],t]//Simplify

Applying the inverse Gaussian distribution to our previous example (c=2, d=-0.35, and A=100), we now display the distribution of time it takes the corresponding Brownian motion to reach zero for the first time.

Module[{c=2,d=-0.35,A=100,Λ,t0},Λ=-d A/c;t0=-A/d;
Plot[(Sqrt[((Λ t0)/(2π (t^3)))])Exp[-((Λ(((t-t0))^2))/(2t t0))],{t,0,600}]]

One last issue: to find the probability density function of T when p=(1/2), we need to take the limit of () as d->0.

Limit[(Sqrt[((A^2)/(2π c (t^3)))])Exp[-((d A/c(((t+A/d))^2))/(2t A/d))],d->0]//PowerExpand

This distribution has an infinite mean and variance.

Conclusion

The study of first passage times and their distributions is an important area of research. The author hopes that this article has provided a useful introduction to the topic.

References

W. Feller, An Introduction to Probability Theory and Its Applications, Vol. I, 3rd ed., New York: Wiley, 1968.

M. S. Bartlett, An Introduction to Stochastic Processes, with Special Reference to Methods and Applications, 2nd ed., Cambridge: Cambridge University Press, 1966.