State-Space Models & Kalman Smoothing

Model the unobserved level $\mu_t$ as a random walk and each reading as the level plus noise:

with all noise terms independent: a Gaussian random walk ( formalML: random-walks ) seen through measurement error. Here $q$ is the variance of the level’s yearly step and $r$ that of the measurement error (Durbin and Koopman’s $\sigma^2_\eta$ and $\sigma^2_\varepsilon$). The maximum-likelihood estimates (§36.8) are $\hat q = 1463.5474$ and $\hat r = 15108.3145$, so $\hat q/\hat r = 0.096870$, with diffuse log-likelihood $\hat\ell = -632.537686$ (Definition 5). Theorem 5 gives the distribution of each $\mu_t$ given all 100 readings; its mean falls from 1038.36 in 1897 to 951.00 in 1899 (Figure 1, left).

Difference the observations: $\Delta X_t = X_t - X_{t-1} = W_t + V_t - V_{t-1}$, which is stationary, with autocovariances

An autocovariance that cuts off after lag 1 is the signature of an MA(1) (Topic 35 §35.5), with $\rho(1) = -r/(q + 2r)$. So the hidden-level model is, from the observation-driven side, an ARIMA(0,1,1). At the MLE, $\rho(1) = -0.476901$; the sample lag-1 autocorrelation of the 99 differences is $-0.402043$, far outside the white-noise band $\pm 0.197$. The MA coefficient is $\theta = K_\infty - 1 = -0.733449$, where $K_\infty = 0.266551$ is the steady-state Kalman gain (§36.9, Example 15); as a check, $\theta/(1 + \theta^2) = -0.476901$.

We write $\mathbb{E}$, $\mathrm{Var}$ and $\mathrm{Cov}$ for expectation, variance and covariance, and $\mathcal{N}(x; \mu, \Sigma)$ for the Normal density with mean $\mu$ and covariance $\Sigma$ at $x$ (not Topic 34’s $\varphi$: $\phi$ is Topic 35’s AR coefficient). Bold letters stack over time. $T$ is a state-space sample length and $n$ an ARMA series length; they coincide in §36.4. Clashes are settled where they arise: the Kalman innovation $e_t$ versus Wold’s $\varepsilon_t$, the innovation variance $S_t$ versus Topic 35’s $v_{t-1}$ (equal, by Theorem 3), the covariance $P_t$ versus Topic 35’s projection $P_n$ (here $\Pi_n$), and $\vartheta$ for the parameter vector, because $\theta$ is the MA coefficient.

Section 36.4 proves that the hidden-state and autocovariance descriptions coincide for every ARMA model, and §36.3 shows that the Kalman filter computes Topic 35’s prediction errors exactly and fast. Then every tool works on both: the smoother answers the 1913 question (§36.6), the exact likelihood costs $O(n)$ (§36.3), missing data and forecasts are free (§36.5), and EM learns the parameters (§36.8).

For data $x_{1:T}$ and the state $Z_t$:

• the predictive distribution is $p(z_t \mid x_{1:t-1}) = \mathcal{N}(m^-_t, P^-_t)$;
• the filtering distribution is $p(z_t \mid x_{1:t}) = \mathcal{N}(m_t, P_t)$;
• the smoothing distribution is $p(z_t \mid x_{1:T}) = \mathcal{N}(m^s_t, P^s_t)$, for $t \le T$;
• the $h$-step forecast distribution is $p(z_{T+h} \mid x_{1:T})$, for $h \ge 1$.

Filtering is what we believe now; smoothing is what we believe in hindsight. Gaussianity of the last two is Theorems 4–5.

A constant-velocity tracker with state (position, velocity) observes position, with step $\Delta = 1$:

Start from $m_{t-1} = (0, 1)^\top$ and $P_{t-1} = \begin{pmatrix} 0.6 & 0.2 \\ 0.2 & 0.4 \end{pmatrix}$.

• Predict. $F$ shears the ellipse (position picks up velocity), then $Q$ inflates it: $m^-_t = (1, 1)^\top$ and $P^-_t = \begin{pmatrix} 1.5 & 0.75 \\ 0.75 & 0.7 \end{pmatrix}$.
• Update. We observe position $x_t = 1.8$. Then $S_t = 1.9$, $K_t = (0.789474, 0.394737)^\top$, $m_t = (1.631579, 1.315789)^\top$, and $P_t = \begin{pmatrix} 0.315789 & 0.157895 \\ 0.157895 & 0.403947 \end{pmatrix}$.
• Area ($\propto \sqrt{\det P}$): 0.4472 → 0.6982 → 0.3204. Correlation: 0.408 → 0.732 → 0.442.

