# Example 1

## Transition matrix

Let us consider the Markov chain with the following transition matrix $$M = \left(\begin{matrix} 0 & 1/2 & 1/2 & 0 & 0 & 0 \ 0 & 0 & 1-p & 0 & 0 & p \ 0 & 1-q & 0 & q & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & 0 \ 0 & 0 & 0 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 1 \ \end{matrix}\right)$$ with $1>p>0$ close to zero, e.g. $p=1\%$ and $1>q>0$, e.g. $q=0.2$.

## R implementation

From the diagram or $M$, it is easy that there are absorbing states (6) and recurrent states (4,5).

library(markovchain)
M<-matrix(c(0,0.5,0.50,0.0,0,0.00,0,0.0,0.99,0.0,0,0.01,0,0.8,0.00,
0.2,0,0.00,0,0.0,0.00,0.0,1,0.00,
0,0.0,0.00,1.0,0,0.00,0,0.0,0.00,0.0,0,1.00), nrow=6, byrow=TRUE)
colnames(M) <- rownames(M) <- paste0("s",1:6)
chain1 <- as(M,"markovchain")
plot(chain1)


## Invariant measure

We need to solve the following system in order to get the invariant measure $$\mu M = \mu \Leftrightarrow \left{\begin{array}{l} 0 = \mu_1 \ 0.5\mu_1 + (1-q)\mu_3 = \mu_2 \ 0.5\mu_1 + (1-p)\mu_2 = \mu_3 \ q\mu_3+\mu_5= \mu_4 \ \mu_4=\mu_5 \ p\mu_2 + \mu_6=\mu_6 \end{array}\right. \Leftrightarrow \left{\begin{array}{l} \mu_1 = 0 \ (1-q)\mu_3 = \mu_2 \ (1-p)\mu_2 = \mu_3 \ q\mu_3+\mu_4= \mu_4 \ \mu_5=\mu_4 \ p\mu_2 =0 \end{array}\right. \Leftrightarrow \left{\begin{array}{l} \mu_1 = 0 \ \mu_2=0 \ \mu_3=0 \ \mu_5=\mu_4 \ \end{array}\right.$$ Using the probability condition $\sum_i \mu_i=1$, we get $2\mu_4+\mu_6=1$. So there exists an infinite number of invariant measures of the form (as long as $p>0$) $$\mu = (0, 0, 0, \nu, \nu, 1-2\nu), \nu\in[0,1/2].$$ Note that there are independent of the value of $p$ and $q$. Indeed, we have

mu <- c(0, 0, 0, 1/4, 1/4, 1/2)
mu %*% M
mu <- c(0, 0, 0, 0, 0, 1)
mu %*% M
mu <- c(0, 0, 0, 1/2, 1/2, 0)
mu %*% M


## Classification of states

Let $T^{x\rightarrow x}$ is the number of periods to go back to state $x$ knowing that the chain starts in $x$.

• A state $x$ is recurrent if $P(T^{x\rightarrow x}<+\infty)=1$ (equivalently $P(T^{x\rightarrow x}=+\infty)=0$).
• A state $x$ is null recurrent if in addition $E(T^{x\rightarrow x})=+\infty$.
• A state $x$ is positive recurrent if in addition $E(T^{x\rightarrow x})<+\infty$.
• A state $x$ is absorbing if in addition $P(T^{x\rightarrow x}=1)=1$.
• A state $x$ is transient if $P(T^{x\rightarrow x}<+\infty)<1$ (equivalently $P(T^{x\rightarrow x}=+\infty)>0$).

On the example above, using the transition matrix, we have $$\begin{array}{l} \forall k>0, P(T^{s1\rightarrow s1}=k) = 0 \Rightarrow P(T^{s1\rightarrow s1}=+\infty)=1, \ P(T^{s4\rightarrow s4}=2) = 1 \Rightarrow P(T^{s4\rightarrow s4}=+\infty)=0,\ P(T^{s5\rightarrow s5}=2) = 1 \Rightarrow P(T^{s5\rightarrow s5}=+\infty)=0,\ P(T^{s6\rightarrow s6}=1) = 1 \Rightarrow P(T^{s6\rightarrow s6}=+\infty)=0. \end{array}$$ So the states $s4, s5, s6$ are recurrent (since $E(T^{s4\rightarrow s4})=2=E(T^{s5\rightarrow s5})$ and $E(T^{s6\rightarrow s6})=1$, they are positive recurrent), whereas the state $s1$ is transient. Furthermore, for the last two states, we have $$\begin{array}{l} P(T^{s2\rightarrow s2}=2k) = 0.8^k(1-p)^k, p>0 \Rightarrow P(T^{s2\rightarrow s2}=+\infty)=0, \ P(T^{s3\rightarrow s3}=2k) = (1-p)^k0.8^k, p>0 \Rightarrow P(T^{s3\rightarrow s3}=+\infty)=0. \ \end{array}$$ So the states $s2, s3$ are recurrent (since $E(T^{s2\rightarrow s2})=E(T^{s3\rightarrow s3})=1/(1-0.8(1-p))<+\infty$, they are positive recurrent).

## Analysis of the eigenvalues and eigenvectors

Let us compute the characteristic polynom of the matrix $M$ \begin{eqnarray} \chi_M(x) &=& \left|\begin{matrix} -x & 1/2 & 1/2 & 0 & 0 & 0 \ 0 & -x & 1-p & 0 & 0 & p \ 0 & 1-q & -x & q & 0 & 0 \ 0 & 0 & 0 & -x & 1 & 0 \ 0 & 0 & 0 & 1 & -x & 0 \ 0 & 0 & 0 & 0 & 0 & 1-x \ \end{matrix}\right| = \left|\begin{matrix} -x & 1/2 & 1/2 & 0 & 0 \ 0 & -x & 1-p & 0 & 0 \ 0 & 1-q & -x & q & 0 \ 0 & 0 & 0 & -x & 1 \ 0 & 0 & 0 & 1 & -x \ \end{matrix}\right| (1-x) \ &=& -(1-x) \left|\begin{matrix} -x & 1/2 & 1/2 & 0 \ 0 & -x & 1-p & 0 \ 0 & 1-q & -x & 0 \ 0 & 0 & 0 & 1 \ \end{matrix}\right| -x(1-x) \left|\begin{matrix} -x & 1/2 & 1/2 & 0 \ 0 & -x & 1-p & 0 \ 0 & 1-q & -x & q \ 0 & 0 & 0 & -x \ \end{matrix}\right| \ &=& -(1-x) \left|\begin{matrix} -x & 1/2 & 1/2 \ 0 & -x & 1-p \ 0 & 1-q & -x \ \end{matrix}\right| +x^2(1-x) \left|\begin{matrix} -x & 1/2 & 1/2 \ 0 & -x & 1-p \ 0 & 1-q & -x \ \end{matrix}\right| \ &=& +x(1-x) \left|\begin{matrix} -x & 1-p \ 1-q & -x \ \end{matrix}\right| -x^3(1-x) \left|\begin{matrix} -x & 1-p \ 1-q & -x \ \end{matrix}\right| \ &=& x(1-x) -x^3(1-x) = x(1-x)^2(1+x)(x-(1-p)(1-q))(x+(1-p)(1-q)) \end{eqnarray} Therefore the eigenvalues are $\pm 1, 0, \pm(1-p)(1-q)$ with 1 of multiplicity 2. Hence the power of $M$ is $$M^n = P \left|\begin{matrix} 1^n & 0 & 0 & 0 & 0 & 0 \ 0 & 1^n & 0 & 0 & 0 & 0 \ 0 & 0 & (-1)^n & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & (1-p)^n(1-q)^n & 0 \ 0 & 0 & 0 & 0 & 0 & (1-p)^n(1-q)^n \ \end{matrix}\right| P^{-1}$$ with an appropriate matrix $P$ changing base in $\mathbb R^6$.

## Numerical consideration

The classification of states is given by dedicated functions

communicatingClasses(chain1)
recurrentClasses(chain1)
absorbingStates(chain1)
transientStates(chain1)


Let us compute power of $M$. Note that even powers of $M$ yield to a matrix different than odd powers since $P(T^{s4\rightarrow s4}=2)=P(T^{s5\rightarrow s5}=2)=1$.

library(expm)
M %^% 1000
M %^% 1001


Invariant measures are the eigenvectors associated to the eigen value 1 for the transposed matrix $M^T$ such that $\mu_i\geq 0$. If in addition $\sum_i \mu_i=1$, then it is a probability measure.

Let us compute the eigenvalue. We don't retrieve the theoretical values $\pm 1, 0, \pm(1-p)(1-q)$.

c(1, -1, 0, .99*.8, -.99*.8)
eigen(t(M))


Keeping those equalling 1, we obtain

eigen(t(M))$vectors[,eigen(t(M))$values >= 1]


So numerically, only one (probability) measure is obtained $(0,0,0,0,0,1)$. That's rather logical.

One can get the "boundary" measures by looking to the subset of recurrent states, extracting positive eigenvectors associated to the eigenvalue 1 and normalize ti. Indeed

Msub <- M[rownames(M) %in% unlist(recurrentClasses(chain1)), colnames(M) %in% unlist(recurrentClasses(chain1))]
eigen(Msub)$vectors[,eigen(Msub)$values >= 1]
t(eigen(Msub)$vectors[,eigen(Msub)$values >= 1]) / colSums(eigen(Msub)$vectors[,eigen(Msub)$values >= 1])


The steadyStates function in [@pkg:markovchain] implements the abovementioned algorithm when negative values are found in the eigenvectors of the full matrix corresponding to unitary eigenvalues.

## Try the markovchain package in your browser

Any scripts or data that you put into this service are public.

markovchain documentation built on Aug. 24, 2018, 1:03 a.m.