simulate a markov chain including the absorbing state "death"

I am simulating a Markov chain process with the transition probability matrix:

P = [0.9044 0.0099 0 0 0.0857;...

0 0.7287 0.1057 0.0799 0.0857;...

0 0.0699 0.8183 0 0.1118; ...

0 0.376 0 0.6199 0.0041; ...

0 0 0 0 1]

where I have five states = single, marriage, divorce, widow and death

When I run the simulations, I always end up with state 5 (which is death).

Should I not include the abosrbing state?

Answer 1

[John D'Errico]

You won't always eventually die? While there are ways for that to not happen eventually, usually they are the plot of some movie. And if you DO know some way to avoid that particular absorbing state, you might consider selling the formula. my guess is there would be a few buyers at any cost.

Seriously, all states eventually allow transition into the absorbing state with some non-zero probability, so the long time state of the system would have you always in the absoring state. This is just an expected behavior for such a process.

[Sehrish Usman]

Yes you are right but I was actually accepting to happen that not so quickly. I am actually using the code you recommended previously. I am running the code by setting the initial condition = 1 (it means non-married). when I run the code for T periods, the individual die within next few years.

I am sure I must be doing some mistake as its the first time I am simulating a Markov process.

Following is my code:

 T = [0.9044 0.0099 0 0 0.0857;0 0.7287 0.1057 0.0799 0.0857;...
             0 0.0699 0.8183 0 0.1118; 0 0.376 0 0.6199 0.0041; 0 0 0 0 1]
CT = [zeros(size(T,1),1),cumsum(T,2)];
numberOfSims = 100;
simTime = 1000;
initialState = 1;
currentState = repmat(initialState,[1,numberOfSims]);
for n = 1:simTime
  r = rand(1,numberOfSims);
  for m = 1:numberOfSims
    % for each sim, where does r(m) lie, relative to the boundaries
    % in the corresponding row of CT?
    [~,~,currentState(m)] = histcounts(r(m),CT(currentState(m),:));
  end
end
finalState = currentState;

[John D'Errico]

At any time state, I would expect if you are not dead already, the probability you will die is very approximately 10%. (Geez, you live in one risky place! Have you considered moving to someplace safer?) I've been married 40 years now. I can't possibly be that much of an outlier. And the last I checked, I think I'm still alive. If I pinch myself, does that work?

Let me look at the expected values for a few time periods out:

P = [1 0 0 0 0;zeros(10,5)];
for i = 1:10
  P(i+1,:) = P(i,:)*T;
end

format short g
P
P =
            1            0            0            0            0
       0.9044       0.0099            0            0       0.0857
      0.81794     0.016168    0.0010464   0.00079101      0.16406
      0.73974      0.02025    0.0025652    0.0017821      0.23566
      0.66902     0.022929    0.0042395    0.0027227      0.30108
      0.60507     0.024652    0.0058927    0.0035198      0.36087
      0.54722     0.025689    0.0074277    0.0041516      0.41551
      0.49491     0.026217    0.0087934    0.0046261      0.46546
      0.44759     0.026358    0.0099668    0.0049625      0.51112
       0.4048     0.026201     0.010942    0.0051823      0.55287
       0.3661     0.025814     0.011723     0.005306      0.59105

Each row of this resulting matrix is the discrete probability distribution you will be in the corresponding state, after that many periods of time. So, after 10 years almost 60% of the time you will be in the absorbing state, given that transition matrix.

You are using a simulation time of 1000 "years". Do you have any idea how unlikely it is to not be in the absorbing state, even after only say 30 such "years"?

We can quickly compute the absorbing state probability after 30, 100, and 1000 time periods as the last element of the corresponding vectors:

>> [1 0 0 0 0]*T^30
ans =
      0.04907    0.0092301     0.007705    0.0023248      0.93167
      
>> [1 0 0 0 0]*T^100
ans =
   4.3257e-05   3.2886e-05   3.4913e-05   8.8373e-06      0.99988

Note that even after 100 years, almost certainly every one of those 100 sims you did will be in the absorbing state.

>> [1 0 0 0 0]*T^1000
ans =
   2.2938e-44   2.9041e-41   3.3015e-41   7.9635e-42            1

I have not looked carefully at the code you show, nor have I verified that it produces consistently the results I have just shown. But what you describe is entirely consistent with my expectations, given that Markov transition matrix.

[Sehrish Usman]

Thank you so much that perfectly makes sense. I just have another question relevant to your code.

You always use transition probability matrix (to simulate). Is there any possibility to simulate using transition rates? I have tried my level best to calculate transition probabilities using various approaches but it either does not come correct or probability outcomes just does not make any sense.

I would be really greatful if you just tell that?

[John D'Errico]

I'm not sure what you are asking. Are you saying that you do not know how to generate a viable transition matrix from data, say from marriage, and death statistics? Or are you talking about a formulation of such a process in the form of a system of ordinary differential equations, where now you would have a death rate, marriage rate, etc., each of which would depend on the state you are currently in.

[Sehrish Usman]

Yes the second case you mentioned. I have a transition rate matrix for four states (nonmarried, married, divorced, widow) and each element of this matrix shows the marriage rate, divorce rate, mortality rate etc. depending on age and current state.

I want to simulate the markov model but I need to calculate the transition probability matrix. I was not sure if I can use directly these rates to calculate probability?

I cannot use simple Probability = expm(Q) ........ as the transitons states are more than two for each node.

any idea? I would really appreciate if you have any suggestion.