In my stochastic processes class, Prof Mike Steele assigned a homework problem to calculate the ruin probabilities for playing a game where you with 1 dollar with probability p and lose 1 dollar with probability 1-p. The probability of winning is not specified, so it can be a biased game. Ruin probabilities are defined to be the probability that in a game you win 10 before losing 10, win 25 before losing 25, and win 50 before losing 50, etc. In total, I found three distinct methods to calculate.
This is a particularly great example to illustrate how to solve a problem using three fundamentally different methods: the first is theoretical calculation, second is simulation to obtain asymptotic values, and third is 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 dollars.
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 up 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 is to simulate sample paths from initial stake of x dollars 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. This is important - in fact, I think the essence of Monte Carlo method is to have a huge number of trials to maintain accuracy, 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 random variable (as in a coin flip) to increment x by 1 with probability p and -1 with probability 1-p.
In Python this becomes:
from numpy import random
import numpy as np def MC(x,a,p):
end_wealth = a
init_wealth = x
list = []
for k in range(0, 1000000):
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 probability of winning the European Roulette with *er'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 has the form of a_i * x_i-1 + b_i * x_i + c_i * x_i+1 = d_i where i's are 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 np
from scipy import sparse
from scipy.sparse.linalg import spsolve n = 100
p = 0.4932
q = 1-p d_main = np.ones(n+1)
d_super = -p * d_main
d_super[1] = 0
d_sub = -q * d_main
d_sub[n-1] = 0 data = [d_sub, d_main, d_super]
print data
A = sparse.spdiags(data, [-1,0,1], n+1, n+1, format='csc') b = np.zeros(n+1)
b[n] = 1
x = spsolve(A, b)
print x