50. Markov chains

Everything so far has treated observations as independent, or at most jointly distributed at a single moment. A stochastic process is a collection of random variables indexed by time — a system that evolves.

The simplest useful kind has a very short memory.

The Markov property

A process X_0, X_1, X_2, \dots is a Markov chain if

P(X_{n+1} = j \mid X_n = i,\; X_{n-1}, \dots, X_0) = P(X_{n+1} = j \mid X_n = i)

The future depends on the present, and on nothing before it. Given where you are now, how you got here is irrelevant.

This is a strong assumption, and it's what makes the mathematics tractable. It's also less restrictive than it first appears: any finite-memory process can be made Markov by enlarging the state to include the recent history. "Yesterday and today's weather" is a state that makes a two-day-memory process Markovian.

The transition matrix

For a chain with states 1, \dots, k, collect the one-step probabilities:

P_{ij} = P(X_{n+1} = j \mid X_n = i)

\mathbf{P} is a stochastic matrix: entries non-negative, and each row sums to 1 (from state i you must go somewhere).

Rows are "from", columns are "to". Getting that backwards is the most common implementation bug in this topic — the other convention exists, and then the columns sum to 1 and the matrix multiplications transpose.

Multi-step transitions

Two steps means summing over the intermediate state:

P(X_2 = j \mid X_0 = i) = \sum_m P_{im}P_{mj} = (\mathbf{P}^2)_{ij}

That's just matrix multiplication — and the law of total probability from §0, partitioning on where you were in between. In general:

P(X_n = j \mid X_0 = i) = (\mathbf{P}^n)_{ij}

Matrix powers give n-step transitions. If the distribution over states at time 0 is a row vector \boldsymbol\pi_0:

\boldsymbol\pi_n = \boldsymbol\pi_0\mathbf{P}^n

Stationary distributions

A distribution \boldsymbol\pi is stationary if it doesn't change:

\boldsymbol\pi\mathbf{P} = \boldsymbol\pi, \qquad \sum_i \pi_i = 1

\boldsymbol\pi is a left eigenvector of \mathbf{P} with eigenvalue 1. Every stochastic matrix has one, because the rows summing to 1 guarantees eigenvalue 1 exists.

Once the chain reaches \boldsymbol\pi it stays there — not frozen in one state, but with the proportions across states holding steady.

When does the chain converge to it?

Two conditions:

Irreducible — every state is reachable from every other. No isolated clusters.

Aperiodic — the chain isn't trapped in a fixed cycle. Formally, the gcd of the return times to a state is 1. A chain that alternates strictly A→B→A→B has period 2 and never settles, even though a stationary distribution exists.

Given both, for a finite chain:

\lim_{n\to\infty}(\mathbf{P}^n)_{ij} = \pi_j \quad \text{for every } i

Every row of \mathbf{P}^n converges to \boldsymbol\pi. The starting state is forgotten entirely — the chain's long-run behaviour is a property of the transition rules alone.

There's also a useful interpretation: \pi_j is the long-run fraction of time spent in state j, and 1/\pi_j is the mean return time to j.

Detailed balance

A sufficient (not necessary) condition for \boldsymbol\pi to be stationary:

\pi_iP_{ij} = \pi_jP_{ji} \quad \text{for all } i, j

The flow from i to j equals the flow back. A chain satisfying this is reversible — run it backwards and it looks statistically identical.

Detailed balance is much easier to check than $\boldsymbol\pi\mathbf P = \boldsymbol\pi$, and it's the design principle behind MCMC: construct a chain whose stationary distribution is the posterior you want to sample from, then run it. That's how modern Bayesian computation works, and §11's sampling material is the foundation for it.

Worked example

Weather: sunny (S) or rainy (R).

• Sunny today → 80% sunny tomorrow.
• Rainy today → 40% sunny tomorrow.

\mathbf{P} = \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix}

Two days ahead, starting sunny:

\mathbf{P}^2 = \begin{pmatrix} 0.72 & 0.28 \\ 0.56 & 0.44 \end{pmatrix}

So P(\text{sunny in 2 days} \mid \text{sunny today}) = 0.72 — down from 0.8, drifting toward the long-run value.

The stationary distribution. Solve $\boldsymbol\pi\mathbf{P} = \boldsymbol\pi$ with \pi_S + \pi_R = 1:

\pi_S = 0.8\pi_S + 0.4\pi_R \implies 0.2\pi_S = 0.4\pi_R \implies \pi_S = 2\pi_R

