Scalar Kalman filtering recursions

The Kalman filtering recursions, written in scalar form (to gain computational speed and insights) are the following (for generality we let also the variance of the measurement error vary over time).

The initial innovation, its variance and the Kalman gains are: $$\begin{aligned} i_1 &= y_1 - a^{(1)}_1 \\ f_1 &= p^{(11)}_1 + \sigma^2_{\varepsilon,1} \\ k^{(1)}_1 &= \left(p^{(11)}_1 + p^{(12)}_1\right) / f_1 \\ k^{(2)}_1 &= p^{(12)}_1 / f_1 \end{aligned}$$ For $t=1, 2, \ldots, n-1$ the recursions are $$\begin{aligned} a^{(1)}_{t+1} &= a^{(1)}_{t} + a^{(2)}_{t} + k^{(1)}_t i_t \\ a^{(2)}_{t+1} &= a^{(2)}_{t} + k^{(2)}_t i_t \\ \\ p^{(11)}_{t+1} &= p^{(11)}_t + 2 p^{(12)}_t + p^{(22)}_t + \sigma^2_{\eta,t} - k^{(1)}_t k^{(1)}_t f_t \\ p^{(12)}_{t+1} &= p^{(12)}_t + p^{(22)}_t - k^{(1)}_t k^{(2)}_t f_t\\ p^{(22)}_{t+1} &= p^{(22)}_t + \sigma^2_{\zeta,t} - k^{(2)}_t k^{(2)}_t f_t\\ \\ i_{t+1} &= y_{t+1} - a^{(1)}_{t+1} \\ f_{t+1} &= p^{11}_{t+1} + \sigma^2_{\varepsilon,t+1} \\ k^{(1)}_{t+1} &= \left(p^{(11)}_{t+1} + p^{(12)}_{t+1}\right) / f_{t+1} \\ k^{(2)}_{t+1} &= p^{(12)}_{t+1} / f_{t+1} \end{aligned}$$

Modifications when missing observations are present

When one or more values of $y_t$ are missing, then the only modifications to the above recursions are the following. $$\begin{aligned} i_{t+1} &= 0 \\ f_{t+1} &= \infty \\ k^{(1)}_{t+1} &= 0 \\ k^{(1)}_{t+1} &= 0 \end{aligned}.$$

Diffuse initial conditions

Since the two state variables are nonstationary, their initialization should be diffuse: $$\begin{bmatrix} \alpha^{(1)}_1 \\ \alpha^{(2)}_1 \end{bmatrix} \sim N\left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} v & 0\\ 0 & v \end{bmatrix} \right),$$ with $v \rightarrow\infty$.

As it will be clear from the computations below, when $v$ is infinite, the mean squared errors of $a^{(1)}_t$ and $a^{(2)}_t$, and the variances of the innovations are infinite for $t=1, 2$, while from $t=3$ on they are finite.

Let us carry out the computations for $t=1, 2, 3$ and then take the limit for $v \rightarrow\infty$.

Smoothing

