Pre-requisites

This document assumes basic familiarity with Markov chains and linear algebra.

Overview

In this note, we illustrate one way of analytically obtaining the stationary distribution for a finite discrete Markov chain.

3x3 example

Assume our probability transition matrix is: $$P = \begin{bmatrix} 0.7 & 0.2 & 0.1 \\ 0.4 & 0.6 & 0 \\ 0 & 1 & 0 \end{bmatrix}$$

Since every state is accessible from every other state, this Markov chain is irreducible. Every irreducible finite state space Markov chain has a unique stationary distribution. Recall that the stationary distribution $\pi$ is the vector such that $$\pi = \pi P$$ .

Therefore, we can find our stationary distribution by solving the following linear system: $$\begin{align*} 0.7\pi_1 + 0.4\pi_2 &= \pi_1 \\ 0.2\pi_1 + 0.6\pi_2 + \pi_3 &= \pi_2 \\ 0.1\pi_1 &= \pi_3 \end{align*}$$ subject to $\pi_1 + \pi_2 + \pi_3 = 1$. Putting these four equations together and moving all of the variables to the left hand side, we get the following linear system: $$\begin{align*} -0.3\pi_1 + 0.4\pi_2 &= 0 \\ 0.2\pi_1 + -0.4\pi_2 + \pi_3 &= 0 \\ 0.1\pi_1 - \pi_3 &= 0 \\ \pi_1 + \pi_2 + \pi_3 &= 1 \end{align*}$$

We will define the linear system in matrix notation: $$\underbrace{\begin{bmatrix} -0.3 & 0.4 & 0 \\ 0.2 & -0.4 & 1 \\ 0.1 & 0 & 1 \\ 1 & 1 & 1 \end{bmatrix}}_A \begin{bmatrix} \pi_1 \\ \pi_2 \\ \pi_3 \end{bmatrix} = \underbrace{\begin{bmatrix} 0 \\ 0 \\ 0 \\ 1 \end{bmatrix}}_b \\ A\pi = b$$

Since this linear system has more equations than unknowns, it is an overdeterminted system. Recall from linear algebra that an overdetermined system is consistent (i.e. we can solve for $\pi$ exactly) when $b$ is in the column space of $A$. We can check this numerically by obtaining the rank of $A$, then obtaining the rank of an augmented matrix with $b$ appended as a column of $A$. (Side note: it can be difficult to find the rank of an arbitrary matrix numerically, espeically for large matrices. But for our small example, we are safe.)

library(Matrix)
A        <- matrix(c(-0.3, 0.2, 0.1, 1, 0.4, -0.4, 0, 1, 0, 1, -1, 1 ), ncol=3,nrow=4)
b        <- c(0,0,0, 1)
rank.A   <- as.numeric(rankMatrix(A))
rank.Ab  <- as.numeric(rankMatrix(cbind(A,b)))
print(paste("The rank of A =", rank.A, "and the rank of the augmented matrix =", rank.Ab))

[1] "The rank of A = 3 and the rank of the augmented matrix = 3"

We see that $A$ has full column rank, and that the rank is unchanged when we add $b$ as a column. Therefore, $b$ is in the column space of $A$, and this system is consistent. We can find $\pi$ by solving the normal equations: $$A^TA\pi = A^Tb$$

We use the solve function in R to solve for the stationary distribution $\pi$:

pi        <- drop(solve(t(A) %*% A, t(A) %*% b))
names(pi) <- c('state.1', 'state.2', 'state.3')
pi 

   state.1    state.2    state.3 
0.54054054 0.40540541 0.05405405 

We find that: $$\begin{align*} \pi_1 \approx 0.54, \pi_2 \approx 0.41, \pi_3 \approx 0.05 \end{align*}$$

Therefore, under proper conditions, we expect the Markov chain to spend more time in states 1 and 2 as the chain evolves.

Session information

sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.5 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] knitr_1.14     MASS_7.3-45    expm_0.999-0   Matrix_1.2-7.1
[5] rmarkdown_1.1 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.7     lattice_0.20-34 gtools_3.5.0    digest_0.6.9   
 [5] assertthat_0.1  mime_0.4        R6_2.1.2        grid_3.3.2     
 [9] xtable_1.8-2    formatR_1.4     magrittr_1.5    evaluate_0.9   
[13] stringi_1.1.1   tools_3.3.2     stringr_1.0.0   shiny_0.13.2   
[17] httpuv_1.3.3    yaml_2.1.13     htmltools_0.3.5 tibble_1.2