Hidden Markov Models

Take two regimes, calm and turbulent. On a calm day the return (in percent) is drawn from $\mathcal{N}(0.05, 0.7^2)$, a Normal distribution with mean 0.05 and standard deviation 0.7; on a turbulent day from $\mathcal{N}(-0.10, 2.0^2)$. A calm day is followed by another calm day with probability 0.98, and a turbulent day by another turbulent day with probability 0.95. Simulating $T = 1000$ days:

• The series is turbulent on a fraction 0.362 of days. The long-run value is $0.02/(0.02 + 0.05) = 0.2857$ (Example 3 derives it), but the regimes are so persistent that 1,000 days contain only 37 switches, and the sample fraction wanders.
• The returns have sample excess kurtosis 2.77: fat tails, from a mixture of two thin-tailed Normals.
• The lag-1 autocorrelation of the returns is −0.018, essentially zero. The lag-1 autocorrelation of the squared returns is 0.206.

That last pair is the fingerprint. Returns look unpredictable, and their magnitudes are strongly predictable. Example 5 derives this in closed form.

A casino rolls one of two dice: a fair die and a loaded die that shows a six with probability 0.5 and each other face with probability 0.1. Between rolls it switches from the fair die to the loaded one with probability 0.05, and back with probability 0.10. We see only the 300 rolls. Which stretches were rolled with the loaded die?

We return to this example throughout §34.4–§34.6. Its emissions are dice, so every posterior probability in it has a concrete reading: “how likely is it that this roll came from the loaded die?”

We write $Z_t \in \{1, \dots, K\}$ for the hidden state at time $t$ and $X_t$ for the observation, with lowercase $x_t$ for an observed value. This matches the EM appendix of Topic 14, where $Z$ was the latent label and $X$ the data. formalML’s sequential Monte Carlo page writes $x_t$ for the hidden state and $y_t$ for the observation; when crossing between the sites, swap the letters. We write $x_{s:t} = (x_s, \dots, x_t)$ for a block of consecutive observations.

Every use of an HMM reduces to one of three computations:

1. Evaluate: compute the likelihood $p(x_{1:T})$ of an observed sequence (§34.4).
2. Decode: recover the hidden states, either one time at a time (§34.5) or as a single best path (§34.6).
3. Learn: estimate the parameters from data (§34.7).

Each has a brute-force solution that costs $K^T$ operations, exponential in the sequence length. The point of this topic is that the HMM’s structure brings each one down to $O(TK^2)$. The picture to hold on to is the trellis: a $K \times T$ grid with one node per (state, time) pair and edges between consecutive columns. Every algorithm here is one dynamic-programming sweep across that grid, and the only thing that changes between them is whether we sum over incoming edges or maximize over them.

A sequence $Z_1, Z_2, \dots$ taking values in $\{1, \dots, K\}$ is a (time-homogeneous) Markov chain if the next state depends on the past only through the present:

The transition matrix $A = (a_{jk})$ has nonnegative entries and rows that sum to 1 (it is row-stochastic): row $j$ is the distribution of the next state given the current state $j$. The initial distribution is the row vector $\nu$ with $\nu_k = P(Z_1 = k)$. This is Topic 26’s transition kernel (§26.2 Def 1) on a finite state space.

For a finite Markov chain with transition matrix $A$ and initial distribution $\nu$:

(i) $P(Z_{t+n} = k \mid Z_t = j) = (A^n)_{jk}$ for all $t, n \ge 1$;

(ii) $A^{m+n} = A^m A^n$, that is, $(A^{m+n})_{jk} = \sum_l (A^m)_{jl}(A^n)_{lk}$;

(iii) the marginal law of $Z_t$ is the row vector $\nu A^{t-1}$.

