Generalized Boundary Problem for One-Dimensional Random Walks

by
David K. Neal
Department of Mathematics
Western Kentucky University
Bowling Green, KY  42101
nealdk@wkuvx1.wku.edu

Mathematica in Education
Vol. 3 No.4
Copyright 1994
TELOS/Springer-Verlag

Abstract

The boundary problem for non-symmetric one-dimensional random walks
is described and simulated graphically.  We estimate the probability of a 
particle reaching the top boundary before the lower boundary and estimate 
the average time to reach either boundary.  Confidence intervals for these
estimates are also given.  A matrix solution to these problems is then given
followed by an application to the Gambler's Ruin problem.

Introduction

A famous problem in probability theory is the “Drunkard's Walk” which can 
be described as follows.  Suppose a particle begins at height m > 0 on the y-axis.
On each unit time interval, the particle moves up 1 step with probability p or
down 1 step with probability q = 1-p.  (i) What is the probability that the particle 
reaches either height n > m or height 0?  (ii) What is the probability P[m] that 
the particle reaches height n before height 0?  (iii) What is the average number 
of steps T[m] needed to reach one of these two boundaries?

Each problem can be solved by one standard technique involving difference 
equations [Chu75, p. 243; Fel57, p. 313].  Most importantly, these solutions yield
closed-form formulas.  Namely, (i) the particle will reach a boundary with 
probability 1; (ii) for p = .5, P[m] = m/n, and for p ≠ .5, P[m] = (1-(q/p)m)/(1-(q/p)n),
and (iii) for p = .5, T[m] =  m(n-m), and for p ≠ .5, T[m] = 
n/(p-q)*(1-(q/p)m)/(1-(q/p)n) - m/(p-q).  These formulas can easily be 
generalized for any lower boundary less than m.

We can generalize the problem further by assuming that the particle moves 
up height h with probability p or down height k with probability q = 1-p, where 
h and k are any positive real numbers.  Unfortunately, the above difference 
equation technique will no longer work.  This problem has been considered in
connection with sequential sampling when h and k are integers [Fel57, p. 330];
however, there is not a known closed-form solution.  In this article, we shall 
describe a probabilistic method of approximating the solutions for given values
of p, m, n, h, and k.  Then we shall give an algebraic method for finding the
exact solution for any given set of initial values provided h and k are integers.
The techniques will be implemented with brief Mathematica  programs.

The material in this article is one of a dozen or so projects that the author
uses as supplemental material for a beginning calculus-based course in
probability and statistics.  For all such projects, students are to read the problem
and solution and then run the program.  Students then work a few exercises
which involve slight modifications of the program.  These students generally
have had no prior experience with Mathematica  or Macintosh machines.  By 
the end of the semester however, they seem to have a good working knowledge
of both.  Moreover, they gain a much greater understanding of probabilistic
concepts and have fun doing it.

Probability Estimation Technique

To approximate the probability of reaching n first, we can apply the Law
of Large Numbers to a sample of random walks and merely count the 
proportion of particles which reach or exceed the top boundary first.  We can
also average the number of steps needed to reach or exceed a boundary for this
sample in order to approximate the true average time.  We can then establish 
bounds for our approximations by constructing confidence intervals. 
Mathematica  can easily generate random walks which fit the parameters of the 
problem and can also display graphs which illustrate the problem.

We shall demonstrate with the following values : p = .40, h = 3.2, k = 2.4, 
m = 6,  and n = 10.  We first enter all the variables:

p = .40;
q = 1-p;
StepsUp = h = 3.2;
StepsDown = k = 2.4;
StartingHeight = m = 6;
UpperBoundary = n = 10;

We shall make our estimations based on, say 300, trials.  Therefore, 
we must construct a loop of 300 random walks which begin at starting height 
m and go up h steps with probability p or down k steps with probability q.  Each 
walk stops when reaching the upper boundary of 10 or the lower boundary of 0.  
(If desired, one can always change the lower boundary to any value below m.)
The height of the ith particle on its jth step is denoted by x[i,j], for i = 1,...,300.
For the program, we first set the number of trials and set all the initial heights 
equal to m.  Within the loop, t[i] counts the number of steps that the ith walk 
needed to reach a boundary and u[i] determines whether or not the walk ended 
at or above height n.

