Stoch Proc - Lab/Tutorial 8 Solutions

Scott Alexander

Q1

(a)

Possible state transitions:
*RR(1) -> RR(1) or RN(2)
*RN(2) -> NR(3) or NN(4)
*NR(3) -> RR(1) or RN(2)
*NN(4) -> NR(3) or NN(4)

We are given:
*P(RR(1)->RR(1))=0.4  ⟹  P(RR(1)->RN(2))=0.6
*P(RN(2)->NR(3))=0.2  ⟹  P(RN(2)->NN(4))=0.8
*P(NR(3)->RR(1))=0.3  ⟹  P(NR(3)->RN(2))=0.7
*P(NN(4)->NR(3))=0.1  ⟹  P(NN(4)->NN(4))=0.9
P ≡ P(1) = (Subscript[[Prob((Subscript[X, t=1])=j|(Subscript[X, t])=i)], i,j]) = (Subscript[p, ij])(1) = (({{(4/10), (6/10), 0, 0}, {0, 0, (2/10), (8/10)}, {(3/10), (7/10), 0, 0}, {0, 0, (1/10), (9/10)}}))

(b)

We use
p(t) = ((((P^T)))^t)p(0) = ((((P^t)))^T)p(0)
and then
𝔼[(Subscript[X, t])] = (∑_(x=1)^4)x*P((Subscript[X, t])=x) = (∑_(x=1)^4)x*(Subscript[p, x])(t).

Calculations

Clear[P1,p0];
P1=(({{(4/10), (6/10), 0, 0}, {0, 0, (2/10), (8/10)}, {(3/10), (7/10), 0, 0}, {0, 0, (1/10), (9/10)}})); (* P(1) *)
p0={(10/100),(15/100),(20/100),(55/100)}; (* p(0) *)
x={1,2,3,4}; (* state vector *)

p4=MatrixPower[Transpose[P1],4].p0; (* p(4) *)
x.p4//N (* E[X(4)] *)

p21=MatrixPower[Transpose[P1],21].p0; (* p(21) *)
x.p21//N (* E[X(21)] *)

p52=MatrixPower[Transpose[P1],52].p0; (* p(52) *)
x.p52//N (* E[X(52)] *)

(c)

We need to solve
(P^T)π = π
which is
((P^T)-I)π = 0
subject to
1·π = 1.    (RHS "1" a vector, LHS a scalar)

Method 1. 
Solve the linear system
(({{(P^T)}, {1}}))π = (({{0}, {1}}))
which is
(({{-(6/10), 0, (3/10), 0}, {(6/10), -1, (7/10), 0}, {0, (2/10), -1, (1/10)}, {0, (8/10), 0, -(1/10)}, {1, 1, 1, 1}}))(({{(Subscript[π, 1])}, {(Subscript[π, 2])}, {(Subscript[π, 3])}, {(Subscript[π, 4])}})) = (({{0}, {0}, {0}, {0}, {1}}))
or
(({{-(6/10), 0, (3/10), 0, 0}, {(6/10), -1, (7/10), 0, 0}, {0, (2/10), -1, (1/10), 0}, {0, (8/10), 0, -(1/10), 0}, {1, 1, 1, 1, 1}})) ~ (({{1, 0, 0, 0, (1/21)}, {0, 1, 0, 0, (2/21)}, {0, 0, 1, 0, (2/21)}, {0, 0, 0, 1, (16/21)}, {0, 0, 0, 0, 0}})).

So    (({{(Subscript[π, 1])}, {(Subscript[π, 2])}, {(Subscript[π, 3])}, {(Subscript[π, 4])}}))=(({{(1/21)}, {(2/21)}, {(2/21)}, {(16/21)}})) ≈ (({{0.047619047619047616`}, {0.09523809523809523`}, {0.09523809523809523`}, {0.7619047619047619`}})).

Method 2.
Find the eigenvector corresponding to the unit eigenvalue
and normalise so elements are positive and sum to one.

Calculations

(* solving with Gauss-Jordan elimination *)

RowReduce[(({{-(6/10), 0, (3/10), 0, 0}, {(6/10), -1, (7/10), 0, 0}, {0, (2/10), -1, (1/10), 0}, {0, (8/10), 0, -(1/10), 0}, {1, 1, 1, 1, 1}}))]//MatrixForm

