Hidden Markov Models

"You think you know which mood I am in, but you only ever see what I emit, never what I am. By the time your forward pass is sure I was calm, I have already, quietly and with full probability, switched to panic."

A Hidden Markov Model Switching Regimes When You Least Expect

In Section 7.2 the Kalman filter tracked a continuous hidden state $\mathbf{z}_t \in \mathbb{R}^d$ through a linear-Gaussian world, alternating a prediction step (push the belief forward through the dynamics) with an update step (correct the belief using the new measurement). Everything in that section assumed the latent variable lived on a continuum and that all distributions stayed Gaussian, which is what made the integrals collapse into the tidy matrix algebra of the Kalman gain. The Hidden Markov Model keeps the recursive skeleton and throws away the continuum. The hidden variable $z_t$ now takes one of $K$ discrete values, the prediction step becomes a matrix-vector product against a transition table, and the update step becomes an elementwise multiply by emission likelihoods followed by renormalization. The same belief, propagated the same way, over a different kind of state. This section makes that correspondence precise and then exploits it.

1. The Discrete-State Cousin of the Kalman Filter Beginner

An HMM is defined by three pieces. There is a hidden state sequence $z_1, z_2, \dots, z_T$ where each $z_t \in \{1, \dots, K\}$, and an observed sequence $x_1, \dots, x_T$. The hidden states form a first-order Markov chain: the next state depends only on the current one, never on the longer past. That assumption is the entire model compressed into one line, the symbol $z_t$ joining the unified notation table of Appendix A on first use:

The three parameter blocks are the initial distribution $\pi_k = p(z_1 = k)$, the transition matrix $A_{jk} = p(z_t = k \mid z_{t-1} = j)$ (a $K \times K$ stochastic matrix whose rows sum to one), and the emission distribution $b_k(x) = p(x_t = x \mid z_t = k)$, which says how each regime generates observations. Write the bundle $\lambda = (\pi, A, B)$. The joint probability of a full hidden path and its observations factorizes along the chain,

which is exactly the discrete analogue of the linear-Gaussian factorization that drove the Kalman recursion. The structural parallel is worth stating as a table in prose: the Kalman state-transition matrix $\mathbf{F}$ becomes the transition table $A$, the Gaussian process noise becomes the stochasticity already baked into $A$, the observation matrix $\mathbf{H}$ and measurement noise become the emission distribution $b_k$, and the integral that marginalizes a continuous state becomes a sum over $K$ discrete states.

Where do discrete hidden states actually arise? Three canonical domains recur throughout this book. In finance, a market alternates between a calm low-volatility regime and a turbulent high-volatility one, and the regime label is never printed on the tape; you infer it from returns, which is the running example of subsection five. In speech, the hidden state is the phoneme being uttered and the observation is the acoustic frame, the application that made HMMs the backbone of recognition systems for two decades. In computational biology, the hidden state labels a region of a DNA sequence (coding versus non-coding, for instance) and the observation is the nucleotide, so that decoding the hidden path annotates the genome. In every case the regime is latent, the emission is observed, and the Markov assumption on the regime is the modeling bet.

2. The Three Classic Problems Intermediate

Rabiner's canonical treatment organizes everything you do with an HMM into three problems, and each is solved by a dynamic program that reuses subresults rather than enumerating the exponentially many hidden paths.

Problem 1, evaluation. Given parameters $\lambda$ and observations $x_{1:T}$, compute the likelihood $p(x_{1:T} \mid \lambda)$. The naive route sums the joint over all $K^T$ hidden paths, which is hopeless. The forward algorithm collapses that sum by caching the forward variable $\alpha_t(k) = p(x_{1:t}, z_t = k)$, the joint probability of the observations seen so far and being in state $k$ now. It obeys a one-line recursion that mirrors Kalman predict-then-update exactly:

The bracketed sum is the prediction step (propagate last step's belief through the transition matrix), and multiplying by $b_k(x_t)$ is the update step (weight by how well state $k$ explains the new observation). The total likelihood is the final column summed, $p(x_{1:T} \mid \lambda) = \sum_k \alpha_T(k)$. Cost is $O(K^2 T)$, not $O(K^T)$.

Problem 2, decoding. Find the single most likely hidden path $z_{1:T}^\star = \arg\max_{z_{1:T}} p(z_{1:T} \mid x_{1:T})$. Replacing the sum in the forward recursion with a maximum gives the Viterbi algorithm, whose variable $\delta_t(k)$ is the probability of the best path ending in state $k$ at time $t$:

The backpointer table $\psi$ remembers which predecessor won, and a backward pass from $\arg\max_k \delta_T(k)$ traces out the optimal path. Viterbi is to the HMM what the most-probable trajectory is to a state-space model; the swap of $\sum$ for $\max$ is the only difference between scoring all paths and finding the best one.

Problem 3, learning. Given only observations, estimate $\lambda$. There is no labeled hidden path, so we cannot just count transitions. The Baum-Welch algorithm, an instance of Expectation-Maximization, alternates an E-step that computes the expected number of times each transition and each emission occurred (using the forward variables above and matching backward variables $\beta_t(k) = p(x_{t+1:T} \mid z_t = k)$) with an M-step that re-estimates $\pi$, $A$, and $B$ from those expected counts. Each iteration provably never decreases the likelihood, climbing to a local optimum.

The E-step needs two posteriors. The first is the smoothed state posterior, the probability of being in state $k$ at time $t$ given the entire sequence, which combines the forward and backward variables:

This $\gamma_t$ is the discrete cousin of the Kalman smoother of Section 7.2: it uses future observations as well as past ones to refine the belief, which is why it beats the forward-only $\alpha_t$ for offline analysis. The second is the smoothed transition posterior, the probability of being in state $j$ at time $t$ and $k$ at time $t+1$:

The M-step then re-estimates each parameter as an expected-count ratio. The new initial distribution is just the smoothed posterior at $t=1$, $\hat{\pi}_k = \gamma_1(k)$. The new transition matrix is expected transitions $j \to k$ over expected visits to $j$, and the new Gaussian emission parameters are posterior-weighted moments of the data:

Read those formulas as soft counting: if every $\gamma_t(k)$ were a hard zero or one, they would reduce to the ordinary maximum-likelihood estimates you would write down from a labeled path, dividing observed transitions by observed visits and averaging the data in each state. Baum-Welch simply replaces the missing labels with their expected values under the current parameters and iterates, which is exactly what EM does for any latent-variable model.

3. Emissions: Categorical, Gaussian, and the Mixture Connection Intermediate

The transition machinery is fixed, but the emission distribution $b_k$ is a modeling choice that adapts the HMM to the data type. Two cases dominate. When observations are discrete symbols (nucleotides, quantized words, weather categories), each state carries a categorical emission, a length-$M$ probability vector $b_k(x = m) = B_{km}$ over the $M$ possible symbols. This is the classic discrete HMM, and Baum-Welch re-estimates $B$ by counting expected symbol emissions per state. When observations are continuous (financial returns, acoustic features, sensor readings), the natural choice is a Gaussian emission, $b_k(x) = \mathcal{N}(x \mid \mu_k, \Sigma_k)$, giving the Gaussian HMM in which each regime is a Gaussian cloud and learning re-estimates each $(\mu_k, \Sigma_k)$ as the posterior-weighted mean and covariance of the data.

A single Gaussian per state is sometimes too rigid: a phoneme's acoustic frames or a market regime's returns may be multimodal. The fix is to let each state emit from a mixture of Gaussians, $b_k(x) = \sum_{m=1}^{M} c_{km} \, \mathcal{N}(x \mid \mu_{km}, \Sigma_{km})$, the GMM-HMM that powered speech recognition before deep learning. This exposes a clean conceptual link. A plain Gaussian mixture model is exactly an HMM with $T = 1$, or equivalently an HMM whose transition matrix has identical rows so the state at each step is drawn independently, forgetting the past entirely. Turn that independence into a Markov chain and the mixture model grows a memory: the GMM clusters points, the HMM clusters points and models how the cluster assignment evolves over time. The HMM is a mixture model that remembers which component it was in last.

4. Regime-Switching Models, Limits, and What Comes After Advanced

Bolt an HMM on top of a forecasting model and you get a regime-switching model. The Markov-switching autoregression of Hamilton (1989) is the canonical case: the observation follows an AR process whose coefficients (and noise variance) depend on the current hidden regime,

so a single series can have, say, a slow mean-reverting calm regime and a violent high-variance crisis regime, with the hidden chain governing when each set of AR dynamics applies. This is the ARIMA family of Chapter 5 made regime-aware, and it is the workhorse of empirical macroeconomics and quantitative finance for exactly the calm-versus-crisis structure of subsection five.

The limitations are real and worth naming, because they are precisely what the rest of the book learns to relax. The number of states $K$ is fixed in advance and must be chosen by hand or by an information criterion; the model cannot grow a new regime it has never budgeted for. The first-order Markov assumption forces the state duration to be geometrically distributed, which often mismatches reality (a market crisis rarely ends with constant per-day probability). The emissions are simple parametric families, so the HMM cannot capture rich, high-dimensional observations such as raw audio or images without a great deal of feature engineering. And the discrete state is a severe information bottleneck: $K$ regimes carry only $\log_2 K$ bits of memory about the past.

Two of these limitations have classical patches worth knowing, because they show where the model bends before it breaks. The geometric-duration problem is fixed by the hidden semi-Markov model (HSMM), which augments each state with an explicit duration distribution so a regime can be made to last, on average, a realistic number of steps rather than decaying geometrically; the price is a heavier inference recursion that tracks how long the chain has dwelt in the current state. The other patch abandons the generative story entirely: a conditional random field models $p(z_{1:T} \mid x_{1:T})$ directly with the same Viterbi-style decoding but arbitrary, overlapping features of the observations, trading the HMM's ability to generate data for far richer emission conditioning. Both patches keep the dynamic-programming backbone; both were eventually subsumed, for high-dimensional data, by the neural sequence models of Part III. Knowing them sharpens your sense of which assumption to relax first when a plain HMM disappoints: duration, emission richness, or the Markov order itself.

5. Worked Example: Recovering Hidden Regimes in a Finance Series Advanced

We now build the whole pipeline on the running finance dataset of Part II: a returns series that secretly switches between a calm regime and a turbulent one. First we synthesize data from a known two-state Gaussian HMM so we have ground-truth regimes to grade against. Code 7.5.1 draws the hidden chain and its emissions.

import numpy as np

rng = np.random.default_rng(7)

# Ground-truth two-regime model: calm (low vol) and turbulent (high vol).
pi_true = np.array([0.8, 0.2])                 # start mostly calm
A_true = np.array([[0.97, 0.03],              # calm is sticky
                   [0.10, 0.90]])             # turbulent is sticky too
mu_true = np.array([0.0005, -0.0010])          # daily mean return per regime
sigma_true = np.array([0.006, 0.022])          # daily volatility per regime

def sample_hmm(T):
    """Draw a hidden Markov chain and its Gaussian (return) emissions."""
    z = np.empty(T, dtype=int)
    x = np.empty(T)
    z[0] = rng.choice(2, p=pi_true)
    for t in range(1, T):
        z[t] = rng.choice(2, p=A_true[z[t - 1]])   # Markov transition
    x = rng.normal(mu_true[z], sigma_true[z])       # regime-dependent emission
    return z, x

T = 1500
z_true, returns = sample_hmm(T)
print("days turbulent (true):", int(z_true.sum()), "of", T)
print("first 8 returns:", np.round(returns[:8], 4))

days turbulent (true): 271 of 1500
first 8 returns: [ 0.0049  0.0009 -0.0007  0.0083  0.0026 -0.0011  0.0046  0.0024]

Now the from-scratch core. Code 7.5.2 implements the forward algorithm (likelihood) and Viterbi decoding (most likely regime path) for a Gaussian HMM, working in log space to avoid the underflow that destroys a naive product of $T$ tiny probabilities. These are the recursions of subsection two written literally.

from scipy.stats import norm
from scipy.special import logsumexp

def gaussian_log_b(x, mu, sigma):
    """T x K matrix of log emission likelihoods log b_k(x_t)."""
    return np.stack([norm.logpdf(x, mu[k], sigma[k]) for k in range(len(mu))],
                    axis=1)

def forward_loglik(x, pi, A, mu, sigma):
    """Forward algorithm in log space; returns total log-likelihood."""
    logB = gaussian_log_b(x, mu, sigma)
    logA = np.log(A)
    log_alpha = np.log(pi) + logB[0]                       # alpha_1
    for t in range(1, len(x)):
        # predict (logsumexp over previous states) then update (+ emission)
        log_alpha = logsumexp(log_alpha[:, None] + logA, axis=0) + logB[t]
    return logsumexp(log_alpha)                            # sum over final states

def viterbi(x, pi, A, mu, sigma):
    """Most likely hidden state path (max-product dynamic program)."""
    logB = gaussian_log_b(x, mu, sigma)
    logA = np.log(A)
    T, K = logB.shape
    delta = np.log(pi) + logB[0]
    psi = np.zeros((T, K), dtype=int)
    for t in range(1, T):
        scores = delta[:, None] + logA          # K(prev) x K(next)
        psi[t] = np.argmax(scores, axis=0)       # backpointer
        delta = scores.max(axis=0) + logB[t]
    path = np.empty(T, dtype=int)
    path[-1] = int(np.argmax(delta))
    for t in range(T - 2, -1, -1):               # trace backpointers
        path[t] = psi[t + 1, path[t + 1]]
    return path

ll = forward_loglik(returns, pi_true, A_true, mu_true, sigma_true)
path_oracle = viterbi(returns, pi_true, A_true, mu_true, sigma_true)
acc = (path_oracle == z_true).mean()
print(f"log-likelihood at true params: {ll:.1f}")
print(f"Viterbi regime accuracy (true params): {acc:.3f}")

log-likelihood at true params: 5413.7
Viterbi regime accuracy (true params): 0.946

So far we supplied the true parameters. The learning problem (Problem 3) is to recover them from the returns alone. Code 7.5.3 implements one Baum-Welch iteration from scratch, the E-step and M-step formulas of subsection two written literally in log space, so you can watch EM climb the likelihood. The backward variable reuses the same log-space pattern as the forward pass.

def backward_loglik(logB, logA):
    """Backward variable beta in log space: log p(x_{t+1:T} | z_t = k)."""
    T, K = logB.shape
    log_beta = np.zeros((T, K))                 # beta_T(k) = 1 -> log 0
    for t in range(T - 2, -1, -1):
        log_beta[t] = logsumexp(logA + logB[t + 1] + log_beta[t + 1], axis=1)
    return log_beta

def baum_welch(x, pi, A, mu, sigma, n_iter=30):
    """Fit a Gaussian HMM from scratch by EM (Baum-Welch), log space."""
    pi, A, mu, sigma = map(np.array, (pi, A, mu, sigma))
    for _ in range(n_iter):
        logB = gaussian_log_b(x, mu, sigma)
        logA = np.log(A)
        T, K = logB.shape
        # E-step: forward, backward, then gamma and xi posteriors.
        log_alpha = np.empty((T, K)); log_alpha[0] = np.log(pi) + logB[0]
        for t in range(1, T):
            log_alpha[t] = logsumexp(log_alpha[t-1][:, None] + logA, 0) + logB[t]
        log_beta = backward_loglik(logB, logA)
        ll = logsumexp(log_alpha[-1])
        log_gamma = log_alpha + log_beta - ll
        gamma = np.exp(log_gamma)                              # T x K
        log_xi = (log_alpha[:-1, :, None] + logA[None]
                  + (logB[1:] + log_beta[1:])[:, None, :] - ll)
        xi = np.exp(log_xi).sum(0)                             # K x K expected
        # M-step: re-estimate pi, A, and per-state mean/variance.
        pi = gamma[0]
        A = xi / xi.sum(1, keepdims=True)
        w = gamma / gamma.sum(0, keepdims=True)
        mu = (w * x[:, None]).sum(0)
        sigma = np.sqrt((w * (x[:, None] - mu) ** 2).sum(0))
    return pi, A, mu, sigma, ll

# Start from deliberately wrong parameters and let EM find the regimes.
pi0, A0 = np.array([0.5, 0.5]), np.array([[0.5, 0.5], [0.5, 0.5]])
mu0, sigma0 = np.array([-0.01, 0.01]), np.array([0.01, 0.03])
pi_f, A_f, mu_f, sigma_f, ll_f = baum_welch(returns, pi0, A0, mu0, sigma0)
print("learned volatilities:", np.round(np.sort(sigma_f), 4))
print("learned transition diag:", np.round(np.diag(A_f), 3))
print(f"final log-likelihood: {ll_f:.1f}")

learned volatilities: [0.0059 0.0215]
learned transition diag: [0.969 0.901]
final log-likelihood: 5418.3

That from-scratch EM is instructive but verbose. In practice we reach for the library that implements the same algorithm, and state the line-count reduction explicitly. Code 7.5.4 fits and decodes with hmmlearn.

# pip install hmmlearn
from hmmlearn.hmm import GaussianHMM

model = GaussianHMM(n_components=2, covariance_type="diag",
                    n_iter=100, random_state=0)
model.fit(returns.reshape(-1, 1))            # Baum-Welch (EM) under the hood

learned_path = model.predict(returns.reshape(-1, 1))   # Viterbi decoding
learned_ll = model.score(returns.reshape(-1, 1))       # forward algorithm

# Label-switching is arbitrary in EM: name the wider-variance state "turbulent".
turbulent = int(np.argmax(model.covars_.ravel()))
pred_turbulent = (learned_path == turbulent).astype(int)
acc = (pred_turbulent == z_true).mean()
print(f"learned log-likelihood: {learned_ll:.1f}")
print(f"recovered volatilities: {np.sqrt(model.covars_.ravel())}")
print(f"regime accuracy vs ground truth: {acc:.3f}")

learned log-likelihood: 5418.2
recovered volatilities: [0.00592 0.02147]
regime accuracy vs ground truth: 0.939

The payoff is interpretive. The learned model hands back not just a regime label per day but the smoothed posterior $\gamma_t = p(z_t = \text{turbulent} \mid x_{1:T})$, a soft probability of turbulence on each day, available through model.predict_proba. Figure 7.5.1 sketches what that posterior looks like over a regime switch: calm returns hug zero with the posterior near zero, then a burst of large-magnitude returns drives the posterior sharply toward one, exactly the belief flip the numeric example of subsection three computed by hand.