The smoothing recursions start from $t=n$ and work backwards down to $t=1$. The following quantities are auxiliar to compute the smoothed values of $\alpha^{(1)}_t$ and their MSE. $$\begin{aligned} r^{(1)}_{n+1} &= 0 \\ r^{(2)}_{n+1} &= 0 \\ n^{(11)}_{n+1} &= 0 \\ n^{(12)}_{n+1} &= 0 \\ n^{(22)}_{n+1} &= 0 \\ e_{n} &= i_{n}/f_{n} \\ d_{n} &= 1/f_{n} \\ \end{aligned}$$ For $t=n, n-1, \ldots, 1$, compute $$\begin{aligned} r^{\left(1\right)}_t &= i_t/f_t + \left(1 - k^{\left(1\right)}_t\right) r^{\left(1\right)}_{t+1} - k^{\left(2\right)}_t r^{\left(2\right)}_{t+1} \\ r^{\left(2\right)}_t &= r^{\left(1\right)}_{t+1} + r^{\left(2\right)}_{t+1} \\ n^{\left(11\right)}_t &= \left(1-k^{\left(1\right)}_t\right)^2 n^{\left(11\right)}_{t+1} - 2\left(1-k^{\left(1\right)}_t\right) k^{\left(2\right)}_t n^{\left(12\right)}_{t+1} + k^{\left(2\right)}_t k^{\left(2\right)}_t n^{\left(22\right)}_{t+1} + 1/f_t \\ n^{\left(12\right)}_t &= \left(1-k^{\left(1\right)}_t\right) \left(n^{\left(11\right)}_{t+1} + n^{\left(12\right)}_{t+1}\right) - k^{\left(2\right)}_t \left(n^{\left(12\right)}_{t+1} + n^{\left(22\right)}_{t+1}\right) \\ n^{\left(22\right)}_t &= n^{\left(11\right)}_{t+1} + 2 n^{\left(12\right)}_{t+1} + n^{\left(22\right)}_{t+1} \\ e_{t-1} &= i_{t-1}/f_{t-1} - k^{(1)}_{t-1} r^{(1)}_{t} - k^{(2)}_{t-1} r^{(2)}_t \\ d_{t-1} &= 1/f_{t-1} + k^{(1)}_{t-1} k^{(1)}_{t-1} n^{(11)}_t + 2 k^{(1)}_{t-1} k^{(2)}_{t-1} n^{(12)}_t + k^{(2)}_{t-1} k^{(2)}_{t-1} n^{(22)}_t \end{aligned}$$ The smoothed values of $\alpha^{(1)}_t$, that is the Hodrick-Prescott filtered time series, and their mean squared errors are given by $$\begin{aligned} a^{(1)}_{t|n} &= a^{(1)}_t + p^{(11)}_t r^{(1)}_t + p^{(12)}_t r^{(2)}_t \\ p^{(11)}_{t|n} &= p^{(11)}_t - p^{(11)}_t p^{(11)}_t n^{(11)}_t - 2 p^{(11)}_t p^{(12)}_t n^{(12)}_t - p^{(12)}_t p^{(12)}_t n^{(22)}_t \end{aligned}$$

Weights for computing the effective degrees of freedom

Since the smoother is linear in the observations, the vector of smoothed $\alpha^{(1)}_t$, say $\mathbf{s}$, is just a linear transformation of the vector of observations, $\mathbf{y}$: $$\mathbf{s} = \mathbf{W} \mathbf{y}.$$ The number of effective degrees of freedom is the trace of the weighting matrix $\mathbf{W}$ (cf. Hastie, Tibshirani and Friedman, 2009, The Elements of Statistical Learning, Section 5.4.1). The formulae for computing such weights in a general state-space form can be found in Koopman and Harvey (2003) Journal of Economic Dynamics and Control, vol. 27. In our framework, the diagonal elements of the matrix $\mathbf{W}$ are given by $$\begin{aligned} w_{tt} &= p^{(11)}_t \left(1/f_t + k^{(1)}_t k^{(1)}_t n^{(11)}_t + 2 k^{(1)}_t k^{(2)}_t n^{(12)}_t + k^{(2)}_t k^{(2)}_t n^{(22)}_t - k^{(1)}_t n^{(11)}_t - k^{(2)}_t n^{(12)}_t\right) \\ &\;\;\;\; - p^{(12)}_t \left(k^{(1)}_t (n^{(11)}_t + n^{(12)}_t) + k^{(2)}_t (n^{(12)}_t + n^{(22)}_t)\right) \end{aligned}$$

Analytical scores

The log-likelihood must be maximised with respect to a very large number of parameters ($n+3$). Thus, providing the numerical optimiser with analytical scores is important for stability and speed. Since all of our parameters are related to quantities in the disturbance covariance matrices, we can adapt the results in Koopman and Shephard (1992, Biometrika vol. 79).

Recall that our (slightly re-parametrised) model is $$\begin{aligned} y_t &= \alpha^{(1)}_t + \varepsilon_t, &\varepsilon_t \sim NID (0, \sigma_\varepsilon^2) \\ \alpha^{(1)}_{t+1} &= \alpha^{(1)}_t + \alpha^{(2)}_t + \eta_t, &\eta_t \sim NID (0, \sigma^2_t) \\ \alpha^{(2)}_{t+1} &= \alpha^{(2)}_t + \zeta_t, &\zeta_t \sim NID (0, \sigma^2 + \gamma^2 \sigma^2_t) \end{aligned}$$ where the parameters to estimate are $\sigma_\varepsilon$, $\sigma$, $\gamma$, and the sequence $\{\sigma_t\}_{t=1,\ldots,n}$, which are all non-negative. Notice that in this parametrisation $\lambda = \sigma_\varepsilon^2 / \sigma^2$.