(* solving linear system *)
Clear[m,b];
m=Append[Transpose[P1]-IdentityMatrix[4],{1,1,1,1}];
b={0,0,0,0,1};
LinearSolve[m,b]//MatrixForm

(* finding eigenvector *)
eigensys=Eigensystem[Transpose[P1]];
eigensys//N//MatrixForm
eigensys[[2,1]]/Total[eigensys[[2,1]]]//MatrixForm

(d)

From (b) we know that all entries of P(4) are positive, so by Theorem 6 page 35 Chapter 5 Notes
(Subscript[X, t]) is ergodic.

We can also see this as
(lim_(t->∞)) p(t) = π
and does not depend on p(0) (see code below).

Note that
(lim_(t->∞)) p(t) = (lim_(t->∞))((((P^T)))^t)p(0)
  = (lim_(t->∞))(((({{(4/10), 0, (3/10), 0}, {(6/10), 0, (7/10), 0}, {0, (2/10), 0, (1/10)}, {0, (8/10), 0, (9/10)}})))^t)(({{(Subscript[p, 1])(0)}, {(Subscript[p, 2])(0)}, {(Subscript[p, 3])(0)}, {1-(∑_(i=1)^3)(Subscript[p, i])(0)}}))
  = (({{(1/21), (1/21), (1/21), (1/21)}, {(2/21), (2/21), (2/21), (2/21)}, {(2/21), (2/21), (2/21), (2/21)}, {(16/21), (16/21), (16/21), (16/21)}}))(({{(Subscript[p, 1])(0)}, {(Subscript[p, 2])(0)}, {(Subscript[p, 3])(0)}, {1-(∑_(i=1)^3)(Subscript[p, i])(0)}})) = (({{(1/21)}, {(2/21)}, {(2/21)}, {(16/21)}})) ≈ (({{0.047619047619047616`}, {0.09523809523809523`}, {0.09523809523809523`}, {0.7619047619047619`}}))
which does not depend on p(0).

This limit is also the solution of
(P^T)π = π    s.t.    1·π=1.

Calculations

Clear[p0,p01,p02,p03,s,t];
p0={p01,p02,p03,1-p01-p02-p03};
Limit[MatrixPower[Transpose[P1],s],s->∞]//MatrixForm
Limit[MatrixPower[Transpose[P1],t].p0,t->∞]//MatrixForm

Q2

(a)

Clear[P1];
P1=(({{0, 0, 0, 1}, {(1/5), 0, (4/5), 0}, {0, 0, 0, 1}, {0, 1, 0, 0}}));
chain=DiscreteMarkovProcess[{1/4,1/4,1/4,1/4},P1];
Graph[chain]
MarkovProcessProperties[chain,"Periods"]

In this case Markov chain is not ergodic as the limiting distribution (lim_(t->∞))p(t) does not exist.
This is becasue all states have period of 3 (see calculations below).

Calculations