We observe only position, yet velocity’s variance shrinks too, from 0.7 to 0.404: the predict step correlated the coordinates, and the update slices the ellipse along the direction the observation sees, pulling the correlated direction along.

Take a regression $y_t = \mathbf{x}_t^\top\beta + V_t$, $V_t \sim \mathcal{N}(0, r)$, with examples arriving one at a time. Make the coefficients a state that never moves ($F = I$, $Q = 0$), with time-varying observation row $H_t = \mathbf{x}_t^\top$ (Remark 5) and prior $Z_1 \sim \mathcal{N}(0, \kappa I)$.

The filter mean after $t$ observations is the Bayesian linear regression posterior mean (Topic 25 §25.9), which is ridge regression with penalty $r/\kappa$; as $\kappa \to \infty$ it becomes ordinary least squares (Topic 21 §21.4). With $n = 200$ examples, three coefficients, $r = 0.25$ and $\kappa = 10^8$, the final filter mean differs from $\hat\beta_{\text{OLS}}$ by at most $1.5 \times 10^{-10}$, and $P_T/r$ matches $(X^\top X)^{-1}$ to $3.1 \times 10^{-12}$.

This is recursive least squares: online learning with exact uncertainty in $O(p^2)$ per example for $p$ coefficients, and no refit. Setting $Q \ne 0$ lets the coefficients drift, the Kalman answer to concept drift.

Topic 34 delivered the filter and the predictive likelihood terms; the rest is new. The smoother of §36.6 is the continuous-state twin of §34.5’s forward–backward, and §36.8’s EM is the twin of §34.7’s Baum–Welch.

The proofs use $P_t = P^-_t - K_tS_tK_t^\top$. Code should use the Joseph form (JOS1968),

The two are algebraically identical (agreement $4.5 \times 10^{-11}$ on the Nile), but the Joseph form is a sum of positive-semidefinite terms, so rounding cannot destroy that property. Square-root filters propagate a Cholesky factor of $P_t$; information filters propagate $P_t^{-1}$, which suits vague priors; UD filters factor $P_t = UDU^\top$ (SAR2013, AND1979).

Every result below holds verbatim with $F_t, H_t, Q_t, R_t$ in place of $F, H, Q, R$, since only §36.9 uses time invariance. Examples 4 and 16 need this.

Let $e_t = x_t - Hm^-_t$ and $S_t$ be the Kalman innovations and innovation covariances, with every $x_t$ observed. Stack $\mathbf{x} = (x_1, \dots, x_T)$ and let $\boldsymbol\mu_X = \mathbb{E}[\mathbf{X}]$ and $\Sigma_X = \mathrm{Cov}(\mathbf{X})$. Then:

(i) The innovations $e_1, \dots, e_T$ are independent, with $e_t \sim \mathcal{N}(0, S_t)$.

(ii) There is a unit lower block-triangular matrix $L$ (identity $d \times d$ blocks on the diagonal) with $\mathbf{x} - \boldsymbol\mu_X = L\mathbf{e}$. Hence $\Sigma_X = LDL^\top$ with $D = \mathrm{diag}(S_1, \dots, S_T)$.

(iii) The log-likelihood is

and the filter evaluates it in $O\bigl(T(k^3 + d^3)\bigr)$ operations.

Theorem 1 holds for any $P_1$. Starting the Nile filter from $P_1 = \kappa = 10^6$ (Definition 5), the $t = 1$ term ($-8.4521$) measures only the vague prior, so we drop it; at $(\hat q, \hat r)$ the diffuse log-likelihood is

• The innovation variance shrinks as the level is learned, from $S_2 = 31455.31$ to $S_{100} = 20599.00$.
• The innovations whiten the series. The lag-1 autocorrelation is 0.4984 for the flows and 0.1151 for the standardized innovations $e_t/\sqrt{S_t}$, and 4 of the 99 standardized innovations fall outside $\pm 1.96$.
• For an ARMA(1,1), the factorization $\Gamma = LDL^\top$ of Topic 35’s $12 \times 12$ Toeplitz $\Gamma_n$ returns the Kalman $S_t$ on the diagonal of $D$, to $1.1 \times 10^{-15}$.

