# Find expected value of time to reach a state in Markov chain, by simulation

Consider a time homogeneous Markov chain $(X_n)_{n=0}$ with state space $E$, initial distribution $p(0)$ and transition probability matrix $P$ given by $E = \{0, 1, 2\}, p(0) = [1\;\; 0\;\; 0]$ and

P= $\begin{bmatrix}1/2 & 1/3 & 1/6\\0 & 2/3 & 1/3 \\0 & 0 & 1 \end{bmatrix}$

respectively. Find by a computer simulations an as good as is possible approximation of the expected value $\Bbb E(T)$ of the time $T = \min\{ n ∈ N : X_n = 2 \}$ it takes the chain to reach the state 2. Someone who can help me in doing some kind of pseudocode for this question, or some matlab code?

• Do you know how to sample from a discrete distribution? For instance, the second value of the process is given by sampling from the discrete distribution $(1/2,1/3,1/6)$. Then once you have the second value you sample from another distribution, etc.
– Ian
Dec 14, 2015 at 22:05
• No sorry. I have not taken any programming course yet. This is the main problem for me at this particulat moment. That is to solve some of these questions and get some kind of intuitive feeling in doing so for it. And taking a programming course after this christmas. I appreciate your'e answer but maybe you could if you have time demonstrate by some matlab code?
– user297559
Dec 14, 2015 at 22:14
• Do you have to use matlab? Is maple or C/C+/C++ also acceptable? Dec 14, 2015 at 22:24
• I prefer matlab in order to understand it myself. And the problem is for myself, to get a starter or intuition of maybe solve similar problems.
– user297559
Dec 14, 2015 at 22:29
• Then I will leave it to others :). Dec 14, 2015 at 22:31

## 3 Answers

In Mathematica:

myP = {{1/2, 1/3, 1/6}, {0, 2/3, 1/3}, {0, 0, 1}};

ListPlot[Transpose[NestList[#.myP &, {1, 0, 0}, 10]],
AxesLabel -> {"step", "Probability"},
Joined -> True,
PlotLegends -> {"P[0]", "P[1]", "P[2]"} ]


The successive differences in probabilities for being in state $2$ are:

difProbs = Differences@NestList[#.myP &, {1, 0, 0}, 100][[All, 3]];


The expected value is:

N@Range[100].difProbs


(* 4 *)

Probably easier for you to understand R than mathematica if you know matlab.

Ask if you dont understand what im doing.

Here is the code in R:

N = 10^4

T<-vector(length=N)

ptm <- proc.time()

for (k in 1:N) {
i<-0
X<-0

while (X==0) {
u=runif(1)

if (1/6<u & u<1/2) {X<-1}
else if (u<1/6) {X<-2}
else {X<-0}

i <- i+1
}

while (X==1) {
u=runif(1)

if (u<1/3) {X<-2}

i <- i+1
}

# Done Save

T[k]<-i

}

print(proc.time()-ptm)

mean(T)

• While X is in state 0 or 1, respectively, i generate uniform random numbers , u, and for the X=1 case i check if this u is less than 1/3. There is a 1/3 probability to go from 1 to 2 so if u fulfills this condition i set X=2 and we have reached our destination. If the condition is not fulfilled it stays in X=1 until i get a u that fulfills the condition. I add one to i for each time i check this. I do the same thing for X=0 and ofc always start in X=0. This procedure i repeat 10^4 times in the for-loop and for each time i save the i's and then i take the mean to get the simulated average time. Dec 15, 2015 at 9:20
• I added this comment before you removed yours. I leave it there if someone else need an explanation. Dec 15, 2015 at 9:22

Let $\psi(i)$ be the expected time to reach state 2 from state i, where i ={0,1,2}.

Then $\psi(2) =0$

$\psi(1) = 1+ \frac{2}{3}\psi(1)+\frac{1}{2}\psi(2)$

$\psi(0) = 1+\frac{1}{2}\psi(0)+\frac{1}{3}\psi(1)+\frac{1}{6}\psi(2)$

Solve the system to find $\psi(0)$

The initial state is 0 and the expected time to reach state 2 from 0 is $\boxed{4}$

• They asked to do it by simulation, though.
– Ian
Dec 14, 2015 at 22:19
• @DavidG.Stork This calculates the expected time to achieve the first hitting of state 2 (possibly doing so by passing through state 1).
– Ian
Dec 14, 2015 at 23:05
• From the transition matrix, I get to understand that state 2 is absorbing with 1 to itself and 0 elsewhere. If it reached 2, it ain't going anywhere Dec 14, 2015 at 23:06