Consider the following Markov process: It takes values on {0,1,2} and has transition matrix \[ P = \begin{bmatrix} 0 & 1 & 0 \\ \frac12 & 0 & \frac12 \\ 0 & 1 & 0 \\ \end{bmatrix} \] Find the n-step transition matrix and the stationary distribution of this process both analytically and via a simulation.

\[ \begin{aligned} &\vert P-\lambda I \vert = \vert \begin{bmatrix} -\lambda & 1 & 0 \\ \frac12 & -\lambda & \frac12 \\ 0 & 1 & -\lambda \\ \end{bmatrix} \vert= \\ &-\lambda \vert \begin{bmatrix} -\lambda & 1/2 \\ 1 & -\lambda \\ \end{bmatrix} \vert - \frac12 \vert \begin{bmatrix} 1 & 0 \\ 1 & -\lambda \\ \end{bmatrix} \vert= \\ &-\lambda(\lambda^2-1/2)-\frac12(-\lambda))=\\ &\lambda\left(-\lambda^2+1\right)= \lambda(\lambda-1)(\lambda+1)=0 \end{aligned} \] so the eigenvalues are -1, 1 and 0. For the eigenvectors we find

\[ \begin{aligned} &(P-(-1)I)x = \begin{bmatrix} 1 & 1 & 0 \\ \frac12 & 1 & \frac12 \\ 0 & 1 & 1 \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2\\ \end{bmatrix} = \begin{bmatrix} x_0 +x_1\\ x_0/2+x_1+x_2/2 \\ x_1+x_2\\ \end{bmatrix} = \begin{bmatrix} 0\\ 0 \\ 0\\ \end{bmatrix}\\ &x_0 = 1;x_1=-1;x_2=1\\ &\sqrt{1^2+(-1)^2+1^2} = \sqrt 3\\ &x=(1/\sqrt 3, -1/\sqrt 3, 1/\sqrt 3)^T \end{aligned} \] \[ \begin{aligned} &(P-1I)x = \begin{bmatrix} -1 & 1 & 0 \\ \frac12 & -1 & \frac12 \\ 0 & 1 & -1 \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2\\ \end{bmatrix} = \begin{bmatrix} -x_0 +x_1\\ x_0/2-x_1+x_2/2 \\ x_1-x_2\\ \end{bmatrix} = \begin{bmatrix} 0\\ 0 \\ 0\\ \end{bmatrix}\\ &x_0 = 1;x_2=1;x_1=1\\ &\sqrt{1^2+1^2+1^2} = \sqrt 3\\ &x=(1/\sqrt 3, 1/\sqrt 3, 1/\sqrt 3)^T \end{aligned} \]

\[ \begin{aligned} &(P-0I)x = \begin{bmatrix} 0 & 1 & 0 \\ \frac12 & 0 & \frac12 \\ 0 & 1 & 0 \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2\\ \end{bmatrix} = \begin{bmatrix} x_1\\ x_0/2+x_2/2 \\ x_1\\ \end{bmatrix} = \begin{bmatrix} 0\\ 0 \\ 0\\ \end{bmatrix}\\ &x_1 =1; x_2=-1;x_1=0\\ &\sqrt{1^2+0^2+(-1)^2} = \sqrt 2\\ &x=(1/\sqrt 2, 0, -1/\sqrt 2)^T \end{aligned} \]

so

\[ U= \begin{bmatrix} 1/\sqrt 3 & 1/\sqrt 3 & -1/\sqrt2 & \\ -1/\sqrt 3 & 1/\sqrt 3 & 0 & \\ 1/\sqrt 3 & 1/\sqrt 3 & -1/\sqrt2 & \\ \end{bmatrix} \]

Using the adjoint we can find the inverse is

\[ U^{-1}= \begin{bmatrix} \sqrt3/4 & -\sqrt3/2 & \sqrt3/4 \\ \sqrt3/4& \sqrt3/2 & \sqrt3/4 \\ 1/\sqrt 2 & 0 & -1/\sqrt 2 \\ \end{bmatrix} \]

and so we find the n-step transition matrix to be

\[ \begin{aligned} &P^{(n)} = U D^nU^{-1}\\ &\begin{bmatrix} 0 & 1 & 0 \\ \frac12 & 0 & \frac12 \\ 0 & 1 & 0 \\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & (-1)^n & 0\\ 0 & 0 & 0\\ \end{bmatrix} \begin{bmatrix} \sqrt3/4 & -\sqrt3/2 & \sqrt3/4 \\ \sqrt3/4& \sqrt3/2 & \sqrt3/4 \\ 1/\sqrt 2 & 0 & -1/\sqrt 2 \\ \end{bmatrix} = \\ &\begin{bmatrix} (1+(-1)^n)/2 & (1+(-1)^{n+1})/2 & (1+(-1)^{n})/2 \\ 1/2 & (1+(-1)^{n})/2 & 1/2 \\ (1+(-1)^{n})/2 & (1+(-1)^{n+1})/2 & (1+(-1)^{n})/2 \\ \end{bmatrix} \end{aligned} \]

notice that this oscillates, so P(2n)=P.

For the stationary distribution we need to solve the equation \(\pi P=\pi\), so

\[ \begin{aligned} &\pi_1/2 = \pi_0\\ &\pi_0+\pi_2 = \pi_1\\ &\pi_1/2 = \pi_2\\ \text{so} \\ &\pi_1/2+\pi_1 + \pi_1/2 = 2\pi_1=1\\ &\pi_0=\pi_2 = \frac14; \pi_1 = \frac12\\ \end{aligned} \]

Now via simulation:

gen.x=function(n=1e4) {
  x=rep(0, n)
  for(i in 2:n) {
    if(x[i-1]==0) x[i]=1
    if(x[i-1]==1) x[i]=sample(c(0,2), size=1)
    if(x[i-1]==2) x[i]=1
  }
  x
}
find.Pn=function(n=1) {
  x=gen.x()
  P=matrix(0, 3, 3)
  for(i in 1:(length(x)-n)) 
    P[x[i]+1, x[i+n]+1]=P[x[i]+1, x[i+n]+1]+1
  for(i in 1:3) P[i, ]=P[i, ]/sum(P[i, ])
  P
}

to verify letโ€™s check n=1,2, 3:

find.Pn(1)
##          [,1] [,2]     [,3]
## [1,] 0.000000    1 0.000000
## [2,] 0.504901    0 0.495099
## [3,] 0.000000    1 0.000000
find.Pn(2)
##           [,1] [,2]      [,3]
## [1,] 0.4875201    0 0.5124799
## [2,] 0.0000000    1 0.0000000
## [3,] 0.5061630    0 0.4938370
find.Pn(3)
##          [,1] [,2]     [,3]
## [1,] 0.000000    1 0.000000
## [2,] 0.505002    0 0.494998
## [3,] 0.000000    1 0.000000

and the stationary distribution:

table(gen.x())/1e4
## 
##      0      1      2 
## 0.2521 0.5000 0.2479