6.2 Maximum Likelihood with the Kalman Filter

The basic idea here is that if we can formulate a time series model as a state space model, then we can use the Kalman filter to compute the log-likelihood of the observed data for a given set of parameters. We can then maximize the log-likelihood in the usual way, using the Kalman filter each time to compute the log-likelihood.

Our general state space model formulation described in the previous section has an observation equation \[ y_t = A x_t + V_t \] and a state equation \[ x_t = \Theta x_{t-1} + W_t \] where \(y_t\) is a \(p\times 1\) vector, \(x_t\) is a \(k\times 1\) vector, \(A\) is a \(p\times k\) matrix and \(\Theta\) is \(k\times k\) matrix. We will assume \(V_t\sim\mathcal{N}(0, S)\) and \(W_t\sim\mathcal{N}(0, R)\).

To evaluate and maximize the likelihood function of the data, we need the joint density of the observed data, \(p(y_1,y_2,\dots,y_n)\), which we can subsequently factor into \[\begin{eqnarray*} p(y_1,y_2,\dots,y_n) & = & p(y_1)p(y_2,\dots,y_n\mid y_1)\\ & = & p(y_1)p(y_2\mid y_1)p(y_3,\dots,y_n\mid y_1,y_2)\\ & \vdots & \\ & = & p(y_1)p(y_2\mid y_1)p(y_3\mid y_1, y_2)\cdots p(y_n\mid y_1,\dots,y_{n-1}) \end{eqnarray*}\]

If we pick apart this factorization, we initially need to compute \(p(y_1)\). We can do this by augmenting it with the state variable \(x_1\) and then integrating it out. Although this sounds harder, it actually makes life easier.

\[\begin{eqnarray*} p(y_1) & = & \int p(y_1, x_1)\,dx_1\\ & = & \int p(y_1\mid x_1)p(x_1)\,dx_1 \end{eqnarray*}\]

If you recall from the previous section, \(p(y_1\mid x_1)\) is the density for the observation equation, which in this case is \(\mathcal{N}(Ax_1, S)\). The density \(p(x_1)\) is (again, from the previous section) \(\mathcal{N}(x_1^0, P_1^0)\), where

\[\begin{eqnarray*} x_1^0 & = & \Theta x_0^0\\ P_1^0 & = & \Theta P_0^0\Theta^\prime + R. \end{eqnarray*}\]

Integrating those two densities gives us \[ p(y_1) = \mathcal{N}(Ax_1^0, AP_1^0A^\prime + S). \] If we assume that \(A\) is known (as previously), then the quantities \(x_1^0\) and \(P_1^0\) are all routinely computed in the implementation of the Kalman filtering algorithm.

The general idea here is that we will start with \(t=1\), run the Kalman filtering algorithm for each \(t\) and compute the necessary quantities for each step of the likelihood function. By the time we reach \(t=n\), we should have everything we need to compute the joint likelihood function.

We can see for \(t=2\), we will need to compute \(p(y_2\mid y_1)\), which can be written

\[\begin{eqnarray*} p(y_2\mid y_1) & = & \int p(y_2, x_2\mid y_1)\,dx_2\\ & = & \int p(y_2\mid x_2)p(x_2\mid y_1)\, dx_2\\ & = & \int \varphi(y_2\mid Ax_2, S)\varphi(x_2\mid x_2^1, P_2^1)\,dx_2\\ & = & \mathcal{N}(Ax_2^1, AP_2^1A^\prime + S). \end{eqnarray*}\]

In general, we will have

\[ p(y_t\mid y_1,\dots,y_{t-1}) = \mathcal{N}(Ax_t^{t-1}, A P_t^{t-1}A^\prime + S). \] If we let \(\boldsymbol{\beta}\) represent the vector of unknown parameters, then computing the log-likelihood would require computing the following sum, \[ \ell(\boldsymbol{\beta}) = \sum_{t=1}^n \log p(y_t\mid y_1,\dots,y_{t-1}) \] where we would define \(p(y_1\mid y_0) = p(y_1)\) because there is no \(y_0\). We could then maximize \(\ell(\boldsymbol{\beta})\) with respect to \(\boldsymbol{\beta}\) using standard non-linear maximization routines like Newton’s method or quasi-Newton approaches.