Evaluating U(l) with Applied Linear Algebra: A Gambler's Demise

[Schwarzschild90]

Homework Statement

I ask for help in solving the exercises in this project on applied linear algebra. The problem outlined in the project is one in which we are tasked with modeling the demise of a gambler.

I need help solving exercise 1 (in red) on page 6. I have pasted the exercise text into the text body below, in block quotes.

Using the equations discussed above, evaluate numerically U(l) as a function of the initial gambler position l. Use e.g. M = 100 or similar. This is the code that I have written

M = 100;

X = ones(M,1); A = diag(-2*X,0); clear X X = ones(M-1,1); B = diag(X,1); C = diag(X,-1); T = A+B+C;

tau = 1; N = -(inv(T)/tau)*n0 n0 = eye;

numeric::int(N, l = 0..infinity)

Homework Equations

The Attempt at a Solution

 

[SteamKing]

There's no attachment or image pasted into the OP. The details of the problem are not clear.

 

[Schwarzschild90]

There's no attachment or image pasted into the OP. The details of the problem are not clear.

Thanks for pointing it out!

https://onedrive.live.com/redir?resid=6AFDD76B502B7B9C!4262&authkey=!AM2th-5i26TrCcI&ithint=file%2Cpdf

 

[Schwarzschild90]

Progress (and bump)

clear
M=10;
X = ones(M,1);
A = diag(-2*X,0);

clear X
X = ones(M-1,1);
B = diag(X,1);
C = diag(X,-1);
T = A+B+C;
clear X B C; % remove these matrices from memory

n0 = randi([1 M-1],1); % Initial gambler position, determined using a random seed generator

tau = 365

N = -inv(T)/tau * n0 % The ultimate probability of Gambler's ruin. R(l) = N(l)

 

[Schwarzschild90]

I can't make sense of what I must do to solve the exercise

 

[Ray Vickson]

I can't make sense of what I must do to solve the exercise

I can't make sense of the gambler's ruin problem as presented in the pdf file you attached. I know the gambler's ruin problem, and none f the versions I have seen have involved continuous time systems and differential equations. It seems like the author of the pdf wants you to treat the time between gambles as exponentially-distributed random variables, so you get a differential equation for the probability of "position" as a function of time. However, that is not involved with computing the probability of ruin.

 

[haruspex]

treat the time between gambles as exponentially-distributed random variables

No, the model is small gambles at very short equal intervals, making it effectively a continuous process.

 

[Schwarzschild90]

Yes thay's correct.

Should I be getting a greater than 1 chance of gambler's ruin occurring as a function of initial position?

I'm talking about this exercise:

• Using the equations discussed above, evaluate numerically U(l) as a function of the initial gambler position l. Use e.g. M = 100 or similar

My code
clear
M=4;
X = ones(M,1);
A = diag(-2*X,0);
clear X
X = ones(M-1,1);
B = diag(X,1);
C = diag(X,-1);
T = A+B+C;

tau=1; %For simplicity, we take tau = 1
n0=zeros(1,M)';
n0(1)=1; % Initial position is taken to be 1
ulP=-(inv(T)/tau)*n0;
plot(ulP,'r*');

 

[Ray Vickson]

Yes thay's correct.

Should I be getting a greater than 1 chance of gambler's ruin occurring as a function of initial position?

I'm talking about this exercise:

• Using the equations discussed above, evaluate numerically U(l) as a function of the initial gambler position l. Use e.g. M = 100 or similar

My code
clear
M=4;
X = ones(M,1);
A = diag(-2*X,0);
clear X
X = ones(M-1,1);
B = diag(X,1);
C = diag(X,-1);
T = A+B+C;

tau=1; %For simplicity, we take tau = 1
n0=zeros(1,M)';
n0(1)=1; % Initial position is taken to be 1
ulP=-(inv(T)/tau)*n0;
plot(ulP,'r*');