trials = 300;
Do[x[i,0] = m, 
  {i, 1, trials}];

Do[j=0;
   While[0 < x[i,j] < n,
        x[i,j+1] = If[Random[] < p, 
                     x[i,j] + h, 
                     x[i,j] - k];
        j++];
   t[i] = j;
   u[i] = If[x[i,t[i]] >= n,
            1, 0],
{i, 1, trials}]

In order to see graphs of the generated random walks, where the x-axis 
denotes time and the y-axis denotes height, we can use the command below. 
Here we display only the the first 2 graphs; but the index limits can be changed, 
for instance to {i,10,20}, to display any consecutive sequence of graphs.

Do[ListPlot[Table[{j, x[i,j]},
                 {j, 0, t[i]}],
           PlotJoined -> True,
           AxesOrigin -> Automatic],
  {i, 1, 2}]

Lastly, we count the proportion to reach the top boundary first and the 
average time to reach a boundary for this sample:

p1 = ProportiontoHitUpper =
   N[Sum[u[i], {i,1,trials}] / trials]

SampleMeanTime =
   N[Sum[t[i], {i,1,trials}] / trials]

When using h = k = 1, we can observe the accuracy of the approximations
 by adding codes which compute the known real values :

RealProb=If[p==.5,m/n,(1-(q/p)^m)/(1-(q/p)^n)]

AveTime=If[p==.5,m(n-m),
n/(p-q)*(1-(q/p)^m)/(1-(q/p)^n) - m/(p-q)]

Confidence Intervals

We can now use the data from this sample of random walks to create a
confidence interval for the probability P[m] of reaching n first when starting
at height m.  The confidence interval for P[m] is centered around the sample 
proportion p1.  Since this measurement is a proportion, the true variance is 
known to be P[m](1-P[m])  [HoTa93, p.151], which has a maximum value of 
.25 since 0 ≤ P[m] ≤ 1.  By taking the square root, we obtain a maximum standard 
deviationof .50.  The variance of the sample proportion however equals the 
true variance divided by sample size [HoTa93, p. 267]; hence, the confidence 
interval is centered at p1 with standard deviation .5/Sqrt[trials].

Using the Statistics`ConfidenceIntervals`  package, we now compute a 95%
confidence interval for P[m] :

<<Statistics`ConfidenceIntervals`

c=.95;

ProbCI = 
   N[NormalCI[p1,.5/Sqrt[trials],
             ConfidenceLevel->c]]

Similarly, we can construct a confidence interval, centered around the 
sample mean xbar, for the average time µ to reach a boundary.  Though we 
would now use the unbiased sample deviation sd as an estimate for the 
unknown standard deviation s of the time to reach a boundary.  However, 
there are some other factors to consider before we do so.

Unlike with the sample proportion, we do not have an upper bound for
the variance of the time it takes to reach a boundary, nor do we know if this
variance is even finite.  If the variance is finite, then we may apply the Central
Limit Theorem to construct a confidence interval based on a normal distribution 
about the sample mean.  However, we must be sure that the sample size trials 
is large enough so that sd is a good estimate of s and so that the the standard 
normal distribution is a good approximation of the distribution of 
(xbar - µ ) / (sd/Sqrt[trials]) .

To test these assumptions, it is a better practice to run our simulations in
batches, say 30 batches of 100 trials each.  We then consider the 30 batch means
as a new random sample.  If these 30 measurements seem normally 
distributed, then it is good evidence that the confidence interval constructed 
will be valid.  The commands below create and display such batch means.

Do[Do[j=0;While[0<x[i,j]<n,
     x[i,j+1]=If[Random[]<p,x[i,j]+h,x[i,j]-k];j++];
     t[i]=j,{i,1,100}];
     b[s]=N[Sum[t[i],{i,1,100}]/100],{s,1,30}]

