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
steadyStates(chain1)


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.