A distribution $\pi$ on $\{1, \dots, K\}$ is stationary for $A$ if $\pi A = \pi$: starting the chain from $\pi$, every $Z_t$ has law $\pi$. Geometrically, $\pi$ is a left eigenvector of $A$ with eigenvalue 1 ( formalCalculus: Eigenvalues & Eigenvectors ). The chain is irreducible if every state can be reached from every other, and aperiodic if for each state $k$ the return times $\{n : (A^n)_{kk} > 0\}$ have greatest common divisor 1.

Let $A$ be the transition matrix of a finite, irreducible, aperiodic chain. Then (a) there is a unique stationary distribution $\pi$, and $\pi_k > 0$ for every $k$; (b) there are constants $C < \infty$ and $\rho \in (0, 1)$ such that $\max_k |(\nu A^n)_k - \pi_k| \le C\rho^n$ for every initial distribution $\nu$; (c) $\pi_k = 1/E[\tau_k \mid Z_1 = k]$, where $\tau_k$ is the time the chain takes to return to $k$.

Write the two-state transition matrix as

so $p$ is the probability of leaving state 1 and $q$ the probability of leaving state 2. Multiplying out shows that $\pi = (q, p)/(p + q)$ satisfies $\pi A = \pi$. The eigenvalues of $A$ are 1 and $\lambda = 1 - p - q$.

To see how fast the chain forgets its start, let $\mu_n = \nu A^n$ and track its first coordinate. By Theorem 1,

The stationary value satisfies the same equation, $\pi_1 = q + \lambda\pi_1$. Subtracting,

and induction gives $(\nu A^n)_1 - \pi_1 = (\nu_1 - \pi_1)\lambda^n$. The distance to stationarity shrinks by the factor $\lambda$ every step. Starting from a known state $j$ (so $\nu$ puts mass 1 on $j$), the same identity in each coordinate gives the $h$-step transition probabilities we will need in Example 5:

where $\mathbf{1}\{j = k\}$ is 1 if $j = k$ and 0 otherwise.

For the market chain, $p = 0.02$ and $q = 0.05$, so $\pi = (0.7143, 0.2857)$ and $\lambda = 0.93$. The influence of the starting regime halves every $\log 2/(-\log 0.93) \approx 9.55$ days.

Suppose the chain enters state $k$ at time $s$ (that is, $Z_s = k$, and either $s = 1$ or $Z_{s-1} \ne k$), and let $D = \min\{d \ge 1 : Z_{s+d} \ne k\}$ be the number of consecutive steps it spends in $k$. If $a_{kk} < 1$, then

so $D \sim \mathrm{Geometric}(1 - a_{kk})$ and $E[D] = 1/(1 - a_{kk})$.

With $K > 2$ states, the second-largest eigenvalue modulus $|\lambda_2|$ of $A$ plays the role of $\lambda$: the chain forgets its starting point at roughly the rate $|\lambda_2|^n$, and $1 - |\lambda_2|$ is called the spectral gap. That is the finite-state face of a mixing-time bound, developed in formalML: random-walks . We won’t need more than the two-state case.

Requiring the next state to depend only on the current one costs less generality than it seems. A chain where the next state depends on the last $m$ states is a first-order chain on the $K^m$ possible blocks of $m$ states, so “first order” restricts only how big the state space is. The genuine restriction is Theorem 3: every state’s dwell time is geometric, so the most likely duration is always 1, whatever the mean. Regimes that tend to last about a month, or phonemes with a typical length, violate that. §34.9 Remark 21 gives the fix.

A hidden Markov model with $K$ states has parameters $\theta = (\nu, A, \{b_k\})$ and satisfies:

(i) $Z_{1:T}$ is a Markov chain with initial distribution $\nu$ and transition matrix $A$;

(ii) given the whole state sequence $Z_{1:T}$, the observations $X_1, \dots, X_T$ are independent, and $X_t$ has distribution $b_{Z_t}$, which depends only on the current state.