TimeData=Table[b[s],{s,1,30}]

One can observe that these sample means do seem to be consistent with 
a rather small variance.  As an exercise, one can perform a goodness of fit test
for normality.  Finally, we construct a 95% confidence interval using the
sample mean and sample deviation of the above batch means :

mu=N[Mean[TimeData]];
sd=N[StandardDeviation[TimeData]];

TimeCI=N[NormalCI[mu,sd/Sqrt[30],
         ConfidenceLevel->c]]

Matrix Solution

We now consider the special case when the upward and downward steps h 
and k are integers.  We shall describe an algebraic solution for the probability 
P[m] of reaching height n before height 0 when starting at height m, for 
0 ≤ m ≤ n.  We first note that P[n] = 1, since if we start at n, then we reach n
before 0 with probability 1.  Thus, we also define P[k] = 1 for k > n.  Moreover, 
P[0] = 0, since we cannot reach n before 0 if we start at 0, and we also define 
P[k] = 0 for k < 0.  So, we must solve for P[1] through P[n-1].  We shall do so by 
solving n-1 equations having these n-1 unknowns.  To avoid consideration of 
numerous cases, we assume that k+h < n-1.  Simple modifications can be made
if k+h ≥ n - 1.

When starting at height m, for 1 ≤ m ≤ n-1, after one step we are at height 
m-k with probability q or at height m+h with probability p.  Then we must 
reach n before 0 when starting  at m-k or when starting at m+h.  Hence,

P[m] = q P[m-k] + p P[m+h].

In particular, the first k equations are

P[1] = q P[1-k] + p P[1+h] = p P[1+h]
                              
                              P[2] = q P[2-k] + p P[2+h] = p P[2+h]
                              
                                                . . . . . . . . . . . . . . .
                                                
                                                . . . . . . . . . . . . . . . 
                                               
                               P[k] = q P[k-k] + p P[n-1] = p P[k+h] .

Since k+h < n-1,  the next n-1-h-k equations are

P[k+1] = q P[1] + p P[k+1+h] 
             
                                                  . . . . . . . . . . . . .
             
                                                  . . . . . . . . . . . . .
                                                  
                                 P[n-1-h] = q P[n-1-h-k] + p P[n-1],

and the last h equations are

P[n-h] = q P[n-h-k] + p P[n] = q P[n-h-k] + p
      
                                                . . . . . . . . . . . . . . .
             
                                                . . . . . . . . . . . . . . .
             
                     P[n-1] = q P[n-1-k] + p P[n-1+h] = q P[n-1-k] + p

The first k equations can be written as follows:

P[1] - p P[1+h] = 0
                
                                            P[2] - p P[2+h] = 0
                
                                                       . . . .
                   
                                            P[k] - p P[k+h] = 0

The next n-1-h-k equations can be written as follows:

-q P[1] + P[k+1] - p P[k+1+h] = 0
                                
                                                  . . . . . . . . . .
                        
                            -q P[n-1-h-k] + P[n-1-h] - p P[n-1] = 0 .

Finally, the last h equations can be written as

- q P[n-h-k] +  P[n-h] = p
                                   
                                                    . . . . . .
                    
                                    - q P[n-1-k] + P[n-1] = p .

Hence, we obtain the following (n-1)x n augmented matrix which represents 
the n-1 equations, having P[1] through P[n-1] as the n-1 unknowns, where the 
-p is h columns to the right of the 1 in each of the first n-1-h rows, the -q is k
columns to the left of the 1 in the last n-1-k rows, and p occurs in the right-
most column in the last h rows.