Evaluating $\mathcal{N}(\mathbf{x}; \boldsymbol\mu_X, \Sigma_X)$ directly needs a Cholesky factor of a $Td \times Td$ matrix: $O(T^3d^3)$ time, $O(T^2d^2)$ memory. The filter needs $O(T(k^3 + d^3))$ time and $O(k^2)$ memory. For a 10,000-step sensor log that is infeasible versus instant; Example 16 shows the same gap for GP regression.

Under a correct model the standardized innovations are i.i.d. $\mathcal{N}(0, 1)$ (Theorem 1(i)): they are the residuals. Check them with the sample ACF and its white-noise band (Topic 35 §35.6), a Ljung–Box statistic, and a Normal Q-Q plot. On the Nile, the run of large values after 1898 flags the level shift as an intervention (Remark 18). This is the residual check for any learned state-space or sequence model.

Let $\phi(B)X_t = \theta(B)\varepsilon_t$ be an ARMA($p$, $q$) with $\varepsilon_t$ i.i.d. $\mathcal{N}(0, \sigma^2)$. Let $k = \max(p, q + 1)$, and pad with $\phi_j = 0$ for $j > p$ and $\theta_j = 0$ for $j > q$.

• $F$ is $k \times k$, with first column $(\phi_1, \dots, \phi_k)^\top$, ones on the superdiagonal, and zeros elsewhere.
• $\mathbf{g} = (1, \theta_1, \dots, \theta_{k-1})^\top$.
• $Z_t = FZ_{t-1} + \mathbf{g}\varepsilon_t$ and $X_t = HZ_t$, with $H = (1, 0, \dots, 0)$.
• So $Q = \sigma^2\mathbf{g}\mathbf{g}^\top$ (rank 1) and $R = 0$: an observation with no noise.
• $P_1$ is the stationary covariance, the solution of $P = FPF^\top + Q$.

For an ARMA(1,1), $k = 2$ and

This is a noise-driven linear difference system ( formalCalculus: linear-systems ; HAR1989).

Let $\phi(B)X_t = \theta(B)\varepsilon_t$ be an ARMA($p$, $q$) with $\varepsilon_t$ i.i.d. $\mathcal{N}(0, \sigma^2)$, and take $F$, $\mathbf{g}$, $H$ from Definition 2, with $k = \max(p, q+1)$. Then:

(i) If $Z_t = FZ_{t-1} + \mathbf{g}\varepsilon_t$ and $X_t = HZ_t$, the process $X_t$ satisfies $X_t = \sum_{j=1}^{p}\phi_jX_{t-j} + \varepsilon_t + \sum_{j=1}^{q}\theta_j\varepsilon_{t-j}$.

(ii) A nonzero $\lambda$ is an eigenvalue of $F$ if and only if $\phi(1/\lambda) = 0$. So the process is causal if and only if the spectral radius $\rho(F)$ is below 1.

(iii) If the model is causal, $P = FPF^\top + Q$ (with $Q = \sigma^2\mathbf{g}\mathbf{g}^\top$) has the unique solution $P_1 = \sum_{j\ge 0}F^jQF^{j\top}$. Started from $Z_1 \sim \mathcal{N}(0, P_1)$, the observations $(X_1, \dots, X_n)$ have the law of $n$ consecutive values of the stationary ARMA, with autocovariance $\gamma(h) = HF^hP_1H^\top$ and $\psi_j = HF^j\mathbf{g}$.

Here $k = 1$, $F = \phi$, and the Lyapunov equation $P = \phi^2P + \sigma^2$ gives $P_1 = \sigma^2/(1 - \phi^2) = \gamma(0)$. Since $R = 0$, the first update pins the state to $x_1$, and the filter then predicts $\phi x_{t-1}$ with variance $\sigma^2$. So $S_1 = \gamma(0)$ and $S_t = \sigma^2$ for $t \ge 2$: the exact likelihood is the conditional likelihood plus one term for $x_1$. For $\phi = 0.7$, $S_1/\sigma^2 = 1.960784$.

Let $\{X_t\}$ be a causal, mean-zero ARMA($p$, $q$) driven by i.i.d. $\mathcal{N}(0, \sigma^2)$ noise, filtered in the form of Definition 2 with the stationary $P_1$ of Theorem 2(iii). Write $\hat X_t = \Pi_{t-1}X_t$ for the best linear predictor of $X_t$ from $X_1, \dots, X_{t-1}$, and $v_{t-1} = \mathbb{E}(X_t - \hat X_t)^2$ (Topic 35 Def 5, with $\hat X_1 = 0$ and $v_0 = \gamma(0)$). Topic 35 wrote this projection $P_{t-1}X_t$; here $P_t$ is a covariance, so the projection is $\Pi_{t-1}$. Then for every $t \ge 1$:

(i) $Hm^-_t = \hat x_t$, so the Kalman innovation is the prediction error: $e_t = x_t - \hat x_t$.

(ii) $S_t = v_{t-1}$.

(iii) If the model is also invertible, then $S_t \ge \sigma^2$, $S_t \to \sigma^2$, and $\mathbb{E}(e_t - \varepsilon_t)^2 = S_t - \sigma^2 \to 0$, where $e_t$ denotes the random innovation $X_t - \hat X_t$.

Topic 35’s ARMA(1,1) fixture ($\phi = 0.7$, $\theta = 0.4$, $n = 1000$), demeaned, with every estimate evaluated under the exact likelihood:

At $n = 1000$ all three are within one standard error of each other; Theorem 8’s asymptotic equivalence shows as contour-level agreement in Figure 4 (left). At the true parameters, the ratios $S_t/\sigma^2$ for $t = 1, \dots, 6$ are 3.372549, 1.112558, 1.016187, 1.002549, 1.000407 and 1.000065, converging at the rate $\theta^2 = 0.16$ per step (Figure 4, right). See mode B of the explorer in §36.8.

Exact Gaussian MLE, mean estimated, on Yule’s sunspots (Topic 35 §35.6); the exact Kalman log-likelihood agrees with SARIMAX to $10^{-9}$:

CSS gives AR(2) $\hat\phi = (1.391812, -0.690282)$ and ARMA(2,1) $(1.472013, -0.756875;\ -0.154109)$; Whittle gives $(1.393809, -0.688528)$ and $(1.469596, -0.751783;\ -0.144834)$.

• AR(2) innovation variances: $S_1/\sigma^2 = 5.910522$ and $S_2/\sigma^2 = 1.901612$, then exactly 1 from $t = 3$ on. CSS throws away precisely these two start-up terms (Remark 10).
• AIC vs BIC: with $n_\vartheta = p + q + 2$ parameters, AIC prefers ARMA(2,1) (2620.28 vs 2622.64) and BIC prefers AR(2) (2637.57 vs 2638.94).
• Reconciliation: Topic 35’s table printed AR(2) $\hat\ell = -1307.319$ from a slightly early-stopped optimizer; the exact optimum is $-1307.318169$.

This is the exact estimator that Topic 35’s Remark 14 deferred.

For invertible $M$, $(MFM^{-1}, HM^{-1}, MQM^\top)$ gives the same observation process. Harvey’s, the companion and Akaike’s forms are similarity transforms of one another, so $F$, $H$, $Q$ are unidentifiable without a canonical form (Remark 20).

Because $P_1$ and $Q$ scale with $\sigma^2$, run the filter at $\sigma^2 = 1$ to get $e_t$ and $s_t = S_t/\sigma^2$. Then $\hat\sigma^2 = n^{-1}\sum_t e_t^2/s_t$ and

This is the surface of Figure 4 and mode B, with one parameter fewer for Nelder–Mead.

CSS conditions on the first $p$ values and sets pre-sample shocks to zero. The exact likelihood keeps the start-up terms with $S_t > \sigma^2$: $O(1)$ in total, but they matter when $n$ is small or roots are near the unit circle. They vanish geometrically, so all three estimators share a limit (Topic 35 Thm 8).

A unit root makes $\rho(F) = 1$, and the Lyapunov equation has no solution. State-space methods handle ARIMA by differencing first or by a diffuse prior on the nonstationary states (Definition 5, §36.7, Remark 17).

Write $\mathcal{O} \subseteq \{1, \dots, T\}$ for the set of times observed and $x_{\mathcal{O}}$ for the observed values. At a time $t \notin \mathcal{O}$ the filter skips the update, setting $m_t = m^-_t$ and $P_t = P^-_t$, and adds no likelihood term. In code, a missing value is NaN.

Let $\mathcal{O} \subseteq \{1, \dots, T\}$ be the observed times, and run the filter of Definition 3: at $t \notin \mathcal{O}$, set $m_t = m^-_t$ and $P_t = P^-_t$ and add no likelihood term. Then for every $t$: