47. The Kalman Filter and Vector Autoregressions#

47.1. Overview#

This lecture derives the Kalman filter for a linear Gaussian state space system and then uses it to construct vector autoregressions (VARs).

It builds on A First Look at the Kalman Filter, where the filter is introduced through filtering and forecasting distributions.

It complements Another Look at the Kalman Filter, where the same recursion is used in a worker-firm learning problem.

This lecture shifts attention from filtering a hidden state to representing and estimating the observable process generated by that hidden state.

Our approach rests on repeated applications of the population linear least squares projection formula, the insight that computing a conditional expectation of a jointly normal random vector is the same as running a population OLS regression.

The lecture covers:

  • deriving the Kalman filter recursions from first principles

  • the matrix Riccati difference equation governing conditional covariance matrices

  • the innovations representation and the Gram-Schmidt whitening property

  • the structure of a hidden Markov model

  • the likelihood function for a state space system and its role in maximum likelihood and Bayesian estimation

  • how the time-invariant Kalman filter generates a vector autoregression

  • why the Kalman filter is an essential tool for interpreting VARs estimated from economic data

47.2. The state space system#

The Kalman filter applies to the state space system for \(t \geq 0\):

(47.1)#\[\begin{split} \begin{aligned} x_{t+1} &= A x_t + C w_{t+1} \\ y_t &= G x_t + v_t \end{aligned} \end{split}\]

where

  • \(x_t\) is an \(n \times 1\) state vector (hidden, unobserved)

  • \(y_t\) is an \(m \times 1\) vector of signals on the hidden state (observed)

  • \(w_{t+1}\) is a \(p \times 1\) IID sequence of normal random variables with mean \(0\) and identity covariance matrix

  • \(v_t\) is an IID sequence of normal random variables with mean zero and covariance matrix \(R\)

  • \(w_{t+1}\) and \(v_s\) are orthogonal for all \(t+1\) and \(s \geq 0\)

The coefficient matrices have the following dimensions: \(A\) is \(n \times n\), \(C\) is \(n \times p\), \(G\) is \(m \times n\), and \(R\) is \(m \times m\).

The initial state satisfies

(47.2)#\[ x_0 \sim N(\hat{x}_0, \Sigma_0) \]

We observe \(y_t, \ldots, y_0\) but not \(x_t, \ldots, x_0\) at time \(t\), and we know all first and second moments implied by (47.1) and (47.2).

47.3. The Kalman filter#

47.3.1. Starting distribution#

Working forward in time starting at \(t = 0\), before observing \(y_0\), the specification (47.1)-(47.2) implies that the marginal distribution of \(y_0\) is

(47.3)#\[ y_0 \sim N(G \hat{x}_0,\; G \Sigma_0 G^\top + R) \]

For \(t \geq 0\) let \(y^t = [y_t, y_{t-1}, \ldots, y_0]\).

We want a convenient recursive representation of the conditional distribution of \(y_t\) given history \(y^{t-1}\).

The Kalman filter attains this by constructing recursive formulas for \(\hat{x}_t\) and \(\Sigma_t\) such that the distribution of \(y_t\) conditional on \(y^{t-1}\) generalises (47.3) to

(47.4)#\[ y_t \sim N(G \hat{x}_t,\; G \Sigma_t G^\top + R) \]

for \(t \geq 1\), where the distribution of \(x_t\) conditional on \(y^{t-1}\) is \(N(\hat{x}_t, \Sigma_t)\).

The objects \(\hat{x}_t\) and \(\Sigma_t\) characterise the population regression

\[ \hat{x}_t = \mathbb{E}[x_t \mid y_{t-1}, \ldots, y_0] \]

and the conditional covariance matrix

\[ \Sigma_t = \mathbb{E}\!\left[(x_t - \hat{x}_t)(x_t - \hat{x}_t)^\top \mid y_{t-1}, \ldots, y_0\right] \]

47.3.2. Derivation#

At each date, our approach is to regress what we do not know on what we know.

Note

Because our assumptions imply that \(\{x_t, y_t\}_{t=0}^\infty\) is a jointly normal stochastic process, linear least squares regressions equal conditional mathematical expectations.

Each step below is an application of Bayes’ law.

Under the weaker assumption that all means and covariances exist without joint normality, the same calculations yield “wide-sense conditional expectations” that coincide with true conditional expectations only when those conditional expectations are linear.

We arrive at \(t = 0\) knowing \(\hat{x}_0\) and \(\Sigma_0\).

The information about \(x_0\) in \(y_0\) that is new relative to \((\hat{x}_0, \Sigma_0)\) is the innovation

\[ a_0 \equiv y_0 - G \hat{x}_0 \]

Let \(L_0\) be the population regression coefficient of the hidden-state error \(x_0-\hat{x}_0\) on the signal surprise \(y_0-G\hat{x}_0\).

The conditional mean \(\mathbb{E}[x_0 \mid y_0] = \hat{x}_0 + L_0(y_0 - G\hat{x}_0)\) satisfies the population regression formula

(47.5)#\[ x_0 - \hat{x}_0 = L_0(y_0 - G\hat{x}_0) + \eta \]

where \(\eta\) is the least squares residual.

Orthogonality of \(\eta\) to \((y_0 - G\hat{x}_0)\) pins down \(L_0\) via the normal equations

\[ \mathbb{E}(x_0 - \hat{x}_0)(y_0 - G\hat{x}_0)^\top = L_0\, \mathbb{E}(y_0 - G\hat{x}_0)(y_0 - G\hat{x}_0)^\top \]