No, you should never, ever, get a probability > 1 for anything, ruin or otherwise. Sometimes, though, when you use approximate solution techniques and finite-wordlength floating-point computations you can have probabilities that exceed 1 slightly.

I do not understand your code, so have no way of assessing what you have done and where the errors (if any) arise. You do not even tell use what software you are using!

 

[Schwarzschild90]

No, you should never, ever, get a probability > 1 for anything, ruin or otherwise. Sometimes, though, when you use approximate solution techniques and finite-wordlength floating-point computations you can have probabilities that exceed 1 slightly.

I do not understand your code, so have no way of assessing what you have done and where the errors (if any) arise. You do not even tell use what software you are using!

I should mention that I use matlab.

 

[haruspex]

For simplicity, we take tau = 1

Is that valid? The file you attached does not explain what tau is, but it smells like a normalisation factor, i.e. its function is to ensure the probabilities of all the mutually exclusive events add up to 1.

 

[Schwarzschild90]

Is that valid? The file you attached does not explain what tau is, but it smells like a normalisation factor, i.e. its function is to ensure the probabilities of all the mutually exclusive events add up to 1.

Tau determines the rate at which ultimate ruin is reached. It can be of any magnitude, i.e. seconds, months or years.

 

[haruspex]

Tau determines the rate at which ultimate ruin is reached. It can be of any magnitude, i.e. seconds, months or years.

Then shouldn't its value affect the elements of T?
As a check, try varying tau. If I read your code correctly, the calculated probability will vary inversely. That would clearly indicate a problem.

 

[Schwarzschild90]

Varying tau varies the calculated probability inversely.

 

[haruspex]

Varying tau varies the calculated probability inversely.

Right, so you cannot arbitrarily plug in tau=1. I would think there is a relationship between tau and T.

 

[Schwarzschild90]

The project description doesn't mention anything about the relationship. By the way, I was informed that tau is a purely mathematical construct, which governs the rate at which gambler's ruin occurs. Thus, the value of tau can be completely arbitrary.

 

[haruspex]

The project description doesn't mention anything about the relationship. By the way, I was informed that tau is a purely mathematical construct, which governs the rate at which gambler's ruin occurs. Thus, the value of tau can be completely arbitrary.

But in that case the matrix T must depend on its value, no?
One approach is to sum the probabilities that should add to 1 then set tau in such a way that they do.

 

[Schwarzschild90]

Agreed.

The professor suggest that I set the value of tau equal to 1.

With tau = 1, I get the following column matrix:

ulP =

0.8000
0.6000
0.4000
0.2000

While sum(ulP) = 2.0000.

Now, none of the probabilites are greater than 1 on their own. Not even when I set the size of the matrix equal to 100.

What do you think each row represents?

 

[haruspex]

Agreed.

The professor suggest that I set the value of tau equal to 1.

With tau = 1, I get the following column matrix:

ulP =

0.8000
0.6000
0.4000
0.2000

While sum(ulP) = 2.0000.

Now, none of the probabilites are greater than 1 on their own. Not even when I set the size of the matrix equal to 100.

What do you think each row represents?

I don't think the sum of these is interesting. It's hard to know what they mean.

 

[Schwarzschild90]

I'm not even sure what significance the column vector for the ultimate probability has.

 

[Ray Vickson]

I'm not even sure what significance the column vector for the ultimate probability has.

Isn't it explained in the attachment?

You have a "ruin" probability R(i) for each starting state i = 1,2, .., M-1. So, you ought to have a vector of ruin probabilities, whose components are R(i) for i = 1,1, ..., M-1. Of course, for each i we must have P(ruin)+P(not ruined) = 1, but there is no reason at all that the R(i) themselves should have any particular sum, such as 1.

 

[Schwarzschild90]

But I just messaged the professor that I get 1.2 along the diagonal and 0.8 at the end points in a 4x4 matrix, but since probability can't be greater than 1, then I do not think it is correct. That is for the matrix that ulP produces.

 

[haruspex]

I'm not even sure what significance the column vector for the ultimate probability has.

The intended meaning is that each represents the probability of eventual ruin given a certain initial cash level. But since we do not know how to set tau, they are actually some fixed multiple of that. The ratios between them will be meaningful, but not the exact values.
I assume T is a transition matrix. If you were to flip it and run the same analysis you should get the probabilities of ending up at the other absorbing barrier. In this way, you could find, for a given starting point, the two probabilities for the eventual state. If these do not add up to 1, adjust tau accordingly.

 

[Ray Vickson]

But I just messaged the professor that I get 1.2 along the diagonal and 0.8 at the end points in a 4x4 matrix, but since probability can't be greater than 1, then I do not think it is correct. That is for the matrix that ulP produces.

Could you please show an actual example of the 4x4 matrix T you get? And, could you please state exactly is the numerical value of the constant "c" in the differential equation dn(t) = c*T*n(t) that you have solved? Can you tell us exactly what formula you use to get u(i)? Please give actual mathematical expressions, not just Matlab code; not all of us here have access to Matlab or are even very familiar with its syntax.

Give the actual numerical values for the matrix T, and please do NOT attach a file---just type it out here. One easy way to type out a matrix T is to write T = [Row1, Row2, Row3, Row4], where Row_i = [a,b,c,d] . For example, we could write [[1,2,3],[4,5,6],[7,8,9]] for the matrix
$$ \pmatrix{1&2&3\\4&5&6\\7&8&9} $$

 

[Schwarzschild90]

The intended meaning is that each represents the probability of eventual ruin given a certain initial cash level. But since we do not know how to set tau, they are actually some fixed multiple of that. The ratios between them will be meaningful, but not the exact values.
I assume T is a transition matrix. If you were to flip it and run the same analysis you should get the probabilities of ending up at the other absorbing barrier. In this way, you could find, for a given starting point, the two probabilities for the eventual state. If these do not add up to 1, adjust tau accordingly.

I talked to the professor again today. We're only interested in the first row of the inverse matrix T; it gives us the ultimate probability of ruin. By multiplying the first row with n0 - a column vector of all zeros, except in one place, which has a one.

Now, the professor also talked about multiplying the inverse matrix of T by an identiy matrix, which should return what we're interested in.Now,
$$T = \begin{pmatrix} -2 & 1 & 0 & 0 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 0 & 0 & 1 & -2 \\ \end{pmatrix}$$

and
$$T^{-1} = \begin{pmatrix} -0.8 & -0.6 & -0.4 & -0.2 \\ -0.6 & -1.2 & -0.8 & -0.4 \\ -0.4 & -0.8 & -1.2 & -0.6 \\ -0.2 & -0.4 & -0.6 & -0.8 \\ \end{pmatrix}$$

while (matlab notation - hence the odd -0)
$$\textbf{I} = T^{-1}*T = \begin{pmatrix} 1 & 0 & 0 & 0 \\ -0 & 1 & 0 & -0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{pmatrix}$$

 

[Schwarzschild90]

Could you please show an actual example of the 4x4 matrix T you get?

Done.

And, could you please state exactly is the numerical value of the constant "c" in the differential equation dn(t) = c*T*n(t) that you have solved?

Which differential equation is that? The professor has given us all the results.

He metions that the formula (not sure which it is) can be derived analytically, by methods we have not discussed.

Can you tell us exactly what formula you use to get u(i)? Please give actual mathematical expressions, not just Matlab code; not all of us here have access to Matlab or are even very familiar with its syntax.

The formula for U(i) was not determined by me, but by the professor. My assignment is to produce results using the equations and then argue that what I have obtained is correct.

Give the actual numerical values for the matrix T, and please do NOT attach a file---just type it out here. One easy way to type out a matrix T is to write T = [Row1, Row2, Row3, Row4], where Row_i = [a,b,c,d] . For example, we could write [[1,2,3],[4,5,6],[7,8,9]] for the matrix
$$ \pmatrix{1&2&3\\4&5&6\\7&8&9} $$