Here $b_k$ is a pmf or density, the emission distribution of state $k$. In words: the state walks; at every step it emits one observation from its own distribution and otherwise leaves no trace. Categorical emissions have $b_k(v) = B_{kv}$ for symbols $v$ in a finite alphabet (the casino). Gaussian emissions have $b_k = \mathcal{N}(\mu_k, \sigma_k^2)$ (the market).

Under Definition 3, for every $T$ and every $(z_{1:T}, x_{1:T})$:

(i) the joint probability factors along the chain,

(ii) for each $t$, the past block $(Z_{1:t-1}, X_{1:t})$ and the future block $(Z_{t+1:T}, X_{t+1:T})$ are conditionally independent given $Z_t$; in particular $p(x_{t+1:T} \mid z_t, z_{1:t-1}, x_{1:t}) = p(x_{t+1:T} \mid z_t)$;

(iii) $p(x_t \mid z_t, z_{1:t-1}, x_{1:t-1}) = b_{z_t}(x_t)$ and $p(z_{t+1} \mid z_{1:t}, x_{1:t}) = a_{z_t z_{t+1}}$.

In the casino, start the chain at its stationary distribution $\nu = \pi = (2/3, 1/3)$ (fair, loaded). Take the rolls $x_{1:3} = (6, 6, 1)$ and the path $z_{1:3} = (\mathrm{L}, \mathrm{L}, \mathrm{F})$. Theorem 4(i) multiplies one factor per edge and node of the path:

Summing this product over all $2^3 = 8$ paths gives the likelihood of the rolls, $p(x_{1:3}) = 0.0116235$. That 8-term sum is the brute-force computation §34.4 replaces.

Take a two-state HMM started at stationarity with Gaussian emissions, and write $\mu_k$ and $\sigma_k^2$ for the mean and variance of state $k$. Fix a lag $h \ge 1$. Given the whole state path, $X_t$ and $X_{t+h}$ are independent (Theorem 4), so their covariance comes entirely from the states. Conditioning on the path and using the tower property,

Now use the $h$-step formula from Example 3:

Subtracting $E[\mu_{Z_t}]E[\mu_{Z_{t+h}}] = (\pi\cdot\mu)^2$ leaves $\lambda^h$ times the variance of $\mu_Z$. For two states that variance is $\pi_1\pi_2(\mu_1 - \mu_2)^2$, so

The same argument applied to $X_t^2$, with $s_k = E[X_t^2 \mid Z_t = k] = \mu_k^2 + \sigma_k^2$ in place of $\mu_k$, gives

In the market model the means barely differ (0.05 vs −0.10), so returns are nearly uncorrelated: the theoretical lag-1 autocorrelation is 0.0029. The second moments differ a lot ($s_1 = 0.49$, $s_2 = 4.01$), so squared returns are strongly autocorrelated, decaying like $0.93^h$. At lags $h = 1, 5, 20$ the theory gives 0.1947, 0.1457, 0.0490; a 200,000-day simulation gives 0.1893, 0.1450, 0.0471, and the 1,000-day series of Example 1 gives 0.206, 0.124, 0.043. That is volatility clustering, produced by nothing more than a persistent hidden state. The model has no feedback from past returns to future variance; GARCH builds that feedback in directly, and Topic 35 fits it to this same data.

If every row of $A$ equals $\pi$, the next state ignores the current one, the labels $Z_t$ are i.i.d. draws from $\pi$, and the HMM is exactly the Gaussian mixture of Topic 14’s EM appendix. Then $\lambda = 0$ and every autocovariance in Example 5 vanishes. The HMM’s extra content over a mixture is persistence: the cluster label at time $t$ predicts the label at time $t + 1$. The “No memory” preset in the simulator below makes this concrete. Mixtures and their EM are developed in formalML: clustering .

Figure 3’s directed graph encodes Theorem 4 by itself. Conditional independence (ii) is d-separation: every path from the past to the future passes through $Z_t$, and observing $Z_t$ blocks it. Naive Bayes reads the same picture one time step at a time: one hidden class, many conditionally independent features. The HMM strings naive-Bayes models together along a Markov chain.