Evaluating the moment matrices and solving for \(L_0\) gives

(47.6)#\[ L_0 = \Sigma_0 G^\top(G \Sigma_0 G^\top + R)^{-1} \]

Thus \(L_0\) updates the estimate of \(x_0\), while \(K_0=A L_0\) updates the forecast of \(x_1\).

To forecast \(x_1\), note that

(47.7)#\[ x_1 = A\hat{x}_0 + A(x_0 - \hat{x}_0) + C w_1 \]

Applying (47.5) gives \(\mathbb{E}[x_1 \mid y_0] = A\hat{x}_0 + AL_0(y_0 - G\hat{x}_0)\), which we write as

(47.8)#\[ \hat{x}_1 = A\hat{x}_0 + K_0(y_0 - G\hat{x}_0) \]

where the Kalman gain at time 0 is

(47.9)#\[ K_0 = A \Sigma_0 G^\top(G \Sigma_0 G^\top + R)^{-1} \]

Subtracting (47.8) from (47.7) yields

(47.10)#\[ x_1 - \hat{x}_1 = A(x_0 - \hat{x}_0) + C w_1 - K_0(y_0 - G\hat{x}_0) \]

Using (47.10) and \(y_0 = G x_0 + v_0\) to evaluate \(\Sigma_1 \equiv \mathbb{E}[(x_1 - \hat{x}_1)(x_1 - \hat{x}_1)^\top \mid y_0]\) gives

(47.11)#\[ \Sigma_1 = (A - K_0 G)\Sigma_0(A - K_0 G)^\top + CC^\top + K_0 R K_0^\top \]

Thus \(f(x_1 \mid y_0) \sim N(\hat{x}_1, \Sigma_1)\).

Collecting the time-\(0\) equations:

(47.12)#\[\begin{split} \begin{aligned} a_0 &= y_0 - G\hat{x}_0 \\ K_0 &= A\Sigma_0 G^\top(G\Sigma_0 G^\top + R)^{-1} \\ \hat{x}_1 &= A\hat{x}_0 + K_0 a_0 \\ \Sigma_1 &= CC^\top + K_0 R K_0^\top + (A - K_0 G)\Sigma_0(A - K_0 G)^\top \end{aligned} \end{split}\]

System (47.12) maps a mean-covariance pair \((\hat{x}_0, \Sigma_0)\) into a new pair \((\hat{x}_1, \Sigma_1)\), with auxiliary outputs \((a_0, K_0)\).

Recognising that “we are in the same situation at the start of period 1 as at the start of period 0” activates a recursion, the Kalman filter.

47.3.3. The Kalman filter recursions#

Iterating system (47.12) yields the Kalman filter for \(t \geq 0\):

(47.13)#\[\begin{split} \begin{aligned} a_t &= y_t - G\hat{x}_t \\ K_t &= A\Sigma_t G^\top(G\Sigma_t G^\top + R)^{-1} \\ \hat{x}_{t+1} &= A\hat{x}_t + K_t a_t \\ \Sigma_{t+1} &= CC^\top + K_t R K_t^\top + (A - K_t G)\Sigma_t(A - K_t G)^\top \end{aligned} \end{split}\]

Here \(K_t\) is the Kalman gain at time \(t\).

47.3.4. The matrix Riccati equation#

Substituting the expression for \(K_t\) from the second line of (47.13) into the fourth line gives an equivalent update formula:

(47.14)#\[ \Sigma_{t+1} = A\Sigma_t A^\top + CC^\top - A\Sigma_t G^\top(G\Sigma_t G^\top + R)^{-1} G\Sigma_t A^\top \]

Equation (47.14) is the matrix Riccati difference equation.

It governs the sequence of conditional covariance matrices \(\{\Sigma_t\}_{t=0}^\infty\) without reference to the observations \(\{y_t\}\).

47.4. The Gram-Schmidt process#

The random vector

\[ a_t = y_t - \mathbb{E}[y_t \mid y_{t-1}, \ldots, y_0] \]

is the innovation of \(y_t\) with respect to \(y^{t-1}\), the part of \(y_t\) that cannot be predicted from past observations.

Note that \(\mathbb{E} a_t a_t^\top = G\Sigma_t G^\top + R\), the matrix whose inverse appears in the Kalman gain formula (47.13).

A direct calculation using \(a_t = G(x_t - \hat{x}_t) + v_t\) shows that \(\mathbb{E} a_t a_{t-1}^\top = 0\) and, more generally, \(\mathbb{E}[a_t \mid a_{t-1}, \ldots, a_0] = 0\).

Note

An alternative argument from first principles: let \(H(y^t)\) denote the closed linear span of \(y^t\).

Since \(a_{t+1} = y_{t+1} - \mathbb{E}[y_{t+1} \mid y^t]\) is a least-squares error, \(a_{t+1} \perp H(y^t)\), and in particular \(a_{t+1} \perp a_t\).

Thus \(\{a_t\}\) is a white-noise process of innovations to \(\{y_t\}\).

Sometimes (47.13) is called a whitening filter: it takes the signal process \(\{y_t\}\) as input and produces the white-noise innovation process \(\{a_t\}\) as output.

With \(H(a^t)\) defined analogously, the linear space \(H(a^t)\) is an orthogonal basis for the linear space \(H(y^t)\).

Rather than computing \(\mathbb{E}[x_t \mid y_{t-1}, \ldots, y_0]\) via one large regression, the Kalman filter performs a sequence of small regressions on successive orthogonal components of the basis \([a_{t-1}, \ldots, a_0]\), an instance of the Gram-Schmidt procedure.

47.5. Hidden Markov model#

System (47.1)-(47.2) is an example of a hidden Markov model.

The observed process \(\{y_t\}_{t=0}^\infty\) is not Markov, but the hidden process \(\{x_t\}_{t=0}^\infty\) is Markov.

So is the process of means and covariances \(\{(\hat{x}_t, \Sigma_t)\}\), which are sufficient statistics for the distribution of \(x_t\) conditional on \([y_{t-1}, \ldots, y_0]\).

47.6. Estimation#

47.6.1. The innovations representation#

The innovations representation emerging from the Kalman filter is

(47.15)#\[\begin{split} \begin{aligned} \hat{x}_{t+1} &= A\hat{x}_t + K_t a_t \\ y_t &= G\hat{x}_t + a_t \end{aligned} \end{split}\]

where \(\hat{x}_t = \mathbb{E}[x_t \mid y^{t-1}]\) for \(t \geq 1\) and \(\mathbb{E}[a_t a_t^\top \mid y^{t-1}] = G\Sigma_t G^\top + R \equiv \Omega_t\).

For \(t \geq 1\), \(\mathbb{E}[y_t \mid y^{t-1}] = G\hat{x}_t\) and the conditional distribution of \(y_t\) given \(y^{t-1}\) is \(N(G\hat{x}_t, \Omega_t)\).

The objects \((G\hat{x}_t, \Omega_t)\) emerging from the Kalman filter recursions therefore completely characterise this conditional distribution.

47.6.2. The likelihood function#

We can factor the likelihood of a sample \((y_T, y_{T-1}, \ldots, y_0)\) as

(47.16)#\[ f(y_T, \ldots, y_0) = f(y_T \mid y^{T-1})\, f(y_{T-1} \mid y^{T-2}) \cdots f(y_1 \mid y_0)\, f(y_0) \]

The log conditional density of the \(m \times 1\) vector \(y_t\) is

(47.17)#\[ \log f(y_t \mid y^{t-1}) = -\frac{m}{2}\log(2\pi) - \frac{1}{2}\log\det(\Omega_t) - \frac{1}{2}\, a_t^\top \Omega_t^{-1} a_t \]

Using (47.17) and (47.13) together, we can evaluate the likelihood (47.16) recursively for any parameter vector \(\theta\) that underlies the matrices \(A, G, C, R\).

Such calculations are at the heart of efficient strategies for computing maximum likelihood estimators of free parameters.

47.6.3. Bayesian inference#

The likelihood function is also central to Bayesian inference.

Where \(\theta\) is the parameter vector, \(y_0^T\) the data, and \(\tilde{p}(\theta)\) a prior density over \(\theta\) before seeing \(y_0^T\), Bayes’ law gives the posterior

\[ \tilde{p}(\theta \mid y_0^T) = \frac{f(y_0^T \mid \theta)\,\tilde{p}(\theta)} {\int f(y_0^T \mid \theta)\,\tilde{p}(\theta)\, d\theta} \]

The denominator is the marginal joint density \(f(y_0^T)\).

47.7. Vector autoregressions and the Kalman filter#

47.7.1. Convergence to a steady state#

Under conditions discussed by Anderson et al. [1996], iterations on the Riccati equation (47.14) converge to a time-invariant matrix \(\Sigma\) from any positive semi-definite starting value \(\Sigma_0\).

A time-invariant fixed point \(\Sigma_t = \Sigma\) of (47.14) is the covariance matrix of \(x_t\) around

\[ \mathbb{E}\!\left[x_t \mid \{y_s\}_{s \leq t-1}\right] \]

where the conditioning extends over the semi-infinite past \(s \leq t-1\).

47.7.2. A time-invariant VAR#

If the fixed point \(\Sigma\) exists and we initialise the filter at \(\Sigma_0 = \Sigma\), the innovations representation (47.15) becomes time-invariant:

(47.18)#\[\begin{split} \begin{aligned} \hat{x}_{t+1} &= A\hat{x}_t + K a_t \\ y_t &= G\hat{x}_t + a_t \end{aligned} \end{split}\]

where \(\mathbb{E} a_t a_t^\top = G\Sigma G^\top + R\) and the steady-state Kalman gain is \(K = A\Sigma G^\top(G\Sigma G^\top + R)^{-1}\).

From (47.18) we obtain \(\hat{x}_{t+1} = (A - KG)\hat{x}_t + K y_t\).

If the eigenvalues of \(A - KG\) are bounded in modulus strictly below unity, we can solve this equation forward to get

(47.19)#\[ \hat{x}_{t+1} = \sum_{j=0}^\infty (A - KG)^j K\, y_{t-j} \]

Substituting (47.19) into the observation equation of (47.18) gives the vector autoregression

(47.20)#\[ y_t = G \sum_{j=0}^\infty (A - KG)^j K\, y_{t-j-1} + a_t \]

where by construction

(47.21)#\[ \mathbb{E}\!\left[a_t\, y_{t-j-1}^\top\right] = 0 \quad \forall\, j \geq 0 \]

The orthogonality conditions (47.21) identify (47.20) as a vector autoregression.

Defining the lag operator \(L\) by \(L x_{t+1} \equiv x_t\), the moving average representation deduced from (47.18) is

\[ y_t = \left[I + G(I - AL)^{-1} KL\right] a_t = \left[I + G\sum_{j=0}^\infty A^j K L^{j+1}\right] a_t \]

47.7.3. Interpreting VARs#

Equilibria of economic models (or their linear or log-linear approximations) typically take the form of state space system (47.1).

This hidden Markov model disturbs the state \(x_t\) by the \(p \times 1\) shock vector \(w_{t+1}\) and perturbs the \(m \times 1\) vector of observables \(y_t\) by the \(m \times 1\) measurement error \(v_t\).

An economic theory typically makes \(w_{t+1}\) and \(v_t\) directly interpretable as shocks to preferences, technologies, endowments, or information sets.

The state space system (47.1) represents \(\{y_t\}\) in terms of these interpretable shocks.

However, in the typical situation these shocks cannot be recovered directly from the \(y_t\)’s, even when \(A, G, C, R\) are known.

The innovations representation (47.18) represents the same stochastic process \(\{y_t\}\) in terms of the \(m \times 1\) vector \(a_t\) of innovations that would be recovered by running an infinite-order population vector autoregression.

Its role in mapping the original representation (47.1) to the VAR (47.20) makes the Kalman filter an indispensable tool for interpreting vector autoregressions.

47.8. Spectral factorization identity#

Because the original state space system (47.1) and the innovations representation (47.18) describe the same stochastic process \(\{y_t\}\), they imply two distinct formulas for the spectral density matrix of \(\{y_t\}\).

Equating those formulas yields the spectral factorization identity.

47.8.1. Two representations of the spectral density#

Consider first the original state space system.

Writing the first line of (47.1) as \(x_t = (zI - A)^{-1} C w_{t+1}\) (using the \(z\)-transform convention \(z^{-1} x_t = x_{t-1}\)), the covariance generating function of \(\{x_t\}\) is

\[ S_x(z) = (zI - A)^{-1} CC^\top (z^{-1}I - A^\top)^{-1}. \]

Since \(v_t\) is orthogonal to \(x_t\), the spectral density of \(\{y_t\}\) is

(47.22)#\[ S_y(z) = G(zI - A)^{-1} CC^\top (z^{-1}I - A^\top)^{-1} G^\top + R. \]

Now consider the innovations representation.

The time-invariant innovations representation (47.18) gives \(y_t = [G(zI - A)^{-1}K + I]\, a_t\).

Since \(a_t\) is white noise with covariance matrix \(G\Sigma G^\top + R\), the spectral density is also

(47.23)#\[ S_y(z) = \bigl[G(zI-A)^{-1}K + I\bigr] \bigl(G\Sigma G^\top + R\bigr) \bigl[K^\top(z^{-1}I - A^\top)^{-1}G^\top + I\bigr]. \]

47.8.2. The spectral factorization identity#

Equating (47.22) and (47.23) gives the spectral factorization identity:

(47.24)#\[ G(zI - A)^{-1} CC^\top (z^{-1}I - A^\top)^{-1} G^\top + R = \bigl[G(zI-A)^{-1}K + I\bigr] \bigl(G\Sigma G^\top + R\bigr) \bigl[K^\top(z^{-1}I - A^\top)^{-1}G^\top + I\bigr]. \]

The left side expresses \(S_y(z)\) in terms of the structural shocks \((w_{t+1}, v_t)\) and the matrices \((A, C, G, R)\).

The right side expresses the same object as a spectral factor built from the innovations \(a_t\) and the steady-state Kalman gain \(K\).

47.8.3. Wold and autoregressive representations#

Starting from the innovations representation (47.18), we can obtain both a Wold moving average representation and an autoregressive representation.

For the Wold representation, iterate the state equation in (47.18) to express the current observation in terms of current and past innovations.

With \(L\) denoting the lag operator, (47.18) implies

(47.25)#\[ y_t = \left[I + G(I - AL)^{-1} K L\right] a_t \]

This is the Wold moving-average representation of \(\{y_t\}\) in terms of its one-step forecast errors.

For the autoregressive representation, invert the moving-average operator in (47.25) and solve for \(a_t\) in terms of current and past observations.

Using the identity

\[ \left[I + G(I - AL)^{-1} K L\right]^{-1} = I - G\left[I - (A - KG)L\right]^{-1} K L \]

gives

(47.26)#\[ y_t = G\bigl[I-(A-KG)L\bigr]^{-1}K\, y_{t-1} + a_t = \sum_{j=1}^\infty G(A-KG)^{j-1}K\, y_{t-j} + a_t, \]

which is the vector autoregression already stated in (47.20).

The key analytical fact is that, under mild stability conditions, the zeros of \(\det[G(zI-A)^{-1}K + I]\) all lie inside the unit circle.

This ensures that the moving-average operator in (47.25) has a causal one-sided inverse.

Hence \(a_t\) lies in the closed linear span of current and past observations \(y^t\), so \(a_t\) is the population forecast error in the VAR.

47.9. Python implementation#

We now illustrate the theory using the quantecon library, which provides LinearStateSpace and Kalman classes that implement everything derived above.

We use the following imports:

import numpy as np
import matplotlib.pyplot as plt
import quantecon as qe

47.9.1. A scalar hidden AR(1) model#

Consider a scalar hidden AR(1) state observed with measurement noise:

\[\begin{split} \begin{aligned} x_{t+1} &= \rho\, x_t + \sigma_w\, w_{t+1} \\ y_t &= x_t + \sigma_v\, v_t \end{aligned} \end{split}\]