Clear[p01,p02,p03];
p0={p01,p02,p03,1-p01-p02-p03};
Table[MatrixPower[Transpose[P1],t].p0//MatrixForm,{t,1,6}]

Table[MatrixPower[P1,s]//MatrixForm,{s,1,6}]

(b)

Clear[P1];
P1=(({{(1/6), (1/2), 0, 0, (1/3)}, {(1/4), (1/4), 0, 0, (1/2)}, {0, 0, (1/2), (1/2), 0}, {0, 0, (6/7), 0, (1/7)}, {(1/10), (3/10), 0, 0, (3/5)}}));
chain=DiscreteMarkovProcess[{1/5,1/5,1/5,1/5,1/5},P1];
Graph[chain]
MarkovProcessProperties[chain,"Periods"]

In this case the Markov chain is ergodic.
Consider the arguments below.

METHOD 1
Irrespective of the intial distribution,the Markov chain will eventually enter
state 5, from which it will never exit the class {1,2,5}.

To see this, consider the chain when (Subscript[X, t])∈{1,2,5}.
In this case (Subscript[X, t+s])∈{1,2,5} for all s>0 as it is in an absorbing class.

If (Subscript[X, t])=3, then it will eventually hit (Subscript[X, t+s])=4.
If (Subscript[X, t+s])=4|(Subscript[X, t])=3 did not occur for all s>0 then this would imply (Subscript[p, 34])=0 which is a contradiction.

If  (Subscript[X, t])=4, then it will eventually hit (Subscript[X, t+s])=5.
If (Subscript[X, t+s])=5|(Subscript[X, t])=4 did not occur for all s>0 then this would imply (Subscript[p, 45])=0 which is a contradiction.

So irrespective of starting distribution p(0), the chain will eventually hit the absorbing class {1,2,5}.

So we can consider the 3-state Markov chain involving the states 1,2 and 5 and ask if this is ergodic.
The answer is yes as its transition probability matrix is
(P^((1,2,5)))(1) = (({{(1/6), (1/2), (1/3)}, {(1/4), (1/4), (1/2)}, {(1/10), (3/10), (3/5)}}))
which has all elements greater than zero.
So by Theorem 6 (Chapter 5 Notes) the chain is ergodic.

METHOD 2
Another way to show this is to use the definition of ergodicity:
that (lim_(t->∞))p(t) exists and does not depend on p(0).

We can see this holds by having Mathematica perform the calculations (see below).

Note that
(lim_(t->∞)) p(t) = (lim_(s->∞))((((P^T)))^t)p(0)
  = (lim_(t->∞))(((({{(1/6), (1/4), 0, 0, (1/10)}, {(1/2), (1/4), 0, 0, (3/10)}, {0, 0, (1/2), (6/7), 0}, {0, 0, (1/2), 0, 0}, {(1/3), (1/2), 0, (1/7), (3/5)}})))^t)(({{(Subscript[p, 1])(0)}, {(Subscript[p, 2])(0)}, {(Subscript[p, 3])(0)}, {(Subscript[p, 4])(0)}, {1-(∑_(i=1)^4)(Subscript[p, i])(0)}}))
  = (({{(3/19), (3/19), (3/19), (3/19), (3/19)}, {(6/19), (6/19), (6/19), (6/19), (6/19)}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {(10/19), (10/19), (10/19), (10/19), (10/19)}}))(({{(Subscript[p, 1])(0)}, {(Subscript[p, 2])(0)}, {(Subscript[p, 3])(0)}, {(Subscript[p, 4])(0)}, {1-(∑_(i=1)^4)(Subscript[p, i])(0)}})) = (({{(3/19)}, {(6/19)}, {0}, {0}, {(10/19)}}))
which does not depend on p(0).

This limit is also the solution of
(P^T)π = π    s.t.    1·π=1.

Calculations

(* check ergodicity *)
Clear[p01,p02,p03,p04,s,t];
p0={p01,p02,p03,p04,1-p01-p02-p03-p04};
Limit[MatrixPower[Transpose[P1],s],s->∞]//FullSimplify//MatrixForm
Limit[MatrixPower[Transpose[P1],t],t->∞].p0//FullSimplify//MatrixForm

(* solve stationary distribution eigenvector problem *)
eigensys=Eigensystem[Transpose[P1]];
eigensys//N//MatrixForm
eigensys[[2,1]]/Total[eigensys[[2,1]]]//FullSimplify//MatrixForm

Q3

(a)

The intensity matrix
A=(({{-4, 3, (4/5), (1/5)}, {(6/4), -5, (11/4), (3/4)}, {1, 2, -7, 4}, {0, 0, 0, 0}})) = (({{-(Subscript[ν, 1]), (Subscript[⋁, 1])(Subscript[p, 1,2]), (Subscript[⋁, 1])(Subscript[p, 1,3]), (Subscript[⋁, 1])(Subscript[p, 1,4])}, {(Subscript[⋁, 2])(Subscript[p, 2,1]), -(Subscript[ν, 2]), (Subscript[⋁, 2])(Subscript[p, 2,3]), (Subscript[⋁, 2])(Subscript[p, 2,4])}, {(Subscript[⋁, 3])(Subscript[p, 3,1]), (Subscript[⋁, 3])(Subscript[p, 3,2]), -(Subscript[ν, 3]), (Subscript[⋁, 3])(Subscript[p, 3,4])}, {0, 0, 0, -(Subscript[ν, 4])}}))
so the transition matrix applying when a jump occurs is
(P^jump) = (({{0, ((Subscript[a, 1,2])/(Subscript[ν, 1])), ((Subscript[a, 1,3])/(Subscript[ν, 1])), ((Subscript[a, 1,4])/(Subscript[ν, 1]))}, {((Subscript[a, 2,1])/(Subscript[ν, 2])), 0, ((Subscript[a, 2,3])/(Subscript[ν, 2])), ((Subscript[a, 2,4])/(Subscript[ν, 2]))}, {((Subscript[a, 3,1])/(Subscript[ν, 3])), ((Subscript[a, 3,2])/(Subscript[ν, 3])), 0, ((Subscript[a, 3,4])/(Subscript[ν, 3]))}, {0, 0, 0, 1}})) = (({{0, (15/20), (4/20), (1/20)}, {(6/20), 0, (11/20), (3/20)}, {(1/7), (2/7), 0, (4/7)}, {0, 0, 0, 1}})).

(b)

Clear[AA];
AA=(({{-4, 3, (4/5), (1/5)}, {(6/4), -5, (11/4), (3/4)}, {1, 2, -7, 4}, {0, 0, 0, 0}}));
chain=ContinuousMarkovProcess[{1/4,1/4,1/4,1/4},AA];
Graph[chain]

METHOD 1
(Subscript[X, t]) is ergodic as it has a single absorbing state that can be reached from all other states (see arguments Q1(b)).
So the limiting distribution will be
(lim_(t->∞))p(t) = (({{0}, {0}, {0}, {1}}))
which is the unique stationary solution.

METHOD 2
Another way to show this is to use the definition of ergodicity:
that (lim_(t->∞))p(t) exists and does not depend on p(0).

We can see this holds by having Mathematica perform the calculations (see below).

Note that
(lim_(t->∞)) p(t) = (lim_(s->∞))((((P^T)))^t)p(0)
  = (({{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 1, 1, 1}}))(({{(Subscript[p, 1])(0)}, {(Subscript[p, 2])(0)}, {(Subscript[p, 3])(0)}, {1-(∑_(i=1)^3)(Subscript[p, i])(0)}})) = (({{0}, {0}, {0}, {1}}))
which does not depend on p(0).

This limit is also the solution of
(A^T)π = 0    s.t.    1·π=1.

Stationary distribution calculations

(* check ergodicity *)
Clear[PP,s,t];
PP[s_]:=MatrixExp[s AA];
Clear[p01,p02,p03,s,t];
p0={p01,p02,p03,1-p01-p02-p03};
Limit[Transpose[PP[s]],s->∞]//MatrixForm
Limit[Transpose[PP[t]],t->∞].p0//MatrixForm

(* solve stationary distribution problem *)
LinearSolve[Append[Transpose[AA],{1,1,1,1}],{0,0,0,0,1}]//Chop//MatrixForm

(c)

chain=ContinuousMarkovProcess[{1/4,1/4,1/4,1/4},AA];
Graph[chain]

MarkovProcessProperties[chain,"CommunicatingClasses"]

MarkovProcessProperties[chain,"TransientClasses"]

MarkovProcessProperties[chain,"AbsorbingClasses"]

MarkovProcessProperties[chain,"RecurrentClasses"]

States 1, 2 and 3 are transient states. 
State 4 is an absorbing and recurrent state.

States {1,2,3} form a class.
State {4} forms a class.

(d)

Plot[{PP[s][[1,4]],PP[s][[2,4]],PP[s][[3,4]]},{s,0,10},PlotLegends->{"(Subscript[p, 1, 4])(s)","(Subscript[p, 2, 4])(s)","(Subscript[p, 3, 4])(s)"},AxesLabel->{"s","(Subscript[p, i, 4])(s)"},PlotRange->All]

(e)

Note that because state 4 is absorbing
P(τ≤12|(Subscript[X, 1])=3) = P((Subscript[X, 12])=4|(Subscript[X, 1])=3) ≈ 1
and
P(τ≤17|(Subscript[X, 3])=2) = P((Subscript[X, 17])=4|(Subscript[X, 3])=2) ≈ 1.

Calculations

{PP[11][[3,4]],PP[14][[2,4]]}//N

Q4

(a)

First note that
𝔼[(Subscript[T, 1])]=2    ⟹    (1/(Subscript[λ, 1]))=2    ⟹    (Subscript[λ, 1])=(1/2)    ⟹    (Subscript[T, 1])~exp((1/2)),
var((Subscript[T, 2]))=4    ⟹    (1/(Subscript[λ, 2]^2))=4    ⟹    (Subscript[λ, 2])=(1/2)    ⟹    (Subscript[T, 2])~exp((1/2))
and
𝔼[(Subscript[T, 3]^2)]=(1/8)    ⟹    (2/(Subscript[λ, 3]^2))=(1/8)    ⟹    (Subscript[λ, 3])=4    ⟹    (Subscript[T, 3])~exp(4).

So
A= (({{-(Subscript[ν, 1]), (Subscript[⋁, 1])(Subscript[p, 1,2]), (Subscript[⋁, 1])(Subscript[p, 1,3])}, {(Subscript[⋁, 2])(Subscript[p, 2,1]), -(Subscript[ν, 2]), (Subscript[⋁, 2])(Subscript[p, 2,3])}, {(Subscript[⋁, 3])(Subscript[p, 3,1]), (Subscript[⋁, 3])(Subscript[p, 3,2]), -(Subscript[ν, 3])}})) = (({{-(1/2), (1/2)(2/3), (1/2)(1/3)}, {(1/2)(1/2), -(1/2), (1/2)(1/2)}, {4(3/4), 4(1/4), -4}})) = (({{-(1/2), (1/3), (1/6)}, {(1/4), -(1/2), (1/4)}, {3, 1, -4}})).

Calculations

Expectation[t,t\[Distributed]ExponentialDistribution[λ]]

Expectation[(((t-(1/λ)))^2),t\[Distributed]ExponentialDistribution[λ]]

Expectation[(t^2),t\[Distributed]ExponentialDistribution[λ]]

(b)

By page 146 Lecture Notes Part 1 ,a stationary distribution can be found by solving
(A^T)π = 0
subject to
1·π = 1.

Calculations

Clear[AA];
AA={{-1/2,1/3,1/6},{1/4,-1/2,1/4},{3,1,-4}};LinearSolve[Append[Transpose[AA],{1,1,1}],{0,0,0,1}]//MatrixForm

(c)

chain=ContinuousMarkovProcess[{1/3,1/3,1/3},AA];
Graph[chain]

We can determine ergodicity by checking if (lim_(t->∞))p(t) exists and does not depend on p(0).
We can see this holds by having Mathematica perform the calculations (see below).

This limit is also the solution of
(A^T)π = 0    s.t.    1·π=1
which was found in (b).

Note that
(lim_(t->∞))p(t) = (lim_(t->∞)) P(((t))^T)p(0) = π.

To see why the dependence on t vanishes for (lim_(t->∞)) P(t) consider for instance
(Subscript[p, 11])(t) = (21/41)-((5 (e^(-((5 t)/2)-(Sqrt[(17/6)]) t)))/(2 (34+5 (Sqrt[102]))))+(((Sqrt[(51/2)]) (e^(-((5 t)/2)-(Sqrt[(17/6)]) t)))/(34+5 (Sqrt[102])))+((5 (e^(-((5 t)/2)+(Sqrt[(17/6)]) t)))/(2 (-34+5 (Sqrt[102]))))+(((Sqrt[(51/2)]) (e^(-((5 t)/2)+(Sqrt[(17/6)]) t)))/(-34+5 (Sqrt[102]))).

Then in the limit the exponential terms vanish as  t->∞.

Clear[PP];
PP[s_]:=MatrixExp[s AA];
PP[s][[1,1]]//MatrixForm

Calculations

(* check ergodicity *)
Clear[PP,s,t];
PP[s_]:=MatrixExp[s AA];
Clear[p01,p02,s,t];
p0={p01,p02,1-p01-p02};
Limit[PP[s],s->∞]//MatrixForm
Limit[Transpose[PP[t]],t->∞].p0//FullSimplify//MatrixForm