Stationary Distributions in Markov Models

What It Means to Be Stationary

A Markov model describes how a probability distribution over states changes from one step to the next. A stationary distribution is a distribution that does not change under this update. Finding one means solving a fixed-point equation, while interpreting one requires care: a distribution can be stationary without every starting distribution converging to it.

We use the column-vector convention for a finite model with \(n\) states. A state distribution is a column \(\mathbf{x}\in\mathbb{R}^n\) whose entries are nonnegative and sum to \(1\). The transition matrix \(\mathbf{P}\in\mathbb{R}^{n\times n}\) is column-stochastic: its entries are nonnegative, and the entries in each column sum to \(1\). Thus \(p_{ij}\) is the probability of moving to state \(i\) from state \(j\). If the current distribution is \(\mathbf{x}\), then the next distribution is \(\mathbf{P}\mathbf{x}\).

The equation \(\mathbf{P}\boldsymbol{\pi}=\boldsymbol{\pi}\) is homogeneous when written as \((\mathbf{P}-\mathbf{I}_n)\boldsymbol{\pi}=\mathbf{0}\). That equation alone is not enough: it admits the zero vector, which is not a probability distribution, and any nonzero solution can be rescaled. The sum-to-one constraint selects the required scale, while nonnegativity must also be checked.

Proof. If \(\boldsymbol{\pi}\) is stationary, its defining conditions give all three stated conditions. Conversely, if a vector satisfies them, then \((\mathbf{P}-\mathbf{I}_n)\boldsymbol{\pi}=\mathbf{0}\) rearranges to \(\mathbf{P}\boldsymbol{\pi}=\boldsymbol{\pi}\), and the other two conditions say it is a probability distribution. Hence it is stationary. \(\square\)

Setting Up the Linear System

For each state, the corresponding equation in \(\mathbf{P}\boldsymbol{\pi}=\boldsymbol{\pi}\) balances the probability flowing into that state against the probability already assigned to it. In components, the \(i\)th equation is \(\sum_{j=1}^n p_{ij}\pi_j=\pi_i\). These equations are not all independent. Since the columns of \(\mathbf{P}\) sum to \(1\), the entries of \(\mathbf{P}\boldsymbol{\pi}\) sum to the entries of \(\boldsymbol{\pi}\). Consequently, summing all \(n\) equations in \((\mathbf{P}-\mathbf{I}_n)\boldsymbol{\pi}=\mathbf{0}\) gives \(0=0\). A system may therefore have a redundant equation; this does not remove the need to impose \(\sum_i\pi_i=1\).

Solving a Three-State System

For a larger system, the same equations can be organized as one linear system. The matrix \(\mathbf{P}-\mathbf{I}_n\) has rows that sum to zero as a collection: because \(\mathbf{P}\) is column-stochastic, \(\mathbf{1}^{\top}(\mathbf{P}-\mathbf{I}_n)=\mathbf{0}^{\top}\), where \(\mathbf{1}\) is the column of \(n\) ones. Thus the sum of the rows of \(\mathbf{P}-\mathbf{I}_n\) is the zero row. One fixed-point equation is redundant. Keep a sufficient set of the fixed-point equations, add the normalization equation, solve, and then check every original equation and the probability conditions.

Existence and Long-Run Interpretation

Proof. Start with any probability distribution \(\mathbf{x}_0\), for example a standard basis vector. Define \(\mathbf{x}_{k+1}=\mathbf{P}\mathbf{x}_k\). Each \(\mathbf{x}_k\) remains a probability distribution: nonnegativity follows from multiplying a nonnegative matrix by a nonnegative vector, and the sum of the entries is preserved because \(\mathbf{1}^{\top}\mathbf{P}=\mathbf{1}^{\top}\). For each positive integer \(N\), form the average \(\mathbf{q}_N=\frac1N\sum_{k=0}^{N-1}\mathbf{x}_k\). Its entries are nonnegative and sum to \(1\), so the sequence \((\mathbf{q}_N)\) is bounded in \(\mathbb{R}^n\). By the Bolzano–Weierstrass theorem, a bounded sequence in finite-dimensional real space has a convergent subsequence; write its limit as \(\boldsymbol{\pi}\). Limits preserve nonnegativity and the finite sum, so \(\boldsymbol{\pi}\) is a probability distribution. Also, by cancellation of consecutive terms, \(\mathbf{P}\mathbf{q}_N-\mathbf{q}_N=(\mathbf{x}_N-\mathbf{x}_0)/N\). Each coordinate of \(\mathbf{x}_N-\mathbf{x}_0\) lies between \(-1\) and \(1\), so the right-hand side tends to \(\mathbf{0}\). Along the convergent subsequence, matrix multiplication by the fixed matrix \(\mathbf{P}\) preserves limits. Taking limits gives \(\mathbf{P}\boldsymbol{\pi}-\boldsymbol{\pi}=\mathbf{0}\). Therefore \(\boldsymbol{\pi}\) is stationary. \(\square\)

This existence result does not say that the stationary distribution is unique, or that the distributions \(\mathbf{x}_k\) from every starting state converge to one. A stationary distribution is invariant: if \(\mathbf{x}_0=\boldsymbol{\pi}\), then \(\mathbf{x}_1=\mathbf{P}\boldsymbol{\pi}=\boldsymbol{\pi}\), and repeated updates leave it unchanged. More generally, if the sequence \(\mathbf{x}_k\) does converge to a limit \(\mathbf{x}\), then \(\mathbf{x}\) must be stationary: \(\mathbf{x}_{k+1}=\mathbf{P}\mathbf{x}_k\), and taking limits gives \(\mathbf{x}=\mathbf{P}\mathbf{x}\). The reverse implication does not follow just from stationarity.

Common Pitfall: Mixing Up Row and Column Conventions

With the convention used here, state distributions are columns and are updated by \(\mathbf{x}_{k+1}=\mathbf{P}\mathbf{x}_k\). A common error is to solve \(\mathbf{P}\mathbf{q}=\mathbf{q}\) using a row vector, or to multiply a row vector on the right without changing the fixed-point equation. These choices describe a different convention and can produce the wrong answer.

For the three-state matrix above, solving \(\mathbf{q}\mathbf{P}=\mathbf{q}\) as though distributions were rows gives \(\mathbf{q}=(\tfrac13,\tfrac13,\tfrac13)\). This is not stationary under the column convention: multiplying the column \((\tfrac13,\tfrac13,\tfrac13)^\top\) by \(\mathbf{P}\) gives \((\tfrac12,\tfrac13,\tfrac16)^\top\), not \((\tfrac13,\tfrac13,\tfrac13)^\top\). The correction is to use the state-vector convention consistently and solve \(\mathbf{P}\boldsymbol{\pi}=\boldsymbol{\pi}\), which gives \((\tfrac25,\tfrac25,\tfrac15)^\top\).

Another error is to stop after finding a nonzero solution of the homogeneous equations. For instance, a scalar multiple of a stationary vector still solves \((\mathbf{P}-\mathbf{I}_n)\boldsymbol{\pi}=\mathbf{0}\), but it is a probability distribution only when its entries are nonnegative and sum to \(1\). Always impose the normalization condition and check the result in the original fixed-point equation.