with \(w_t, v_t \sim N(0, 1)\) IID.

Here \(\rho\) is the persistence parameter, while \(\sigma_w\) and \(\sigma_v\) are shock standard deviations.

The LinearStateSpace class parameterizes measurement noise by a matrix \(H\) such that \(R = HH^\top\).

# Model parameters
ρ = 0.9
σ_w = 0.5
σ_v = 1.0

# State-space matrices
A = np.array([[ρ]])
C = np.array([[σ_w]])
G = np.array([[1.0]])
R = np.array([[σ_v**2]])

# Build a LinearStateSpace and a Kalman filter object
H = np.array([[σ_v]])   # measurement noise factor: R = H @ H.T
lss = qe.LinearStateSpace(
  A, C, G, H, mu_0=np.zeros(1), Sigma_0=np.eye(1) * 10.0)
kf = qe.Kalman(lss)
kf.set_state(np.zeros(1), np.eye(1) * 10.0)  # diffuse prior

We first simulate a sample path of the true hidden state and noisy observations.

T = 200
x_path, y_path = lss.simulate(ts_length=T, random_state=42)

# Shapes: x_path is (n, T+1), y_path is (m, T)
x_true = x_path[0, :T]
y_obs = y_path[0, :]

We then run the Kalman filter manually, step by step, to collect the filtered estimates.

x_hats = np.zeros(T)
Sigmas = np.zeros(T)

for t in range(T):
    kf.update(y_obs[t:t+1])          # one full filter cycle
    x_hats[t] = kf.x_hat.item()
    Sigmas[t] = kf.Sigma.item()

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)

t_range = np.arange(T)

axes[0].plot(t_range, x_true, lw=2,
        label='true state $x_t$')
axes[0].plot(t_range, x_hats, lw=2,
        linestyle='--', label=r'$\hat{x}_t$ (Kalman)')
axes[0].plot(t_range, y_obs, alpha=0.35, lw=2, label='observation $y_t$')
axes[0].set_title('state and observation')
axes[0].legend(fontsize=9)

axes[1].plot(t_range, Sigmas, color='C1', lw=2,
             label=r'conditional variance $\Sigma_t$')
axes[1].axhline(kf.Sigma_infinity[0, 0], ls='--', color='k',
                label=r'steady-state $\Sigma_\infty$')
axes[1].set_title('conditional variance')
axes[1].legend(fontsize=9)

axes[2].plot(t_range, y_obs - x_hats, color='C2', lw=2, alpha=0.7,
             label=r'innovation $a_t = y_t - G\hat{x}_t$')
axes[2].set_title('innovation')
axes[2].set_xlabel('time $t$')
axes[2].legend(fontsize=9)
fig.tight_layout()
plt.show()
_images/3758a4e406d65c7d05d8682f4b74a7ef9468cf081b0b15cbfbf82aad99413e84.png

Fig. 47.1 Scalar Kalman filtering#

After a short adjustment, the Kalman estimate follows the hidden state much more closely than the raw noisy observations.

The conditional variance drops quickly from the diffuse prior and then settles onto its steady-state value.

The innovation series fluctuates around zero, as a one-step forecast error should.

47.9.2. Convergence of the Riccati equation#

The Kalman class computes the steady-state covariance \(\Sigma_\infty\) by solving the discrete algebraic Riccati equation directly.

Sigma_inf, K_inf = kf.stationary_values()

print(f"Steady-state covariance  Σ_inf = {Sigma_inf[0, 0]:.6f}")
print(f"Kalman filter converged to Σ_t = {Sigmas[-1]:.6f}")
print(f"Steady-state Kalman gain K  = {K_inf[0, 0]:.6f}")

A_minus_KG = A - K_inf @ G
eigval = np.linalg.eigvals(A_minus_KG)[0]
print(f"\nEigenvalue of (A - KG)      = {eigval:.6f}")
print(f"Stable VAR: {np.abs(eigval) < 1}")

Steady-state covariance  Σ_inf = 0.530899
Kalman filter converged to Σ_t = 0.530899
Steady-state Kalman gain K  = 0.312110

Eigenvalue of (A - KG)      = 0.587890
Stable VAR: True

47.9.3. The VAR representation#

Using (47.20), the coefficients in the infinite-order VAR representation are \(G(A - KG)^j K\) for \(j = 0, 1, 2, \ldots\)

We retrieve them via stationary_coefficients:

J = 30
var_coeffs = kf.stationary_coefficients(J, coeff_type='var')

# Lag j+1 coefficient matrices
lags = np.arange(1, J + 1)
coeff_values = np.array([var_coeffs[j][0, 0] for j in range(J)])

fig, ax = plt.subplots()
ax.stem(lags, coeff_values, basefmt=' ')
ax.set_xlabel('lag $j$')
ax.set_ylabel(r'VAR coefficient $G(A{-}KG)^{j-1}K$')
fig.tight_layout()
plt.show()
_images/5882993d2a1f139f2620ebffe76a68828b09fd1d84ae9a036ffced492c763bd5.png

Fig. 47.2 VAR coefficients from the innovations representation#

Most of the autoregressive weight is concentrated in the first few lags.

The later coefficients are nearly zero, so a short finite-lag VAR captures most of the infinite-order representation in this example.

47.9.4. Likelihood evaluation#

We use (47.17) to evaluate the log-likelihood of the simulated sample.