With \pi_S + \pi_R = 1:

\boldsymbol\pi = \left(\tfrac{2}{3}, \tfrac{1}{3}\right)

Two-thirds of days are sunny in the long run, whatever the weather today. And the mean time between rainy days is 1/\pi_R = 3 days.

Check detailed balance: \pi_S P_{SR} = \frac23(0.2) = 0.1333 and \pi_R P_{RS} = \frac13(0.4) = 0.1333. ✓ This chain is reversible.

Doing it in Python

Transitions, powers and convergence:

import numpy as np

P = np.array([[0.8, 0.2],
              [0.4, 0.6]])
states = ["sunny", "rainy"]

print("row sums (must all be 1):", P.sum(axis=1))

print(f"\n{'n':>4} {'P^n row 0':>22} {'P^n row 1':>22}")
Pn = np.eye(2)
for n in range(0, 21):
    if n in (1, 2, 5, 10, 20):
        print(f"{n:>4} {str(Pn[0].round(5)):>22} {str(Pn[1].round(5)):>22}")
    Pn = Pn @ P

print("\nBoth rows converge to the same vector -- the start is forgotten.")

Finding the stationary distribution three ways:

import numpy as np

P = np.array([[0.8, 0.2],
              [0.4, 0.6]])

# 1. By matrix power
pi_power = np.linalg.matrix_power(P, 200)[0]

# 2. As a left eigenvector with eigenvalue 1
vals, vecs = np.linalg.eig(P.T)
v = np.real(vecs[:, np.argmin(np.abs(vals - 1))])
pi_eig = v / v.sum()

# 3. By solving the linear system directly
A = np.vstack([(P.T - np.eye(2))[:-1], np.ones(2)])
b = np.array([0, 1])
pi_solve = np.linalg.solve(A, b)

print("by matrix power :", pi_power.round(6))
print("by eigenvector  :", pi_eig.round(6))
print("by linear solve :", pi_solve.round(6))
print("\nexact: (2/3, 1/3) =", np.array([2/3, 1/3]).round(6))
print("stationary check pi @ P == pi:", np.allclose(pi_solve @ P, pi_solve))

Simulating the chain, and confirming the time-average interpretation:

import numpy as np

rng = np.random.default_rng(0)
P = np.array([[0.8, 0.2],
              [0.4, 0.6]])

def simulate(P, start, steps, rng):
    state = start
    visits = np.zeros(len(P), dtype=int)
    for _ in range(steps):
        visits[state] += 1
        state = rng.choice(len(P), p=P[state])
    return visits

for start in (0, 1):
    visits = simulate(P, start, 200_000, rng)
    print(f"starting {'sunny' if start == 0 else 'rainy'}: "
          f"time fractions {(visits / visits.sum()).round(4)}")

print("\ntheory: [0.6667 0.3333] -- and the starting state doesn't matter.")
print("mean return time to rainy =", round(1 / (1/3), 4), "days")

A chain that never converges, because it's periodic:

import numpy as np

# Strict alternation: A -> B -> A -> B ...
P = np.array([[0.0, 1.0],
              [1.0, 0.0]])

print("A stationary distribution still EXISTS:")
A = np.vstack([(P.T - np.eye(2))[:-1], np.ones(2)])
pi = np.linalg.solve(A, [0, 1])
print("  pi =", pi, " check:", np.allclose(pi @ P, pi))

print("\nBut P^n never settles -- it oscillates forever:")
for n in (1, 2, 3, 10, 11, 100, 101):
    print(f"  P^{n:<4} row 0 = {np.linalg.matrix_power(P, n)[0]}")

print("\nThis chain is irreducible but PERIODIC (period 2), so the")
print("convergence theorem does not apply.")

And detailed balance as a design tool — building a chain to target a chosen distribution, the core idea behind MCMC:

import numpy as np

rng = np.random.default_rng(1)

# The distribution we WANT to sample from
target = np.array([0.1, 0.2, 0.5, 0.2])
k = len(target)

# Metropolis-Hastings on a symmetric proposal: propose a neighbour, accept
# with probability min(1, target[new] / target[old]) -- this construction
# satisfies detailed balance by design.
def step(i, rng):
    j = rng.integers(k)                       # symmetric proposal
    return j if rng.random() < min(1, target[j] / target[i]) else i

state, visits = 0, np.zeros(k, dtype=int)
for _ in range(400_000):
    visits[state] += 1
    state = step(state, rng)

