1 Markov Chains
1.1 Markov Chains
Markov property
For all n = 0, 1, 2, \ldots, and all states s^{k} \in \mathcal{S}, k = 0, 1, \ldots, n, we have
\mathbb{P}(X_{n} = s^{n} | X_{n-1} = s^{n-1}, \ldots, X_0 = s^{0}) = \mathbb{P}(X_{n} = s^{n} | X_{n-1} = s^{n-1})
one-step transition probability
p_{ij} := \mathbb{P}(X_{n+1} = j | X_n = i)
n-step transition probability
p^{(n)}_{ij} := \mathbb{P}(X_{k+n} = j | X_k = i)
Since p^{(n)}_{ij} are conditional probabilities, for all i, j \in \mathcal{S} and n = 0, 1, 2, \ldots, they satisfy the properties
p^{(n)}_{ij} \geq 0, \quad \sum_{j \in \mathcal{S}} p^{(n)}_{ij} = 1.
Let \mathbf{P} be the transition matrix of a Markov chain, that is, \mathbf{P} is a |\mathcal{S}| \times |\mathcal{S}| matrix with (i, j)-th entry p_{ij}. We use \mathbf{P}^{(n)} to denote the transition matrix of the n-step transition probabilities. Then, \mathbf{P}^{(0)} = \mathbf{I}, the identity matrix, and \mathbf{P}^{(n)} have entries p^{(n)}_{ij}, for all i, j \in \mathcal{S}.
Let m = |\mathcal{S}| - 1. Then, we can write \mathbf{P}^{(n)} as
\mathbf{P}^{(n)} = \begin{pmatrix} p^{(n)}_{00} & p^{(n)}_{01} & \cdots & p^{(n)}_{0m} \\ p^{(n)}_{10} & p^{(n)}_{11} & \cdots & p^{(n)}_{1m} \\ \vdots & \vdots & \ddots & \vdots \\ p^{(n)}_{m0} & p^{(n)}_{m1} & \cdots & p^{(n)}_{mm} \end{pmatrix}.
Given the initial distribution \rho_0(s) := \rho(s) := \mathbb{P}(X_0 = s) for s \in \mathcal{S}, the distribution of X_n is given by
\rho^n(s) = \sum_{j \in \mathcal{S}} \rho(j) p^{(n)}_{sj}.
We can write the above equation in matrix form as
\boldsymbol{\rho}^{(n)} = \boldsymbol{\rho} \mathbf{P}^{(n)}.
Note that \boldsymbol{\rho}^{(n)} is a row vector, and \boldsymbol{\rho}^{(0)} = \boldsymbol{\rho}.
1.2 Chapman-Kolmogorov equations
The n-step transition probabilities can be computed recursively using the Chapman-Kolmogorov equations:
p^{(n)}_{ij} = \sum_{k \in \mathcal{S}} p^{(m)}_{ik} p^{(n-m)}_{kj}, \quad \text{for all } i, j \in \mathcal{S}, \text{ and } m = 0, 1, \ldots, n.
Here, we introduce some special cases of the Chapman-Kolmogorov equations. For n = 2, we have
p^{(2)}_{ij} = \sum_{k \in \mathcal{S}} p_{ik} p_{kj}, \quad \text{for all } i, j \in \mathcal{S}.
We use matrix \mathbf{P}^2 to denote the two-step transition probabilities. This matrix can be computed as the product of the transition matrix with itself:
\mathbf{P}^{(2)} = \mathbf{P} \mathbf{P}.
Some special cases of the Chapman-Kolmogorov equations are
\mathbf{P}^{(n)} = \mathbf{P}^{n-1} \mathbf{P} = \mathbf{P} \mathbf{P}^{n-1} = \mathbf{P}^n, \quad \text{for all } n = 1, 2, \ldots
As a result, we can compute the n-step transition probabilities by multiplying the transition matrix with itself n times.
1.3 Classification of states
State j is said to be accessible from state i if there exists n \geq 0 such that p^{(n)}_{ij} > 0.
If state j is accessible from state i and state i is accessible from state j, then states i and j are said to communicate with each other.
Let A(i) denote the set of states that are accessible from state i. Then, a state i is said to be recurrent if for all j \in A(i), we have i \in A(j).
If a state i is not recurrent, then it is said to be transient. That is, there are some states j \in A(i) such that i \notin A(j).
If p_{ii} = 1, then state i is said to be absorbing. That means that once the Markov chain enters state i, it will stay in state i forever.
Example 1.1 (A chain with all three kinds of states) Consider a Markov chain on \mathcal{S} = \{1, 2, 3, 4, 5\} with transition matrix
\mathbf{P} = \begin{pmatrix} 1/2 & 1/2 & 0 & 0 & 0 \\ 1/4 & 0 & 1/4 & 0 & 1/2 \\ 0 & 0 & 1/2 & 1/2 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix},
whose transition diagram is shown in Figure 1.1. An arrow from i to j is drawn whenever p_{ij} > 0, and it is labelled with p_{ij}.
Reading off accessibility. State j is accessible from state i exactly when the diagram contains a directed path from i to j. For example, state 3 is accessible from state 1: the path 1 \to 2 \to 3 gives
p^{(2)}_{13} = \sum_{k \in \mathcal{S}} p_{1k} p_{k3} = \tfrac{1}{2} \times 0 + \tfrac{1}{2} \times \tfrac{1}{4} = \tfrac{1}{8} > 0.
Squaring the whole matrix makes the pattern of zeros visible:
\mathbf{P}^{(2)} = \mathbf{P}^2 = \begin{pmatrix} 3/8 & 1/4 & 1/8 & 0 & 1/4 \\ 1/8 & 1/8 & 1/8 & 1/8 & 1/2 \\ 0 & 0 & 3/4 & 1/4 & 0 \\ 0 & 0 & 1/2 & 1/2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix}.
A zero entry of \mathbf{P}^{(2)} does not mean the state is unreachable; accessibility asks for some n \geq 0. Indeed p^{(2)}_{14} = 0, but going one step further along 1 \to 2 \to 3 \to 4 gives p^{(3)}_{14} = \tfrac{1}{8} \times \tfrac{1}{2} = \tfrac{1}{16} > 0. Recall also that \mathbf{P}^{(0)} = \mathbf{I}, so every state is accessible from itself and i \in A(i) always holds. Collecting all directed paths,
\begin{aligned} A(1) &= \{1, 2, 3, 4, 5\}, \\ A(2) &= \{1, 2, 3, 4, 5\}, \\ A(3) &= \{3, 4\}, \\ A(4) &= \{3, 4\}, \\ A(5) &= \{5\}. \end{aligned}
Classifying the states. We now apply the definitions to each state.
- States 3 and 4 are recurrent. We have \mathcal{S}(3) = \mathcal{S}(4) = \{3, 4\}, and both 3 and 4 belong to this set. So every state accessible from 3 can reach 3 again, and likewise for 4. Note that this is visible in the matrix: rows 3 and 4 put all their mass on \{3, 4\}, so once the chain enters \{3, 4\} it never leaves. States 3 and 4 also communicate with each other, since p_{34} > 0 and p_{43} > 0.
- State 5 is recurrent, and moreover absorbing. Here p_{55} = 1, so \mathcal{S}(5) = \{5\} and the chain stays at 5 forever. Absorbing states are a special case of recurrent states.
- States 1 and 2 are transient. State 5 is accessible from state 1, but 1 \notin \mathcal{S}(5) = \{5\}; similarly state 3 is accessible from state 1, but 1 \notin \mathcal{S}(3) = \{3, 4\}. The same argument applies to state 2. Once the chain leaves \{1, 2\} it can never come back, which is exactly the block of zeros in the lower-left corner of \mathbf{P}.
Why “transient”. The name refers to the fact that the chain spends only a finite amount of time in \{1, 2\}. Restricting \mathbf{P} to the transient states gives
\mathbf{Q} = \begin{pmatrix} 1/2 & 1/2 \\ 1/4 & 0 \end{pmatrix}, \qquad \mathbf{Q}^2 = \begin{pmatrix} 3/8 & 1/4 \\ 1/8 & 1/8 \end{pmatrix},
and the row sums of \mathbf{Q}^2 are 5/8 and 1/4. The row sums of \mathbf{Q}^n give the probability that the chain is still inside \{1, 2\} after n steps, so starting from either state,
\mathbb{P}(X_{2k} \in \{1, 2\} \mid X_0 = i) \leq \left(\tfrac{5}{8}\right)^{k} \xrightarrow{\;k \to \infty\;} 0, \quad i \in \{1, 2\}.
With probability one the chain eventually settles either into the recurrent pair \{3, 4\}, where it alternates and lingers forever, or into the absorbing state 5.
1.4 Classes
If state i is recurrent, then the set of states A(i) are all communicating with each other, and they form a recurrence class or class.
A Markov chain can be partitioned into classes, and a set of transient states.
1.4.1 Periodicity
The period of a state i is defined as
d(i) := \gcd\{n \geq 1 : p^{(n)}_{ii} > 0\}.
Here, \gcd is the greatest common divisor. If d(i) = 1, then state i is said to be aperiodic; otherwise, it is said to be periodic.
The periodicity is also a class property.
Example 1.2 (The espresso machine in the lab) The graduate students’ lounge has exactly one espresso machine, and every morning its condition is recorded as good (\mathsf{G}), wobbly (\mathsf{W}: it still pours, but the pressure is off and it is overdue for descaling), or broken (\mathsf{B}: out of order, waiting for the technician). After years of careful observation over morning coffee, the students concluded that tomorrow’s condition depends on the past only through today’s condition:
- from \mathsf{G}: stays good with probability 0.8, turns wobbly with probability 0.15, breaks down with probability 0.05;
- from \mathsf{W}: someone finally descales it and it is good again with probability 0.4, stays wobbly with probability 0.4, breaks down with probability 0.2;
- from \mathsf{B}: the technician arrives and restores it to good with probability 0.5, gets it running but wobbly with probability 0.3, and with probability 0.2 the spare part has not arrived and it is still broken.
The condition X_n on morning n is therefore a Markov chain on \mathcal{S} = \{\mathsf{G}, \mathsf{W}, \mathsf{B}\}. Ordering the states as (\mathsf{G}, \mathsf{W}, \mathsf{B}), the transition matrix is
\mathbf{P} = \begin{pmatrix} 0.8 & 0.15 & 0.05 \\ 0.4 & 0.4 & 0.2 \\ 0.5 & 0.3 & 0.2 \end{pmatrix},
where the (i, j)-th entry is p_{ij}, so each row is a probability distribution over tomorrow’s condition and sums to one.
On the first morning of the semester, call it day 0, nobody knows how the machine was left over the break. Past semesters suggest
\boldsymbol{\rho} = \begin{pmatrix} 0.6 & 0.3 & 0.1 \end{pmatrix},
that is, \mathbb{P}(X_0 = \mathsf{G}) = 0.6, \mathbb{P}(X_0 = \mathsf{W}) = 0.3, and \mathbb{P}(X_0 = \mathsf{B}) = 0.1. What is the chance of getting coffee on day 2?
Propagating the distribution one step at a time. For day 1 we compute \boldsymbol{\rho}^{(1)} = \boldsymbol{\rho}\mathbf{P} entry by entry, each entry being a sum over all the ways of arriving at that state:
\begin{aligned} \rho^{(1)}(\mathsf{G}) &= 0.6 \times 0.8 + 0.3 \times 0.4 + 0.1 \times 0.5 = 0.65, \\ \rho^{(1)}(\mathsf{W}) &= 0.6 \times 0.15 + 0.3 \times 0.4 + 0.1 \times 0.3 = 0.24, \\ \rho^{(1)}(\mathsf{B}) &= 0.6 \times 0.05 + 0.3 \times 0.2 + 0.1 \times 0.2 = 0.11, \end{aligned}
so \boldsymbol{\rho}^{(1)} = \begin{pmatrix} 0.65 & 0.24 & 0.11 \end{pmatrix}. Feeding this back into the same matrix gives day 2:
\boldsymbol{\rho}^{(2)} = \boldsymbol{\rho}^{(1)} \mathbf{P} = \begin{pmatrix} 0.65 & 0.24 & 0.11 \end{pmatrix} \begin{pmatrix} 0.8 & 0.15 & 0.05 \\ 0.4 & 0.4 & 0.2 \\ 0.5 & 0.3 & 0.2 \end{pmatrix} = \begin{pmatrix} 0.671 & 0.2265 & 0.1025 \end{pmatrix}.
Using the two-step transition matrix. Alternatively, square the transition matrix first,
\mathbf{P}^{(2)} = \begin{pmatrix} 0.725 & 0.195 & 0.08 \\ 0.58 & 0.28 & 0.14 \\ 0.62 & 0.255 & 0.125 \end{pmatrix},
whose entries are the two-step probabilities p^{(2)}_{ij}. For instance, the (\mathsf{G}, \mathsf{B}) entry says that a machine that is good today is broken the day after tomorrow with probability 0.08, collecting the three intermediate possibilities \mathsf{G} \to \mathsf{G} \to \mathsf{B}, \mathsf{G} \to \mathsf{W} \to \mathsf{B}, and \mathsf{G} \to \mathsf{B} \to \mathsf{B}. Applying it to the initial distribution,
\boldsymbol{\rho}^{(2)} = \boldsymbol{\rho}\mathbf{P}^{(2)} = \begin{pmatrix} 0.6 & 0.3 & 0.1 \end{pmatrix} \begin{pmatrix} 0.725 & 0.195 & 0.08 \\ 0.58 & 0.28 & 0.14 \\ 0.62 & 0.255 & 0.125 \end{pmatrix} = \begin{pmatrix} 0.671 & 0.2265 & 0.1025 \end{pmatrix},
the same answer, as it must be since \boldsymbol{\rho}\mathbf{P}^{(2)} = (\boldsymbol{\rho}\mathbf{P}^{(1)})\mathbf{P}. On the third morning of the semester the machine is out of order with probability 0.1025, and pours something drinkable with probability 0.8975.
Repeating the multiplication, the distribution settles down to \begin{pmatrix} 28/41 & 9/41 & 4/41 \end{pmatrix} \approx \begin{pmatrix} 0.683 & 0.220 & 0.098 \end{pmatrix}, and it does so no matter what \boldsymbol{\rho} we started from: in the long run the students face a coffee-less morning about one day in ten.
1.5 Chain Structure
- regular: \mathcal{S} is a single closed and aperiodic class;
- recurrent: \mathcal{S} is a single closed class. It may be periodic or aperiodic;
- unichain: A single recurrent class and some transient states;
- multichain: More than one recurrent class and possibly some transient states.
1.6 Limiting behavior
When \mathbf{P} is regular,
\mathbf{P}^* = \lim_{n \to \infty} \mathbf{P}^n .
When \mathbf{P} is not regular, the limit does not exist (periodic chain). Instead, the limit
\mathbf{P}^* := \lim_{n \to \infty} \frac{1}{n} \sum_{n=0}^{N-1} \mathbf{P}^n
exist.
\mathbf{P}^* \mathbf{P} = \mathbf{P} \mathbf{P}^* = \mathbf{P}^* \mathbf{P}^* = \mathbf{P}^*
When the Markov chain is recurrent,
\mathbf{P}^* = \mathbf{e} \mathbf{q}^\top, where \mathbf{q} is the stationary distribution of the Markov chain. As a result, \mathbf{P}^* has the following structure:
\mathbf{P}^* = \begin{pmatrix} q_0 & q_1 & \cdots & q_m \\ q_0 & q_1 & \cdots & q_m \\ \vdots & \vdots & \ddots & \vdots \\ q_0 & q_1 & \cdots & q_m \end{pmatrix},
and \sum_{j=0}^{m} q_j = 1.
The lemmas below are the matrix analogues of the scalar geometric series \sum_{n=0}^{\infty} q^n = (1 - q)^{-1}, valid for |q| < 1. They all rest on the following elementary identity.
Lemma 1.1 For every \mathbf{Q} \in \mathbb{R}^{|\mathcal{S}| \times |\mathcal{S}|} and every N = 0, 1, 2, \ldots,
(\mathbf{I} - \mathbf{Q}) \sum_{n=0}^{N} \mathbf{Q}^n = \left( \sum_{n=0}^{N} \mathbf{Q}^n \right) (\mathbf{I} - \mathbf{Q}) = \mathbf{I} - \mathbf{Q}^{N+1}.
Proof. Expanding the product and shifting the index of the second sum,
\begin{aligned} (\mathbf{I} - \mathbf{Q}) \sum_{n=0}^{N} \mathbf{Q}^n &= \sum_{n=0}^{N} \mathbf{Q}^n - \sum_{n=0}^{N} \mathbf{Q}^{n+1} \\ &= \sum_{n=0}^{N} \mathbf{Q}^n - \sum_{n=1}^{N+1} \mathbf{Q}^n \\ &= \mathbf{Q}^0 - \mathbf{Q}^{N+1} \\ &= \mathbf{I} - \mathbf{Q}^{N+1}, \end{aligned}
since all the intermediate terms cancel. The same computation with the factors in the other order gives the second equality, because \mathbf{Q} commutes with every power of itself.
Since \mathcal{S} is finite, all norms on \mathbb{R}^{|\mathcal{S}| \times |\mathcal{S}|} are equivalent, so the statement \mathbf{Q}^N \rightarrow \mathbf{0} is unambiguous: it holds entrywise if and only if it holds in any matrix norm. Below, \|\cdot\| denotes any submultiplicative matrix norm, so that \|\mathbf{A} \mathbf{B}\| \leq \|\mathbf{A}\| \|\mathbf{B}\|.
Lemma 1.2 Let \mathbf{Q} \in \mathbb{R}^{|\mathcal{S}| \times |\mathcal{S}|}. Then \mathbf{Q}^N \rightarrow \mathbf{0} as N \rightarrow \infty if and only if \sigma(\mathbf{Q}) < 1, where
\sigma(\mathbf{Q}) := \max \left\{ |\lambda| : \lambda \text{ is an eigenvalue of } \mathbf{Q} \right\}
is the spectral radius of \mathbf{Q}.
Proof. Suppose first that \sigma(\mathbf{Q}) < 1. Write \mathbf{Q} = \mathbf{T} \mathbf{J} \mathbf{T}^{-1} in Jordan form, so that \mathbf{Q}^N = \mathbf{T} \mathbf{J}^N \mathbf{T}^{-1} and it suffices to show that \mathbf{J}^N \rightarrow \mathbf{0} block by block. A Jordan block of size m associated with the eigenvalue \lambda is \mathbf{J}_m(\lambda) = \lambda \mathbf{I}_m + \mathbf{N}, where \mathbf{N} is nilpotent with \mathbf{N}^m = \mathbf{0}. The two summands commute, so the binomial theorem gives
\mathbf{J}_m(\lambda)^N = \sum_{k=0}^{\min(N,\, m-1)} \binom{N}{k} \lambda^{N-k} \mathbf{N}^k ,
where only the finitely many indices k \leq m - 1 survive because \mathbf{N}^k = \mathbf{0} beyond them. Each coefficient obeys
\left| \binom{N}{k} \lambda^{N-k} \right| \leq \frac{N^k}{k!} \, |\lambda|^{N-k} \xrightarrow{\;N \to \infty\;} 0 ,
since |\lambda| \leq \sigma(\mathbf{Q}) < 1 and geometric decay dominates polynomial growth. Hence every entry of \mathbf{J}_m(\lambda)^N tends to 0, and therefore \mathbf{Q}^N \rightarrow \mathbf{0}.
Conversely, suppose \sigma(\mathbf{Q}) \geq 1, and let \lambda be an eigenvalue with |\lambda| \geq 1 and eigenvector \mathbf{v} \neq \mathbf{0}, both possibly complex. Then \mathbf{Q}^N \mathbf{v} = \lambda^N \mathbf{v}, so
\left\| \mathbf{Q}^N \mathbf{v} \right\| = |\lambda|^N \|\mathbf{v}\| \geq \|\mathbf{v}\| > 0 \qquad \text{for all } N,
which rules out \mathbf{Q}^N \rightarrow \mathbf{0}. Contraposition gives the claim.
Lemma 1.3 Let \mathbf{Q} \in \mathbb{R}^{|\mathcal{S}| \times |\mathcal{S}|} be a matrix which \mathbf{Q}^N \rightarrow \mathbf{0} as N \rightarrow \infty. Then, the (\mathbf{I} - \mathbf{Q})^{-1} exists and
\sum_{n=0}^{\infty} \mathbf{Q}^n = (\mathbf{I} - \mathbf{Q})^{-1}.
Proof. We first show that \mathbf{I} - \mathbf{Q} is nonsingular. Let \mathbf{x} satisfy (\mathbf{I} - \mathbf{Q}) \mathbf{x} = \mathbf{0}, that is, \mathbf{Q} \mathbf{x} = \mathbf{x}. Iterating gives \mathbf{Q}^N \mathbf{x} = \mathbf{x} for every N, and letting N \rightarrow \infty yields \mathbf{x} = \mathbf{0}, because \mathbf{Q}^N \rightarrow \mathbf{0}. So the null space of \mathbf{I} - \mathbf{Q} is trivial and, the matrix being square, (\mathbf{I} - \mathbf{Q})^{-1} exists.
Write \mathbf{S}_N := \sum_{n=0}^{N} \mathbf{Q}^n for the partial sums. Multiplying the identity of Lemma 1.1 on the left by (\mathbf{I} - \mathbf{Q})^{-1},
\mathbf{S}_N = (\mathbf{I} - \mathbf{Q})^{-1} \left( \mathbf{I} - \mathbf{Q}^{N+1} \right) = (\mathbf{I} - \mathbf{Q})^{-1} - (\mathbf{I} - \mathbf{Q})^{-1} \mathbf{Q}^{N+1} .
The matrix (\mathbf{I} - \mathbf{Q})^{-1} is fixed, so the remainder satisfies
\left\| (\mathbf{I} - \mathbf{Q})^{-1} \mathbf{Q}^{N+1} \right\| \leq \left\| (\mathbf{I} - \mathbf{Q})^{-1} \right\| \left\| \mathbf{Q}^{N+1} \right\| \xrightarrow{\;N \to \infty\;} 0 .
Therefore \mathbf{S}_N converges, and \sum_{n=0}^{\infty} \mathbf{Q}^n = \lim_{N \rightarrow \infty} \mathbf{S}_N = (\mathbf{I} - \mathbf{Q})^{-1}.
Lemma 1.4 Let \mathbf{Q} \in \mathbb{R}^{|\mathcal{S}| \times |\mathcal{S}|}. The series \sum_{n=0}^{\infty} \mathbf{Q}^n converges if and only if \sigma(\mathbf{Q}) < 1, where \sigma(\mathbf{Q}) is the spectral radius of \mathbf{Q}. In that case, the (\mathbf{I} - \mathbf{Q})^{-1} exists and
\sum_{n=0}^{\infty} \mathbf{Q}^n = (\mathbf{I} - \mathbf{Q})^{-1}.
Proof. Suppose \sigma(\mathbf{Q}) < 1. Then \mathbf{Q}^N \rightarrow \mathbf{0} by Lemma 1.2, and Lemma 1.3 gives both the existence of (\mathbf{I} - \mathbf{Q})^{-1} and the value of the sum.
Conversely, suppose the partial sums \mathbf{S}_N = \sum_{n=0}^{N} \mathbf{Q}^n converge to some limit \mathbf{S}. Then
\mathbf{Q}^N = \mathbf{S}_N - \mathbf{S}_{N-1} \xrightarrow{\;N \to \infty\;} \mathbf{S} - \mathbf{S} = \mathbf{0},
which is the familiar statement that the terms of a convergent series vanish. Applying Lemma 1.2 once more gives \sigma(\mathbf{Q}) < 1.
Remark. The biconditional in Lemma 1.4 is a statement about the series, not about the inverse alone. Invertibility of \mathbf{I} - \mathbf{Q} places no bound whatsoever on the spectral radius: for \mathbf{Q} = 2 \mathbf{I} we have (\mathbf{I} - \mathbf{Q})^{-1} = -\mathbf{I}, yet \sigma(\mathbf{Q}) = 2 and \sum_{n=0}^{\infty} \mathbf{Q}^n diverges. Read together, Lemma 1.2 and Lemma 1.4 say that the three conditions
\sigma(\mathbf{Q}) < 1, \qquad \mathbf{Q}^N \rightarrow \mathbf{0}, \qquad \sum_{n=0}^{\infty} \mathbf{Q}^n \text{ converges}
are equivalent, and that Lemma 1.3 is the implication from the second to the third.
Lemma 1.5 Let \mathbf{Q} \in \mathbb{R}^{|\mathcal{S}| \times |\mathcal{S}|} be a matrix which \lim_{N \rightarrow \infty} \frac{1}{N} \sum_{n=0}^{N-1} \mathbf{Q}^n = \mathbf{0}. Then, the (\mathbf{I} - \mathbf{Q})^{-1} exists and
\lim_{N \rightarrow \infty} \frac{1}{N} \sum_{k=1}^{N} \sum_{n=0}^{k} \mathbf{Q}^n = (\mathbf{I} - \mathbf{Q})^{-1}.
Proof. Write \mathbf{A}_N := \frac{1}{N} \sum_{n=0}^{N-1} \mathbf{Q}^n for the Cesàro averages, so that \mathbf{A}_N \rightarrow \mathbf{0} by hypothesis, and \mathbf{S}_k := \sum_{n=0}^{k} \mathbf{Q}^n for the partial sums.
We first show that \mathbf{I} - \mathbf{Q} is nonsingular. Let \mathbf{Q} \mathbf{x} = \mathbf{x}. Then \mathbf{Q}^n \mathbf{x} = \mathbf{x} for every n, so that
\mathbf{A}_N \mathbf{x} = \frac{1}{N} \sum_{n=0}^{N-1} \mathbf{Q}^n \mathbf{x} = \frac{1}{N} \cdot N \mathbf{x} = \mathbf{x} \qquad \text{for all } N,
and letting N \rightarrow \infty gives \mathbf{x} = \mathbf{0}. Hence 1 is not an eigenvalue of \mathbf{Q} and (\mathbf{I} - \mathbf{Q})^{-1} exists. Note that the hypothesis here does not give \mathbf{Q}^N \rightarrow \mathbf{0}, so Lemma 1.3 cannot be invoked directly.
By Lemma 1.1 we have \mathbf{S}_k = (\mathbf{I} - \mathbf{Q})^{-1} \left( \mathbf{I} - \mathbf{Q}^{k+1} \right), and averaging over k = 1, \ldots, N,
\frac{1}{N} \sum_{k=1}^{N} \mathbf{S}_k = (\mathbf{I} - \mathbf{Q})^{-1} \left[ \mathbf{I} - \frac{1}{N} \sum_{k=1}^{N} \mathbf{Q}^{k+1} \right].
It remains to show that the average inside the bracket vanishes. Reindexing with m = k + 1 and completing the sum,
\sum_{k=1}^{N} \mathbf{Q}^{k+1} = \sum_{m=2}^{N+1} \mathbf{Q}^{m} = \sum_{m=0}^{N+1} \mathbf{Q}^{m} - \mathbf{I} - \mathbf{Q} = (N+2) \mathbf{A}_{N+2} - \mathbf{I} - \mathbf{Q},
so that
\frac{1}{N} \sum_{k=1}^{N} \mathbf{Q}^{k+1} = \frac{N+2}{N} \mathbf{A}_{N+2} - \frac{1}{N} \left( \mathbf{I} + \mathbf{Q} \right) \xrightarrow{\;N \to \infty\;} 1 \cdot \mathbf{0} - \mathbf{0} = \mathbf{0},
using \frac{N+2}{N} \rightarrow 1 and \mathbf{A}_{N+2} \rightarrow \mathbf{0}. Multiplication by the fixed matrix (\mathbf{I} - \mathbf{Q})^{-1} preserves limits, and therefore
\lim_{N \rightarrow \infty} \frac{1}{N} \sum_{k=1}^{N} \sum_{n=0}^{k} \mathbf{Q}^n = (\mathbf{I} - \mathbf{Q})^{-1} \left[ \mathbf{I} - \mathbf{0} \right] = (\mathbf{I} - \mathbf{Q})^{-1}.
Remark. The hypothesis of Lemma 1.5 is strictly weaker than that of Lemma 1.3, which is exactly why it is the right tool for periodic chains. Take |\mathcal{S}| = 1 and \mathbf{Q} = (-1). Its powers (-1)^N oscillate and never converge, so \sum_{n=0}^{\infty} \mathbf{Q}^n diverges and Lemma 1.3 does not apply, yet the averages \frac{1}{N} \sum_{n=0}^{N-1} (-1)^n do tend to 0. The partial sums \mathbf{S}_k alternate between 1 and 0, and their averages converge to \frac{1}{2} = \left( 1 - (-1) \right)^{-1}, in agreement with the lemma.
Theorem 1.1 Let d \in D^{\text{MR}}. Then
\mathbf{h}^{d^{\infty}} = \mathbf{H}_d \mathbf{r}^d,
where \mathbf{H}_d := (\mathbf{I} - (\mathbf{P}_d - \mathbf{P}^*_d))^{-1}(\mathbf{I} - \mathbf{P}^*_d).
Proof. Given \mathbf{h}^{d^{\infty}} = \sum_{n=1}^{\infty} (\mathbf{P}_d^{n-1} - \mathbf{P}_d^*) \mathbf{r}^d, we need to show that \mathbf{H}_d = \sum_{n=1}^{\infty} (\mathbf{P}_d^{n-1} - \mathbf{P}_d^*).
Throughout we use the identities \mathbf{P}_d^* \mathbf{P}_d = \mathbf{P}_d \mathbf{P}_d^* = \mathbf{P}_d^* \mathbf{P}_d^* = \mathbf{P}_d^* recorded earlier in this section. Applying \mathbf{P}_d \mathbf{P}_d^* = \mathbf{P}_d^* repeatedly gives
\mathbf{P}_d^{n} \mathbf{P}_d^* = \mathbf{P}_d^{n-1} \left( \mathbf{P}_d \mathbf{P}_d^* \right) = \mathbf{P}_d^{n-1} \mathbf{P}_d^* = \cdots = \mathbf{P}_d^* , \qquad n \geq 1,
and symmetrically \mathbf{P}_d^* \mathbf{P}_d^{n} = \mathbf{P}_d^*.
When n \geq 1, we have \mathbf{P}_d^{n} - \mathbf{P}_d^* = (\mathbf{P}_d - \mathbf{P}_d^*)^{n}. This is trivial for n = 1, and if it holds for some n \geq 1, then
\begin{aligned} \left( \mathbf{P}_d - \mathbf{P}_d^* \right)^{n+1} &= \left( \mathbf{P}_d^{n} - \mathbf{P}_d^* \right) \left( \mathbf{P}_d - \mathbf{P}_d^* \right) \\ &= \mathbf{P}_d^{n+1} - \mathbf{P}_d^{n} \mathbf{P}_d^* - \mathbf{P}_d^* \mathbf{P}_d + \mathbf{P}_d^* \mathbf{P}_d^* \\ &= \mathbf{P}_d^{n+1} - \mathbf{P}_d^* - \mathbf{P}_d^* + \mathbf{P}_d^* \\ &= \mathbf{P}_d^{n+1} - \mathbf{P}_d^* , \end{aligned}
so the identity holds for all n \geq 1 by induction, so that
\begin{aligned} \sum_{n=1}^{\infty} \left( \mathbf{P}_d^{n-1} - \mathbf{P}_d^* \right) &= \sum_{n=0}^{\infty} \left( \mathbf{P}_d^{n} - \mathbf{P}_d^* \right) \\ &= \left( \mathbf{I} - \mathbf{P}_d^* \right) + \sum_{n=1}^{\infty} \left( \mathbf{P}_d - \mathbf{P}_d^* \right)^{n} . \end{aligned}
The term n = 0 has to be separated here, because \mathbf{P}_d^0 - \mathbf{P}_d^* = \mathbf{I} - \mathbf{P}_d^* while \left( \mathbf{P}_d - \mathbf{P}_d^* \right)^0 = \mathbf{I}; this is the one index at which the identity just proved fails.
Write \mathbf{Q} := \mathbf{P}_d - \mathbf{P}_d^*, so that \mathbf{H}_d = (\mathbf{I} - \mathbf{Q})^{-1} (\mathbf{I} - \mathbf{P}_d^*). Since \mathbf{Q}^0 = \mathbf{I}, the display above reads
\sum_{n=1}^{\infty} \left( \mathbf{P}_d^{n-1} - \mathbf{P}_d^* \right) = \sum_{n=0}^{\infty} \mathbf{Q}^n - \mathbf{P}_d^* .
The inverse exists. By the definition of \mathbf{P}_d^* as a Cesàro limit,
\begin{aligned} \frac{1}{N} \sum_{n=0}^{N-1} \mathbf{Q}^n &= \frac{1}{N} \left[ \mathbf{I} + \sum_{n=1}^{N-1} \left( \mathbf{P}_d^{n} - \mathbf{P}_d^* \right) \right] \\ &= \frac{1}{N} \sum_{n=0}^{N-1} \mathbf{P}_d^{n} - \frac{N-1}{N} \mathbf{P}_d^* \;\xrightarrow{\;N \to \infty\;}\; \mathbf{P}_d^* - \mathbf{P}_d^* = \mathbf{0} . \end{aligned}
So \mathbf{Q} satisfies the hypothesis of Lemma 1.5, and (\mathbf{I} - \mathbf{Q})^{-1} exists. In particular \mathbf{H}_d is well defined for every d \in D^{\text{MR}}, periodic chains included.
The series equals the inverse. Writing the premise as a matrix series presupposes that \sum_{n=0}^{\infty} \mathbf{Q}^n converges. By Lemma 1.4 that is equivalent to \mathbf{Q}^N \rightarrow \mathbf{0}, that is, to \mathbf{P}_d^N \rightarrow \mathbf{P}_d^*, which holds precisely when the chain induced by \mathbf{P}_d is aperiodic. Under this assumption Lemma 1.3 gives
\sum_{n=0}^{\infty} \mathbf{Q}^n = (\mathbf{I} - \mathbf{Q})^{-1} .
Absorbing the \mathbf{P}_d^* term. Finally,
\mathbf{Q} \mathbf{P}_d^* = \left( \mathbf{P}_d - \mathbf{P}_d^* \right) \mathbf{P}_d^* = \mathbf{P}_d \mathbf{P}_d^* - \mathbf{P}_d^* \mathbf{P}_d^* = \mathbf{P}_d^* - \mathbf{P}_d^* = \mathbf{0},
so (\mathbf{I} - \mathbf{Q}) \mathbf{P}_d^* = \mathbf{P}_d^*, and multiplying on the left by (\mathbf{I} - \mathbf{Q})^{-1} gives \mathbf{P}_d^* = (\mathbf{I} - \mathbf{Q})^{-1} \mathbf{P}_d^*. Combining the last three displays,
\begin{aligned} \sum_{n=1}^{\infty} \left( \mathbf{P}_d^{n-1} - \mathbf{P}_d^* \right) &= \sum_{n=0}^{\infty} \mathbf{Q}^n - \mathbf{P}_d^* \\ &= (\mathbf{I} - \mathbf{Q})^{-1} - (\mathbf{I} - \mathbf{Q})^{-1} \mathbf{P}_d^* \\ &= (\mathbf{I} - \mathbf{Q})^{-1} \left( \mathbf{I} - \mathbf{P}_d^* \right) \\ &= \left( \mathbf{I} - \left( \mathbf{P}_d - \mathbf{P}_d^* \right) \right)^{-1} \left( \mathbf{I} - \mathbf{P}_d^* \right) = \mathbf{H}_d , \end{aligned}
which is the required identity. Multiplying on the right by \mathbf{r}^d gives \mathbf{h}^{d^{\infty}} = \mathbf{H}_d \mathbf{r}^d.