← Home

Find the Steady State Distribution of a Markov Process in R.

A Markov process consists of states and probabilities, where the probability of transitioning from one state to another depends only on the current state and not on the past; it is memoryless. A Markov process is often depicted as a transition probability matrix, where the $(i, j)$ entry of this matrix is the probability that the Markov process transitions to state $j$ at the next time-step given that the process is currently in state $i$.

Consider the following Markov process, where each state is my mood on a given day (assume these are mutually exclusive).

p = matrix(c(0.4,0.0,0.1,
            nrow = 3,
            ncol = 3)

states = c("Happy", "Gloomy", "Sleepy")
colnames(p) = states
rownames(p) = states

##        Happy Gloomy Sleepy
## Happy    0.4    0.4    0.2
## Gloomy   0.0    0.5    0.5
## Sleepy   0.1    0.3    0.6

$P$ is the one-time-step transition probability matrix (in this example, a time-step is a day), where the $(i, j)$ entry indicates the probability that my mood one day from now will be in state $j$ given that my mood is currently in state $i$. For example, if I am gloomy today, then tomorrow there is a 50% chance that I will be gloomy and a 50% chance that I will be sleepy.

To get the two-time-step transition probability matrix, you raise the one-time-step transition probability matrix, $P$, to the second power ($P^2$).


p %^% 2
##        Happy Gloomy Sleepy
## Happy   0.18   0.42   0.40
## Gloomy  0.05   0.40   0.55
## Sleepy  0.10   0.37   0.53

$P^2$ is the two-time-step transition probability matrix. The $(i, j)$ entry indicates the probability that my mood two days from now will be in state $j$ given that my mood is currently in state $i$. For example, if I am happy today, then two days from now there is an 18% chance that I will be happy, a 42% chance that I will be gloomy, and a 40% chance that I will be sleepy.

To get the $k$-time-step transition probability matrix, you raise the one-time-step probability matrix, $P$, to the $k^{th}$ power ($P^k$).

As you continue raising $P$ to higher powers, you’ll start to notice that the probabilities are converging.

p %^% 5
##          Happy  Gloomy  Sleepy
## Happy  0.08886 0.38766 0.52348
## Gloomy 0.08675 0.38530 0.52795
## Sleepy 0.08824 0.38617 0.52559
p %^% 20
##            Happy    Gloomy    Sleepy
## Happy  0.0877193 0.3859649 0.5263158
## Gloomy 0.0877193 0.3859649 0.5263158
## Sleepy 0.0877193 0.3859649 0.5263158

For Markov processes that are finite, irreducible, and aperiodic, there is a long-run equilibrium that is reached despite the starting state. For example, when looking at $P^{20}$, the chance that I will be sleepy 20 days from now is 52.63% regardless of what my mood is today.

This equilibrium is called the steady-state distribution of the Markov process and is typically denoted by $\mu$ or $\pi$. I prefer $\mu$. The steady-state distribution can be interpreted as the fraction of time the Markov process spends in state $i$ in the long-run. In this example, after raising $P$ to a high enough power, we can see that the steady-state distribution is

Instead of raising $P$ to sufficiently high powers, this steady-state distribution can be found more easily in a few ways, all of which come down to solving

for $\mu$ with the constraint that $\sum_{i} \mu_{i} = 1$.

Solve with System of Equations

$ \mu P = \mu $ can be written as a system of equations.

with the constraint that

After some re-arranging, we can get this into the form $Ax = b$ where we can solve for $x$ with R’s qr.solve.

It should be noted that $A$ is equivalent to $(P - I)^T$ (where $I$ is the identity matrix) with a row of 1s appended to the bottom.

In R, we can solve for $x$ with the following:

n = ncol(p)
A = t(p - diag(n))
A = rbind(A, rep(1, n))
b = c(rep(0, n), 1)

##        Happy Gloomy Sleepy
## Happy   -0.6    0.0    0.1
## Gloomy   0.4   -0.5    0.3
## Sleepy   0.2    0.5   -0.4
##          1.0    1.0    1.0
## [1] 0 0 0 1
mu = qr.solve(A, b)

##     Happy    Gloomy    Sleepy 
## 0.0877193 0.3859649 0.5263158

Solve with Eigenvectors

A row vector that remains unchanged (or is a constant multiple of itself) when multiplied by a square matrix is a left-eigenvector. That is, a vector $v$ that satisfies the equation

is a left-eigenvector of $A$ where $\lambda$ are the eigenvalues.

Because $\mu$ is a solution to $ \mu P = \mu $, $\mu$ is a left-eigenvector of the square matrix $P$ corresponding to the eigenvalue $\lambda = 1$. In R, the left-eigenvectors can be found by executing eigen on the transpose of $P$.

e = eigen(t(p))
first = Re(e$vectors[ ,1])

# Normalize.
mu = first / sum(first)

## [1] 0.0877193 0.3859649 0.5263158

Solve with Null Spaces

My favorite way to get the steady-state distribution is by finding the basis of the null space of $(P - I)^T$, where $I$ is the identity matrix. Note that the null space of $(P - I)^T$ is equivalent to the left null space of $P - I$.

Rearranging $\mu P = \mu$ yields

such that

This might have been clear when you re-arranged the system of equations and set everything to 0.

Conveniently, the pracma package has a function nullspace for finding the basis of a null space.


A = t(p - diag(n))
basis = nullspace(A)

# Normalize.
mu = basis / sum(basis)

##           [,1]
## [1,] 0.0877193
## [2,] 0.3859649
## [3,] 0.5263158
comments powered by Disqus