The forward variable is the joint probability of the observations so far and the current state,

The filtering distribution is $P(Z_t = k \mid x_{1:t})$: what we believe about the state now, having seen everything up to now and nothing after.

The forward variables satisfy $\alpha_1(k) = \nu_k b_k(x_1)$ and, for $t \ge 2$,

Consequently $p(x_{1:T}) = \sum_k \alpha_T(k)$, and computing it takes $O(TK^2)$ operations.

The forward algorithm costs about $2K^2(T - 1) + 2K$ multiply-adds: 2,396 for the 300 casino rolls. Brute force costs about $TK^T \approx$ 6.1 × 10⁹². Both compute the same number. For the first 10 casino rolls (faces 1, 2, 5, 6, 4, 4, 2, 2, 1, 4), forward and brute-force enumeration of all $2^{10}$ paths both give $\log p(x_{1:10}) =$ −18.488737947800, identical to double precision.

Define $\hat\alpha_1(k) = \nu_k b_k(x_1)/c_1$ with $c_1 = \sum_k \nu_k b_k(x_1)$, and for $t \ge 2$

Then:

(i) $c_t = p(x_t \mid x_{1:t-1})$, with $c_1 = p(x_1)$;

(ii) $\hat\alpha_t(k) = \alpha_t(k)/p(x_{1:t}) = P(Z_t = k \mid x_{1:t})$;

(iii) $\log p(x_{1:T}) = \sum_{t=1}^T \log c_t$.

Every normalized quantity lies in $[0, 1]$, so nothing underflows.

Run the unscaled recursion of Theorem 5 on the 1,000-day market series in double precision. The raw $\alpha_t$ becomes exactly 0 at t = 542; after that the “likelihood” is $\log 0$. The scaled recursion runs to the end without trouble and gives $\log p(x_{1:1000}) =$ −1522.1208729, which matches an independent log-sum-exp implementation to 1.4 × 10⁻¹² (about 15 significant digits). Figure 5’s right panel is the filter itself: $\hat\alpha_t(\text{turbulent})$, jumping toward 1 at the first large return of each turbulent stretch.

Scaling (Theorem 6) is one fix; the other is to run the recursion on log probabilities, replacing each sum with log-sum-exp, $\log\sum_j e^{u_j} = m + \log\sum_j e^{u_j - m}$ with $m = \max_j u_j$. Both are exact. The site’s implementation uses scaling, because it produces the filtering distribution for free and keeps the backward pass of §34.5 in the same units.

Predict, then update: this is the same two-step recursion formalML’s sequential Monte Carlo page writes for general state spaces, with integrals over a continuous state in place of our sums over $K$ states. When the state is continuous and the model non-linear, those integrals have no closed form, and the bootstrap particle filter replaces the exact sum with a weighted sample of states ( formalML: sequential-monte-carlo ). The finite-state forward pass is the exact case it approximates. The same recursion appears in probabilistic programming: Stan and PyMC can’t sample discrete hidden states with gradient-based samplers, so they marginalize them out with exactly this forward pass ( formalML: probabilistic-programming ).

The backward variable is the probability of the future observations given the current state,

The smoothing marginals are $\gamma_t(k) = P(Z_t = k \mid x_{1:T})$, the posterior probability of state $k$ at time $t$ given all the data, and $\xi_t(j,k) = P(Z_t = j, Z_{t+1} = k \mid x_{1:T})$, the posterior probability that the transition from $t$ to $t + 1$ was $j \to k$.

The backward variables satisfy $\beta_T(j) = 1$ and, for $t < T$,

For every $t$:

(i) $\alpha_t(k)\beta_t(k) = p(x_{1:T}, Z_t = k)$, hence $\sum_k \alpha_t(k)\beta_t(k) = p(x_{1:T})$ for every $t$;