def log_likelihood(A, C, G, R, y_data, x_hat_0, Sigma_0):
    """Evaluate the log-likelihood using the Kalman filter recursions."""
    H_ = np.linalg.cholesky(R)   # R = H_ @ H_.T
    lss_ = qe.LinearStateSpace(A, C, G, H_, mu_0=x_hat_0, Sigma_0=Sigma_0)
    kf_ = qe.Kalman(lss_)
    kf_.set_state(x_hat_0, Sigma_0)

    T_, m_ = y_data.shape
    loglik = 0.0

    for t in range(T_):
        x_h = kf_.x_hat
        Sig = kf_.Sigma
        Omega = G @ Sig @ G.T + R        # innovation covariance
        a_t = y_data[t] - (G @ x_h).flatten()

        sign, logdet = np.linalg.slogdet(Omega)
        loglik += -0.5 * (m_ * np.log(2 * np.pi) + logdet
                          + float(a_t @ np.linalg.solve(Omega, a_t)))
        kf_.update(y_data[t])

    return loglik


y_data_col = y_obs.reshape(-1, 1)
ll = log_likelihood(A, C, G, R,
                    y_data_col,
                    np.zeros(1), np.eye(1) * 10.0)
print(f"Log-likelihood of sample: {ll:.4f}")

Log-likelihood of sample: -325.2335

47.10. An example#

We now work through a structured example that shows how a bivariate VAR(2) fits naturally into the state space framework and how the Kalman filter delivers a Wold (innovations) representation.

47.10.1. A linear state-space system and its filter#

The state and observation equations are

(47.27)#\[ x_{t+1} = A x_t + C w_{t+1} \]
(47.28)#\[ y_t = G x_t + v_t \]

with initial condition and shock distributions

\[ x_0 \sim N(\hat{x}_0, \Sigma_0), \quad w_{t+1} \sim N(0, I), \quad v_t \sim N(0, R). \]

The steady-state error covariance matrix \(\Sigma\) satisfies the Riccati equation

(47.29)#\[ \Sigma = A \Sigma A^\top + CC^\top - A \Sigma G^\top \bigl(G \Sigma G^\top + R\bigr)^{-1} G \Sigma A^\top \]

and the associated steady-state Kalman gain is

(47.30)#\[ K = A \Sigma G^\top \bigl(G \Sigma G^\top + R\bigr)^{-1} \]

Starting from an initial estimate \(\hat{x}_0\), the Kalman filter updates the state estimate via

(47.31)#\[ \hat{x}_{t+1} = A \hat{x}_t + K a_t \]

where the innovation is

(47.32)#\[ a_t = y_t - G \hat{x}_t \]

Substituting (47.32) into (47.31) and expanding:

(47.33)#\[ \hat{x}_{t+1} = A \hat{x}_t + K(y_t - G\hat{x}_t) = (A - KG)\hat{x}_t + K y_t = (A - KG)\hat{x}_t + K G x_t + K v_t \]

47.10.2. Impulse responses of \(y_t\) to the innovations \(a_t\)#

It is useful to compute the ordinary impulse response functions of the observable vector \(y_t\) to its own innovations \(a_t\), the moving-average (Wold) representation that is the mirror image of the VAR (47.20).

From the time-invariant innovations representation (47.18)

\[ \hat{x}_{t+1} = A\hat{x}_t + K a_t, \qquad y_t = G\hat{x}_t + a_t, \]

the moving-average representation (47.25) is

\[ y_t = \bigl[I + G(I - AL)^{-1} K L\bigr]\, a_t = a_t + \sum_{h=1}^{\infty} G A^{h-1} K\, a_{t-h}. \]

Hence the impulse response of \(y_t\) to a unit innovation \(a_t\) is

(47.34)#\[ \Psi_0 = I, \qquad \Psi_h = G A^{h-1} K \quad (h \ge 1). \]

These coefficients decay at the rate governed by the eigenvalues of \(A\).

We can read the coefficients (47.34) directly off a quantecon LinearStateSpace object.

We build a state-space system whose state is the filtered estimate \(\hat{x}_t\), whose single “shock” is the innovation \(a_t\) loaded through \(C = K\), and whose observation matrix is \(G\).

The impulse_response method of that object returns the sequence \(G A^{j} K\) for \(j = 0, 1, 2, \ldots\), which are exactly the \(\Psi_h\) for \(h \ge 1\); we prepend \(\Psi_0 = I\) to capture the contemporaneous feed-through \(y_t = G\hat{x}_t + a_t\).

The array returned below has entry [h, i, j] equal to the response of observable i at horizon h to innovation component j.

def y_to_a_irf(A, K, G, T=40):
    """
    Return Wold IRFs of y_t to its own innovations a_t.
    """
    n, m = A.shape[0], G.shape[0]
    lss = qe.LinearStateSpace(A, K, G, np.zeros((m, m)), mu_0=np.zeros(n))
    _, ycoef = lss.impulse_response(j=T - 2)      # [GK, GAK, GA^2K, ...]
    Psi = np.empty((T, m, m))
    Psi[0] = np.eye(m)                          # contemporaneous response
    for h in range(1, T):
        Psi[h] = ycoef[h - 1]
    return Psi

47.10.3. Bivariate VAR(2) in state-space form#

Consider two observable series \(r_t\) and \(z_t\).

Stack them into the state vector \(x_t = (r_t,\; r_{t-1},\; z_t,\; z_{t-1})^\top\).

We posit the VAR(2) state-transition equation:

(47.35)#\[\begin{split} \begin{pmatrix} r_{t+1} \\ r_t \\ z_{t+1} \\ z_t \end{pmatrix} = \begin{pmatrix} d_1 & d_2 & d_3 & d_4 \\ 1 & 0 & 0 & 0 \\ \delta_1 & \delta_2 & \delta_3 & \delta_4 \\ 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} r_t \\ r_{t-1} \\ z_t \\ z_{t-1} \end{pmatrix} + \begin{pmatrix} c_{11} & c_{12} \\ 0 & 0 \\ c_{21} & c_{22} \\ 0 & 0 \end{pmatrix} \begin{pmatrix} w_{1,t+1} \\ w_{2,t+1} \end{pmatrix} \end{split}\]

We consider two possible observation equations.

The first is a bivariate observation of \(r_t\) and \(z_t\):

(47.36)#\[\begin{split} \begin{pmatrix} r_t \\ z_t \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} r_t \\ r_{t-1} \\ z_t \\ z_{t-1} \end{pmatrix} + \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} v_{1t} \\ v_{2t} \end{pmatrix} \end{split}\]

The second is a univariate observation of \(r_t\):

(47.37)#\[\begin{split} y_t = \begin{pmatrix} 1 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} r_t \\ r_{t-1} \\ z_t \\ z_{t-1} \end{pmatrix} + v_{1t} \end{split}\]

We now compare the Wold impulse responses generated by these two observation systems.

System 1 observes both \(r_t\) and \(z_t\), so its innovation \(a_t\) is \(2 \times 1\).

System 2 observes only \(r_t\), so its innovation \(u_t\) is scalar.

The transition matrices are the same in the two systems, but the observation matrix changes.

Consequently, the steady-state Kalman gain changes too, and so do the Wold responses of the observables to their own innovations.

47.10.4. Numerical example: impulse responses to innovations#

The parameter values are:

\[\begin{split} \begin{aligned} d_1 &= 0.80,\quad d_2 = 0.05,\quad d_3 = 0.75,\quad d_4 = -0.72 \\ \delta_1 &= 0.00,\quad \delta_2 = 0.00,\quad \delta_3 = 0.75,\quad \delta_4 = 0.20 \\ c_{11} &= 1.0,\quad c_{12} = 0.0,\quad c_{21} = 0.0,\quad c_{22} = 1.0 \\ R &= 0.0001 \times I_2 \quad \text{(bivariate case)}, \qquad R = 0.0001 \quad \text{(univariate case)}. \end{aligned} \end{split}\]

These give the \(4 \times 4\) transition matrix and \(4 \times 2\) shock-loading matrix

\[\begin{split} A = \begin{pmatrix} 0.80 & 0.05 & 0.75 & -0.72 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0.75 & 0.20 \\ 0 & 0 & 1 & 0 \end{pmatrix}, \qquad C = \begin{pmatrix} 1 & 0 \\ 0 & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix}. \end{split}\]

System 1 uses the bivariate observation equation (47.36), so \(G\) selects \((r_t, z_t)^\top\) from the state and the innovation \(a_t\) is \(2 \times 1\).

System 2 uses the univariate observation equation (47.37), so the row vector in that equation selects only \(r_t\) and the innovation \(u_t\) is scalar.

# Parameters
d1, d2, d3, d4 = 0.80, 0.05, 0.75, -.72
δ1, δ2, δ3, δ4 = 0.00, 0.00, 0.75, 0.20
c11, c12, c21, c22 = 1.0,  0.0,  0.0,  1.0
σ_v = 0.01  # sqrt(0.0001)

# Shared matrices
A_var = np.array([[d1,     d2,     d3,     d4    ],
                  [1.0,    0.0,    0.0,    0.0   ],
                  [δ1, δ2, δ3, δ4],
                  [0.0,    0.0,    1.0,    0.0   ]])

C_var = np.array([[c11, c12],
                  [0.0, 0.0],
                  [c21, c22],
                  [0.0, 0.0]])

# System 1: bivariate observation
G_biv = np.array([[1.0, 0.0, 0.0, 0.0],
                  [0.0, 0.0, 1.0, 0.0]])
H_biv = σ_v * np.eye(2)          # H @ H.T = 0.0001 * I_2

lss_biv = qe.LinearStateSpace(A_var, C_var, G_biv, H_biv,
                               mu_0=np.zeros(4), Sigma_0=np.eye(4))
kf_biv = qe.Kalman(lss_biv)
_, K_biv = kf_biv.stationary_values()

print("System 1 - steady-state Kalman gain K (4x2):")
print(np.round(K_biv, 5))

# System 2: univariate observation
G_uni = np.array([[1.0, 0.0, 0.0, 0.0]])
H_uni = np.array([[σ_v]])         # H @ H.T = 0.0001

lss_uni = qe.LinearStateSpace(A_var, C_var, G_uni, H_uni,
                               mu_0=np.zeros(4), Sigma_0=np.eye(4))
kf_uni = qe.Kalman(lss_uni)
_, K_uni = kf_uni.stationary_values()

print("\nSystem 2 - steady-state Kalman gain K (4x1):")
print(np.round(K_uni, 5))

System 1 - steady-state Kalman gain K (4x2):
[[7.9987e-01 7.4987e-01]
 [9.9990e-01 0.0000e+00]
 [1.0000e-05 7.4994e-01]
 [0.0000e+00 9.9990e-01]]