Done.

~~~~
Please excuse me double posting.

 

[Schwarzschild90]

OK, so I see from this that my initial description of your problem (back in post #6) was exactly right: you do have exponentially-distributed times between successive gambles of size 1. That all follows (essentially as a theorem) from the differential equation itself (but the details of that are a bit lengthy, so I will skip the proof here).

I am not sure what is the source of your worries about probabilities > 1; there are no such things in this problem. The elements of T are not probabilities, they are transition rates. More precisely, the off-diagonal elements are rates of transition to other states, and the absolute values of the diagonal elements are the exponential holding-time parameters, so that a state remains unchanged for an exponentially-distributed amount of time with mean 1/2.

It follows from your T that the "movements" (rather than the "timings") of fhe transitions are given by
$$P = \pmatrix{1&0&0&0\\1/2 & 0 &1/2 & 0\\ 0& 1/2 & 0 & 1/2 \\ 0&0&0&1} $$
In other words, for states i =1 and i = 2 we have
$$\Pr\{ \text{next state} = i \pm 1 |\text{this state} = 1 \} = 1/2 $$
and for states i = 0 or i = 3 we have
$$\Pr\{ \text{no further transitions} | \text{this state} = i \} = 1.$$ Note that the elements of $P$ are probabilities, and the rows sum to 1.

The absorption probabilities are the same as in a discrete-time Markov chain with 1-step transition probability matrix $P$. As a function of the initial state they are (1, 2/3, 1/3, 0).

It is just a bit more involved than what I have worked with in this course. Where can I read more about this? Will you provide me with keywords that I should know about?

But when I told the professor that I had calculated a matrix element with a probability greater than 1, he didn't correct me, save to say that I was right in acknowledging that a matrix element cannot have a probability greater than 1. So I'm lead to assume that the inverse matrix is in fact a probability matrix. Or does U(l) give the ultimate probability of ruin, with T being a transition matrix?

 

[Ray Vickson]

It is just a bit more involved than what I have worked with in this course. Where can I read more about this? Will you provide me with keywords that I should know about?

I did not state that the elements of T are probabilities, but that the elements of the inverse of T are. Incidentally, I see how the same reasoning can be applied here.

The elements of $T^{-1}$ are NOT probabilities. The inverse of T has little to do with anything in this type of problem

I left out two states accidentally in the above; there should be 6 states 0,1,2,3,4,5, and so the 1-step transition probability is 6x6. I have already deleted the initial post, but not quickly enough.

The true matrix T in the state-probability differential equation should be 6x6, with first and last rows = all zeroes:
$$T = \pmatrix{0&0&0&0&0&0\\
1&-2&1&0&0&0\\0&1&-2&1&0&0\\0&0&1&-2&1&0\\0&0&0&1&-2&1\\0&0&0&0&0&0}$$
Note that I am following the usual convention of including the taboo states 0 and 5, but are making them absorbing. Your instructor seems to be following an alternative (less common) approach, which is to omit the taboo states altogether. That can be done, but it requires extra bookkeeping that is handled automatically by including the absorbing state in the state space.

Now the 1-step transition probability matrix governing the movements is
$$P = \pmatrix{1&0&0&0&0&0\\1/2 & 0 &1/2 & 0&0&0\\ 0& 1/2 & 0 & 1/2&0&0 \\ 0&0&1/2&0&1/2&0\\ 0&0&0&1/2&0&1/2\\0&0&0&0&0&1} $$
In other words, for states i =1,2,3,4 we have
$$\Pr\{ \text{next state} = i \pm 1 |\text{this state} = 1 \} = 1/2 $$
and for states i = 0 or i = 5 we have
$$\Pr\{ \text{no further transitions} | \text{this state} = i \} = 1.$$ Note that the elements of $P$ are probabilities, and the rows sum to 1.

The ruin probabilities are the same as in a discrete-time Markov chain with 1-step transition probability matrix $P$. As a function of the initial state they are (1, 4/5,3/5,2/5,1/5,0).

 

[Schwarzschild90]

Computing U(l) using MATLAB I get the following ruin probabilities: (0.8000, 0.6000, 0.4000, 0.2000). I.e. the same values that you do, but I have not built into my matrix the absorbing boundary conditions.

It's comforting to know my approach works. Also, thank you for computing the matrix.

Taking the first and last row equal to zero returns a computational error when using the matrix T to compute U(l). Also, I can no longer take the inverse of T. How do I work around this? (I'm trying to implement the absorbing boundary conditions, i.e. all zeros in the first and last row)

Matlab code:
clear
M=6;
X = ones(M,1);
A = diag(-2*X,0);
clear X
X = ones(M-1,1);
B = diag(X,1);
C = diag(X,-1);
T = A+B+C;
T(1,:)=0;
T(M,:)=0;
clear X B C;

$$T = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix}$$

 

[Ray Vickson]

Computing U(l) using MATLAB I get the following ruin probabilities: (0.8000, 0.6000, 0.4000, 0.2000). I.e. the same values that you do, but I have not built into my matrix the absorbing boundary conditions.

It's comforting to know my approach works. Also, thank you for computing the matrix.

Taking the first and last row equal to zero returns a computational error when using the matrix T to compute U(l). Also, I can no longer take the inverse of T. How do I work around this? (I'm trying to implement the absorbing boundary conditions, i.e. all zeros in the first and last row)

Matlab code:
clear
M=6;
X = ones(M,1);
A = diag(-2*X,0);
clear X
X = ones(M-1,1);
B = diag(X,1);
C = diag(X,-1);
T = A+B+C;
T(1,:)=0;
T(M,:)=0;
clear X B C;

$$T = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix}$$

Your T is a 4x4 submatrix of my T. Your T does have an inverse, but mine does not. As I said before, I don't really see any point in computing the inverse, because there are much easier, standard ways of solving such problems.

One final point to clear up: the ruin probabilities do not depend in any way on the value of $\tau$ you choose. The value of $\tau$ affects the speed at which the limits are attained, but not the actual limits themselves. This is easy to see. From the DE $dP(t)/dt = \tau T P(t)$ it follows that for small $\Delta t > 0$ we have for states $X(t) = 0,1,2, \ldots, 5$ and transition rate matrix $T$ that
$$ \Pr \{ X(t + \Delta t) = j| X(t) = i \} = \begin{cases}
1- 2 \tau \Delta t,& i \neq 0,5, j = i \\
1 \tau \Delta t, & i \neq 0,5, j = i \pm 1\\
1 , & i = 0, 5, j = i \\
0, & i = 0, 5, j \neq i
\end{cases}
$$
Thus, for states i = 1,2,3,4 we have $\Pr \{ X(t+\Delta t) \neq i | X(t) = i \} = 2 \tau \Delta t$. That means that the probability of leaving state $i$ in the time interval from $t$ to $t+\Delta t$ is $\tau 2 \Delta t$ to first order in small $\Delta t.$ That is why the holding time in state $i$ (i.e., the time spent in state $i$ is exponentially-distributed with rate parameter $2 \tau$ and mean $\frac{1}{2\tau}$. However, we have
$$\Pr \{ X(t + \Delta t) =i \pm 1 | X(t + \Delta t) \neq i\: \&\: X(t) = i \} = \frac{1 \tau \Delta a}{ 2 \tau \Delta t} = \frac{1}{2} $$
and
$$\Pr \{ X(t + \Delta t) \neq i \pm 1 | X(t + \Delta t) \neq i\: \&\: X(t) = i \} = 0 $$.
Note that these transition probabilities (for the next "position") are independent of the value of $\tau$. That is why the exact same 1-step discrete-time matrix $P$ is obtained without any reference to the value of $\tau$.

Another way to see this is that (using my $T$ rather than yous) we have
$$P(t) \equiv (P_{ij}(t) ) = \exp(t \tau T), $$
which can be written as $P(t) = \exp( \sigma T)$, where $\sigma = t \tau$. The ruin probabilities are $R(i) = \lim_{t \to \infty} P_{i0}(t)$ (because $P_{i0}(\infty)$ is the probability we are in state 0 at time $\infty$, and that means that we must have hit state 0 at some finite time---and stay there forever more. Anyway, $\lim_{t \to \infty} \cdots = \lim_{\sigma \to \infty} \cdots$, so the ruin probabilities can be computed as though $\tau = 1$.

 

[Schwarzschild90]

Okay. I didn't recognize that using the inverse should be avoided.

But would this count as solving the problem numerically?

 

[Ray Vickson]

Okay. I didn't recognize that using the inverse should be avoided.

But would this count as solving the problem numerically?

I don't know what "this" refers to.

Also, I do not know Matlab syntax, so I have no idea what you are doing. Perhaps if you just stated what you are calculating (in standard mathematical language) I could comment knowledgably.

 

[Schwarzschild90]

I don't know what "this" refers to.

Also, I do not know Matlab syntax, so I have no idea what you are doing. Perhaps if you just stated what you are calculating (in standard mathematical language) I could comment knowledgably.

By "this" I mean the method you employ to calculate the transition probabilities.

Sorry, I misread your post as your stating that most of you are very familiar with MATLAB and its syntax. Hency why I didn't take care to explain the syntax.

The MATLAB code (see below) produces a transition matrix T of size M x M. The matrix is slightly more difficult to produce in matlab, as you must "compose" it of a variety of "constituent" matrices. Irregardless, the following code produces a matrix T equal to the one you produced manually.

clear
M=6;
X = ones(M,1);
A = diag(-2*X,0);
clear X
X = ones(M-1,1);
B = diag(X,1);
C = diag(X,-1);
T = A+B+C;
T(1,:)=0;
T(M,:)=0;
clear X B C;

$$U(l) = \frac{ T^{-1} }{\tau} * n_0$$ correctly produces the ultimate probability, but only if the first and last row are not equal to zero, in which case MATLAB cannot determine the inverse matrix of T.

 

[Schwarzschild90]

I want to expound on the ultimate probability as determined by matlab.

Taking tau = 1, we get the following expression for the ultimate probability, taking l = 1. (So the first position in the n_0 column vector is = 1.)

$$U(l) = - \frac{T^{-1}}{\tau} n_0 = \left[\begin{array}{cc} 0.8 & 0.6 & 0.4 & 0.2 \\ 0.6 & 1.2 & 0.8 & 0.4 \\ 0.4 & 0.8 & 1.2 & 0.6 \\ 0.2 & 0.4 & 0.6 & 0.8 \\ \end{array} \right] \left[\begin{array}{cc} 1 \\ 0 \\ 0 \\ 0 \\ \end{array} \right] = \left[\begin{array}{cc} 0.8 \\ 0.6 \\ 0.4 \\ 0.2 \\ \end{array} \right]$$

Where $n_0$ is the initial state. Do I use n_0 to "pluck" out the first row? As it is, I don't use n_0 for any purpose whatsoever.

But taking out the first row determines the ultimate probability of ruin as a function of the initial state?

$$Pr = \left[\begin{array}{cc} 0.8 & 0.6 & 0.4 & 0.2 \\ \end{array} \right]$$

 

[haruspex]

the ruin probabilities do not depend in any way on the value of ττ\tau you choose. The value of ττ\tau affects the speed at which the limits are attained

Yet the equation used for U(l) in post #34 and earlier divides by tau, and taking different values of tau is producing different probabilities.
I cannot understand why tau features there in conjunction with T. T should already be rate-independent.