(ii) $\gamma_t(k) = \alpha_t(k)\beta_t(k)/p(x_{1:T})$;

(iii) $\xi_t(j,k) = \alpha_t(j)\,a_{jk}\,b_k(x_{t+1})\,\beta_{t+1}(k)/p(x_{1:T})$ for $t < T$.

With the scaled forward variables of Theorem 6 and $\hat\beta_t = \beta_t/\prod_{s > t} c_s$, the same formulas hold with $\hat\alpha, \hat\beta$ in place of $\alpha, \beta$: the denominator $p(x_{1:T})$ disappears from (ii), so $\gamma_t = \hat\alpha_t \odot \hat\beta_t$ (elementwise product), and $1/c_{t+1}$ replaces $1/p(x_{1:T})$ in (iii).

On the 300 casino rolls (loaded 40.7% of the time), threshold each posterior at 0.5 and compare with the truth. The filter labels 75.7% of rolls correctly and the smoother 85.0%. The mean log-loss of the probabilities against the true states falls from 0.478 to 0.354. At the loaded run that starts at roll 85, the smoother crosses 0.5 at roll 88 and the filter only at roll 91: the smoother has already seen the run of sixes that follows. The filter also raises more false alarms. A cluster of sixes near roll 38 pushes the filtered probability to 0.75, and the smoother, seeing the fair rolls that come after, holds it to 0.58.

On the 1,000-day market series, the true regime path switches 37 times. Thresholding the filtered probability of turbulence at 0.5 switches 53 times: the filter flickers whenever a single large calm-day return looks turbulent. The smoother switches 29 times. Its accuracy is 94.8% against the filter’s 89.6%.

The ML reading is a warning. A backtest that conditions trading decisions on smoothed regimes uses information from the future: at day $t$ the smoother has seen days $t + 1, \dots, T$. A live system gets only the filter. Using $\gamma_t$ where $\hat\alpha_t$ belongs is a textbook leakage bug, and the 5-point accuracy gap above is roughly how much it can flatter a backtest.

The filter is online and the smoother is offline. In between is the fixed-lag smoother $P(Z_t \mid x_{1:t+L})$: wait $L$ steps, then commit to a belief about time $t$. It costs a delay of $L$ and recovers most of the smoother’s accuracy when $L$ is a few multiples of the typical regime length. Which one to use is a latency decision, not a statistical one.

Forward–backward is belief propagation on the HMM’s chain-shaped graph. The forward variables are messages passed left to right and the backward variables are messages passed right to left. On any tree-shaped graph the same messages give exact marginals. On graphs with loops, running them anyway gives “loopy” belief propagation, an approximation that works surprisingly well in practice. The chain is the cleanest case: every message is exact and the whole computation is two sweeps.

The MAP path is the single most probable hidden path,

Posterior (pointwise) decoding picks the most probable state at each time separately: $\hat z_t = \arg\max_k \gamma_t(k)$. The first is “the single most probable story”; the second is “the most probable state at each moment, considered one moment at a time.”

Let $\delta_t(k) = \max_{z_{1:t-1}} p(z_{1:t-1}, Z_t = k, x_{1:t})$, the probability of the best path that ends in state $k$ at time $t$, with $\delta_1(k) = \nu_k b_k(x_1)$. Then for $t \ge 2$

Set $z^*_T \in \arg\max_k \delta_T(k)$ and backtrack with $z^*_t = \psi_{t+1}(z^*_{t+1})$ for $t = T-1, \dots, 1$. Then $z^*_{1:T}$ maximizes $p(z_{1:T} \mid x_{1:T})$ over all $K^T$ paths, and $\max_z p(z_{1:T}, x_{1:T}) = \max_k \delta_T(k)$. The cost is $O(TK^2)$ time and $O(TK)$ memory for the back-pointers $\psi$.