System 2 - steady-state Kalman gain K (4x1):
[[0.72306]
 [0.99994]
 [0.31829]
 [0.30984]]

We now apply the helper y_to_a_irf defined above to compute the ordinary impulse responses (47.34) of the observable \(y_t\) to its own innovations \(a_t\), for both System 1 (bivariate, so \(a_t\) is \(2 \times 1\)) and System 2 (univariate, so \(u_t\) is scalar).

T_irf = 40
horizons = np.arange(T_irf)

Psi_biv = y_to_a_irf(A_var, K_biv, G_biv, T_irf)   # System 1: (T, 2, 2)
Psi_uni = y_to_a_irf(A_var, K_uni, G_uni, T_irf)   # System 2: (T, 1, 1)

obs_labels = [r'$r_t$', r'$z_t$']
innov_labels = [r'$a_{1,t}$', r'$a_{2,t}$']

# System 1 responses
fig, axes = plt.subplots(2, 2, figsize=(10, 6), sharex=True)
for i, obs in enumerate(obs_labels):
    for j, inn in enumerate(innov_labels):
        ax = axes[i, j]
        ax.plot(horizons, Psi_biv[:, i, j], lw=2)
        ax.axhline(0, color='k', lw=0.6, ls='--')
        ax.set_title(fr'{obs} to {inn}', fontsize=9)
        if i == 1:
            ax.set_xlabel('horizon $h$')
        if j == 0:
            ax.set_ylabel('response')
fig.tight_layout()
plt.show()
_images/e4ce0c2a19cd7305229c11be4bcfaba750677e416edd70b72a06a405f3bb22c4.png

Fig. 47.3 System 1 responses to own innovations#

Own innovations move their own observables one-for-one on impact and then fade.

The diagonal panels start at one and the off-diagonal panels start at zero because \(\Psi_0 = I\).

The cross response of \(r_t\) to \(a_{2,t}\) is sizeable at short horizons, while the response of \(z_t\) to \(a_{1,t}\) is tiny on the displayed scale.

fig, ax = plt.subplots()
ax.plot(horizons, Psi_uni[:, 0, 0], lw=2)
ax.axhline(0, color='k', lw=0.6, ls='--')
ax.set_xlabel('horizon $h$')
ax.set_ylabel('response')
fig.tight_layout()
plt.show()
_images/a6963517fd7d2cecae78f5cee86ecacdcc3c87524daca96467289823fbd72688.png

Fig. 47.4 System 2 response to own innovation#

With only \(r_t\) observed, the single innovation moves \(r_t\) one-for-one on impact and then decays monotonically toward zero.

For \(h \ge 1\) the responses propagate through the state matrix \(A\) and decay geometrically, tracing out the Wold moving-average representation of the bivariate (System 1) and univariate (System 2) processes.

These are forecast-error responses from the Wold representation, not structural shock responses.

With the example in place, we can now pull together the main lessons of the lecture.

47.11. Summary#

The Kalman filter recursively updates beliefs about a hidden state by combining a prior forecast with the new information contained in the current observation.

The Riccati equation tracks how the filter’s conditional covariance evolves, and its steady state makes the filter time invariant.

At that steady state, the innovations representation writes the observable process in terms of one-step forecast errors that are white noise by construction.

Solving the innovations representation backward gives an infinite-order VAR, while solving it forward gives the Wold moving-average representation.

The numerical examples show that changing the observed variables changes the Kalman gain and therefore changes the Wold responses, even when the underlying state dynamics are the same.

47.12. Exercises#

Exercise 47.1

Consider the scalar AR(1) state space system used above with \(\rho = 0.9\), \(\sigma_w = 0.5\), \(\sigma_v = 1.0\).

Derive an algebraic expression for the steady-state conditional variance \(\Sigma_\infty\) by solving the scalar Riccati equation (47.14) at its fixed point \(\Sigma_{t+1} = \Sigma_t = \Sigma\).

Show that \(\Sigma\) satisfies a quadratic equation, find its positive root, and verify numerically that your formula matches kf.Sigma_infinity.

Exercise 47.2

This exercise considers a two-dimensional state with a one-dimensional observation:

\[\begin{split} A = \begin{pmatrix} 0.9 & 0.1 \\ 0 & 0.8 \end{pmatrix}, \quad C = \begin{pmatrix} 0.4 \\ 0.1 \end{pmatrix}, \quad G = \begin{pmatrix} 1 & 0 \end{pmatrix}, \quad R = [0.5] \end{split}\]

  1. Simulate \(T = 500\) observations from this system starting from a diffuse prior.

  2. Run the Kalman filter and plot both components of \(\hat{x}_t\) against the true hidden state path.

  3. Compute and report the steady-state covariance \(\Sigma_\infty\) and Kalman gain \(K_\infty\).

  4. Check that the eigenvalues of \(A - K_\infty G\) lie strictly inside the unit circle, confirming that the VAR representation (47.20) is stable.

Exercise 47.3

This exercise studies likelihood and parameter estimation using the scalar model from the main text with true parameters \((\rho, \sigma_w, \sigma_v) = (0.9, 0.5, 1.0)\):

  1. Simulate \(T = 300\) observations.

  2. Write a function that evaluates the log-likelihood as a function of \(\rho \in (0, 1)\), holding \(\sigma_w = 0.5\) and \(\sigma_v = 1.0\) fixed, and plot the log-likelihood against \(\rho\) for a grid of values.

  3. Locate the maximum numerically and check that it is close to the true value \(\rho = 0.9\).