print("target   :", target)
print("simulated:", (visits / visits.sum()).round(4))
print("\nThe chain was never told the answer -- detailed balance made its")
print("stationary distribution equal the target.")

Your turn

1. A chain has \mathbf{P} = \begin{pmatrix}0.5 & 0.5\\ 0.2 & 0.8\end{pmatrix}. Find the stationary distribution.

2. Why must every row of a transition matrix sum to 1?

3. A chain alternates strictly between two states. Does a stationary distribution exist? Does the chain converge to it?

Solutions

1. Solve \boldsymbol\pi\mathbf P = \boldsymbol\pi:

\pi_1 = 0.5\pi_1 + 0.2\pi_2 \implies 0.5\pi_1 = 0.2\pi_2 \implies \pi_2 = 2.5\pi_1

With \pi_1 + \pi_2 = 1:

\pi_1(1 + 2.5) = 1 \implies \pi_1 = \frac{1}{3.5} = \frac{2}{7}, \qquad \pi_2 = \frac{5}{7}

So \boldsymbol\pi = (2/7, 5/7) \approx (0.286, 0.714).

Sanity check: state 2 is "stickier" (0.8 chance of staying versus 0.5), so it should hold more of the long-run mass. It does.

2. Because from state i the chain must go somewhere at the next step. Row i lists the probabilities of every possible destination, and those destinations are exhaustive and mutually exclusive:

\sum_j P_{ij} = \sum_j P(X_{n+1} = j \mid X_n = i) = 1

That's Axiom 2 applied to the conditional distribution given X_n = i (§0).

A practical consequence: a row that doesn't sum to 1 is a bug, and it's the first thing to check when a simulation misbehaves. It also guarantees the vector of all ones is a right eigenvector with eigenvalue 1 — which is why a stationary distribution (a left eigenvector for the same eigenvalue) always exists.

3. A stationary distribution exists; the chain does not converge to it.

For \mathbf P = \begin{pmatrix}0&1\\1&0\end{pmatrix}, solving \boldsymbol\pi\mathbf P = \boldsymbol\pi gives \boldsymbol\pi = (0.5, 0.5), which is genuinely stationary — start there and you stay there.

But starting in state A, the chain is in A at every even time and B at every odd time, forever. \mathbf P^n alternates between the identity and the swap matrix and never approaches anything.

The chain is irreducible but periodic (period 2), so the convergence theorem doesn't apply. The time-average interpretation survives — the chain does spend half its time in each state — but the distribution at time n never settles.

The standard fix is to add a small self-loop (a "lazy" chain): with probability \varepsilon stay put. That breaks the period, makes the chain aperiodic, and convergence follows, with the same stationary distribution.

Check yourself in code

Find the stationary distribution of the weather chain and confirm that P^n converges to it from either starting state.

Print exactly this:

stationary [0.6667 0.3333]
is stationary: True
P^50 row 0 [0.6667 0.3333]
P^50 row 1 [0.6667 0.3333]
rows agree: True

Round every vector to 4 decimal places. Solve for the stationary distribution as a linear system rather than by taking a power.

import numpy as np

P = np.array([[0.8, 0.2],
              [0.4, 0.6]])

A = np.vstack([(P.T - np.eye(2))[:-1], np.ones(2)])
pi = np.linalg.solve(A, [0, 1])
print("stationary", pi.round(4))

# Confirm pi @ P == pi, then show both rows of P^50 match it.

import numpy as np

P = np.array([[0.8, 0.2],
              [0.4, 0.6]])

A = np.vstack([(P.T - np.eye(2))[:-1], np.ones(2)])
pi = np.linalg.solve(A, [0, 1])
print("stationary", pi.round(4))
print("is stationary:", bool(np.allclose(pi @ P, pi)))

P50 = np.linalg.matrix_power(P, 50)
print("P^50 row 0", P50[0].round(4))
print("P^50 row 1", P50[1].round(4))
print("rows agree:", bool(np.allclose(P50[0], P50[1])))

---

A Markov chain forgets everything except its current state. One-step probabilities go in a transition matrix whose rows sum to 1; n-step probabilities are matrix powers. If the chain is irreducible and aperiodic it converges to a unique stationary distribution regardless of where it started — and detailed balance lets you design a chain whose stationary distribution is whatever you want, which is the entire basis of MCMC.

Next: a process where events arrive continuously in time.