1  0  0  0  . -p  0  .  .  .  .  0  0  0 | 0
       0  1  0  0  .  0 -p  0  .  .  .  0  0  0 | 0
       0  0  1  0  .  0  0 -p  0  .  .  0  0  0 | 0
       0  0  0  1  .  0  0  0 -p  0  .  0  0  0 | 0
       .  .  .  .  .  .  .  .  .  .  .  .  .  . | 0
      -q  0  0  .  .  1  0  0  0  . -p  .  0  0 | 0
       0 -q  0  0  .  .  1  0  0  0  . -p  0  0 | 0
       0  . -q  0  0  .  .  1  0  0  0  . -p  0 | 0
       .  .  .  .  .  .  .  .  .  .  .  .  .  . | 0
       0  .  .    -q  0  0  .  .  1  .  .  .  0 | p
       0  .  .  .    -q  0  0  .  .  1  .  .  0 | p
       0  .  .  .  .  . -q  0  0  .     1  .  0 | p
       0  .  .  .  .  .  . -q  0  0  .  .  1  . | p
       0  .  .  .  .  .  .    -q  0  0  .  .  1 | p

If we let A denote the (n-1)x(n-1) matrix of coefficients on the left side of 
the augmented matrix, and let B denote the right-most column, then the  
system can be represented by the equation AX = B, where X denotes the  
column of unknowns P[1],. . ., P[n-1].  The solution for X is then given by  
X = A-1 B.

Of course it is not practical to solve the system of equations by hand; thus,        
we can solve the system numerically for specific p and q or theoretically for 
arbitrary p and q by using a short program.

Again, we first enter the necessary parameters and we note that this 
technique will provide a solution for all starting heights m; so, we can clear 
this variable.  Moreover, we have changed the values of h and k to integers.

Clear[m,a,b];
p = .40;
q = 1-p;
StepsUp = h = 3;
StepsDown = k = 2;
UpperBoundary = n = 10;

The following commands create the (n-1)x(n-1) matrix of coefficients A:

a=IdentityMatrix[n-1];
Do[a[[i,i-k]]=-q,{i,k+1,n-1}];
Do[a[[i,i+h]]=-p,{i,1,n-1-h}];

a

The next command creates the (n-1)x 1 matrix B :

b=Table[If[i<=n-h-1,0,p],{i,1,n-1},{j,1,1}]

Lastly, we solve the system AX = B for X using the LinearSolve command.
We then display a table of values {m,P[m]}, for 1 ≤ m ≤ n-1, which gives the
probabilities of reaching n first when starting at height m.  We have displayed 
the output for the last command which occurs using our initial values.

Soln = LinearSolve[a, b];

Table[{m, Soln[[m]]}, {m,1,n-1}]

Other Matrix Solutions

We now let R[m] be the probability of reaching (or exceeding) a boundary of 0
or n when starting at height m.  With only one modification, we can use the 
same matrix technique as above to verify that R[m] = 1 for all m.  We note that
R[0] = R[n] = 1, since if we start at 0 or n then we reach 0 or n.  Likewise R[k] = 1 
for k < 0  and k > n.  The only modification then is that the first k equations can 
be written with the right hand side equaling q rather than 0.  We need only 
create a new right-most column for the augmented matrix and then solve for 
R[1],. . .,R[m] :

b2 = Table[If[i<=k,q,If[i<=n-h-1,0,p]],
          {i,1,n-1}, {j,1,1}]

Soln2 = LinearSolve[a, b2]

With another modification, we can find the average number of steps needed 
to reach a boundary.  If we let T[m] denote the average number when starting 
at height m, then T[0] = T[n] = 0 since it does not require a step to reach the
boundary if we start at a boundary.  Likewise, T[k] = 0 for k < 0 and k > n.  For
0 < m < n however, we take 1 step and then we proceed from height m-k with 
probability q or from height m+h with probability p.  Thus,
T[m] = 1 + q T[m-k] + p T[m+h].  Since, T[k] = 0 for k ≤ 0 and k ≥ n, we will obtain 
the same system of equations as before; but now, each entry in the right-most
column will be 1.  Again we need only create a new column and then solve
for each T[m].

b3=Table[1,{i,1,n-1},{j,1,1}]

Soln3 = LinearSolve[a, b3];

Table[{m,Soln3[[m]]},
     {m, 1, n-1}]

