in my stochastic processes class, Prof Mike Steele assigned a homework problem to calculate the ruin probabilities for Playing a game where you have with $ probability p and lose $ with probability 1-p. The probability of winning is no specified, so it can be a biased game. Ruin probabilities is defined to being the probability that's a game you win ten before losing, win before Lo Sing, and win before losing, etc. In total, I found three distinct methods to calculate.
This was a particularly great example to illustrate what to solve a problem using three fundamentally different methods:the First is theoretical calculation, second are simulation to obtain asymptotic values, and third are numerical linear algebra (Matrix algorithm) which also gives exact values.
Method 1:first Step Analysis and Direct computation of ruin probabilities
Let h (x) be the probability of winning $n before losing stake of $x.
First step analysis gives us a system of three equations:h (0) = 0; H (n) = 1; H (x) = P*h (x+1) + (1-p) *h (x-1).
How to solve this system of equations? We need the "one" trick and the telescoping sequence.
The trick is: (P + (1-P)) * H (x) = h (x) = P*h (x+1) + (1-p) *h (x-1) = p* (H (x+1)-H (x)) = (1-p) * (H (X)-H (x-1)) + H ( x+1)-h (x) = (1-p)/p * (H (X)-H (x-1))
Denote h (1)-h (0) = C, which is unknown yet, we have a telescoping sequence:h (1)-h (0) = C; H (2)-h (1) = (1-p)/p * C; H (3)-H (2) = ((1-p)/p) ^2 * C ... h (n)-h (n-1) = ((1-p)/p) ^ (n-1) * C.
Now, add to the telescoping sequence and use the initial conditions, we get:1 = h (n) = c* (1+ ((1-p)/P) + ((1-p)/p) ^2 +: . + ((1-p)/p) ^ (n-1)) = + c = (1-(1-P)/P)/(1-((1-p)/P) ^n-1). So h (x) = c * (((1-p)/p) ^ x-1)/((1-P)/p)-1) = (((1-p)/p) ^ x-1)/((((1-p)/p) ^n-1)
Method 2:monte Carlo Simulation of Ruin probabilities
The idea was to simulate sample paths from initial stake of $x and stop when it either hits 0 or targeted wealth of N.
We can specify the number of trials and get the percentage of trials which eventually hit 0 and which Eventuallyhit N. Thi S is important-in fact, I think the essence of Monte Carlo method are to has a huge number of trials to maintain Accurac Y, and to get a percentage of the number of successful trials in the total number of trials.
In each step of a trial, we need a Bernoulli the random variable (as in a coin flip) to increment x by 1 with probability p an D-1 with probability 1-p.
In Python this becomes:
From NumPy import Randomimport numpy as Npdef MC (x,a,p): End_wealth = a init_wealth = x list = [] for k in range (0, 10 00000): While x!= end_wealth and x!= 0:if np.random.binomial (1,p,1) = = 1:x + 1 Else:x-= 1 if x = = A:list.append (1) else:list.append (0) x = init_wealth print float (sum (list))/len (list) MC (10,20, 0.4932) MC (25,50,0.4932) MC (50,100,0.4932)
You can see the result of this simulation by plugging in P = 0.4932 = (18/37) *.5 +. 5*.5 = 0.4932, which is the Probabilit Y of winning the European roulette with prisoner ' s rule. As the number of trials get bigger and bigger, the result gets closer and closer to the theoretical value calculated under Method 1.
Method 3:tridiagonal System
According to wiki, a tridiagonal system have the form of a_i * x_i-1 + b_i * x_i + c_i * x_i+1 = d_i where I ' s is indices.
It is clear that the ruin problem exactly satisfies this form, i.e. H (x): = probability of winning n starting from I, h (x) = (1-p) *h (x-1) + p*h (x+1) =-(1-p) *h (x-1) + H (x)-p*h (x+1) = 0, H (0) = 0, h (n) = 1.
and therefore, for the tridiagonal matrix, the main diagonal consists of 1 ' s, and the upper diagonal consists of-(1-p) ' s, And the lower diagonal consists of-p ' s.
In Python this becomes:
Import NumPy as Npfrom scipy import sparsefrom scipy.sparse.linalg Import spsolven = 100p = 0.4932q = 1-pd_main = Np.ones ( n+1) D_super =-P * d_maind_super[1] = 0d_sub =-Q * d_maind_sub[n-1] = 0data = [D_sub, d_main, d_super]print dataa = spars E.spdiags (data, [ -1,0,1], n+1, n+1, format= ' csc ') b = Np.zeros (n+1) B[n] = 1x = Spsolve (A, b) Print X
Gambler ' s ruin problem and 3 Solutions