For h = k = 1, we could also verify the matrix solutions for P[m] and  
T[m] by  adding commands which list the previously known formulas:

Table[{m,If[p==.5,m/n,(1-(q/p)^m)/(1-(q/p)^n)]},
{m,1,n-1}]

Table[{m,If[p==.5,m(n-m),
n/(p-q)*(1-(q/p)^m)/(1-(q/p)^n) - m/(p-q)]},{m,1,n-1}]

Application To Gambling

A famous application of the “Drunkard's Walk” is the “Gambler's Ruin” 
problem.  We first find the probability Q[m] = 1 - P[m] of reaching 0 before 
reaching n.  This setting can now be interpreted as the probability of “going 
broke” before reaching a goal of $n if starting with $m when one either wins 
$h with probability p or loses $k with probability q on each bet.  For h = k = 1, 
we of course have a closed-form expression: For p = .5, Q[m] = (n-m)/n and for
p ≠ .5, Q[m] = ((q/p)^m - (q/p)^n)/(1 - (q/p)^n).  Then for p ≤ q (or
equivalently p ≤ .50),  lim Q[m] = 1 for all m.
                                       n->∞

This limit is the famous Gambler's Ruin: When one's chance of winning 
on a single bet is no more than .50, then the persistent gambler will eventually 
go broke with probability 1.

In real gambling situations though, one usually expects to win more than 
$1 for each dollar bet unless the probability p of winning is close to .50.  For a $1
bet with a payoff of $h, one's average earnings are -1*q + h*p.  For a fair game, 
the average earnings should be 0; hence, h = q/p.  Since most games are in fact 
unfair, so that the house can make a profit, we shall assume a payoff of $h 
which equals the next integer below q/p.  For p = q though, we should assume
a fair game with h = 1.

For example, if p = .20 and q = .80, then q/p = 4 and we take h = 3.  We then
let k = 1 in the original matrix program.  In other words, we either win $3 with                                       
probability .20 or lose $1 with probability .80.  Now suppose we wish to reach
$30 before going broke with probability .60.  Is  it possible?  If so how much 
money do we initially need?  By running the program with p = .20, h = 3, k = 1, 
and n = 30, we obtain the following final output :

Hence, if we start with $27, we have just slightly more than a 60% chance 
of reaching $30 before going broke.  By experimenting with larger and larger 
values of n, it becomes evident that there is an upper bound to the probability 
of reaching $n before going broke.  With the above figures, this upper bound 
appears to be .75.  Moreover, we cannot contradict the Gambler's Ruin. That is,
for a fixed initial amount m, if n-> ∞, then P[m]-> 0 and hence Q[m] -> 1.  One
will still eventually go broke if one persists in gambling.  However, if one 
decides to “quit while you're ahead”, one can determine the precise probability 
of reaching the goal of $n before going broke!

Exercise

Find an estimation and confidence interval for the average “area under the
curve” A[m] when starting at height m and stopping upon reaching height n
or height 0.  Adjust the matrix solution to find the exact value of A[m] when
the steps h and k are integers.

Acknowledgments

The author would like to thank the referees for their suggestions and 
pointers which made the programming more efficient and elegant.

References

[Chu75]   K.L. Chung. Elementary Probability Theory with Stochastic 
Processes . Springer-Verlag, New York, 1975.

[Fel57]  W. Feller.  An Introduction to Probability Theory and Its Applications 
Volume I, Second Edition.  John Wiley & Sons, Tokyo, 1957.

[HoTa93]  R.V. Hogg and E.A. Tanis. Probability and Statistical Inference  
Fourth Edition.  Macmillan Publishing Co., New York, 1993.

About the Author

David K. Neal is an associate professor of mathematics at Western Kentucky
University where he has taught for six years.  His mathematical interests are
in stochastic analysis and general probability.  In recent years, he has enjoyed
enriching his courses with a variety computer projects using Mathematica .
His other hobbies include music, sports, and felines.

David K. Neal
Department of Mathematics
Western Kentucky University
Bowling Green, KY  42101
nealdk @wkuvx1.wku.edu