#+OPTIONS:   H:3 num:nil toc:nil \n:nil @:t ::t |:t ^:nil -:t f:t *:t TeX:t LaTeX:t skip:t d:t tags:not-in-toc creator:t timestamp:nil author:nil title:nil html5-fancy:t
#+HTML_DOCTYPE: html5
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../css/stylesheet-blog.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../css/theorems-en.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../fontawesome/css/all.css" />
#+HTML_HEAD: <link rel="stylesheet" href="../../../highlight/styles/base16/solarized-light.min.css">
#+HTML_HEAD: <script src="../../../highlight/highlight.min.js"></script>
#+HTML_HEAD: <script>hljs.highlightAll();</script>
#+LANGUAGE: en-GB
#+STARTUP: latexpreview
#+TITLE: The Kalman filter and the Hodrick-Prescott filter
#+DATE: September 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-kalman/bin/python

#+BEGIN_QUOTE
The Kalman filter recursively computes the distribution of an unobserved
state vector given past observations, in a linear Gaussian model written in
state space form. This note derives in detail the filter, the smoother,
which also exploits future observations, and the likelihood, which comes out
as a by-product of the filter. It studies the convergence of the filter
towards a steady state and, for the local level model, its link with
exponential smoothing. The application shows that the Hodrick-Prescott
filter is the smoother of a particular trend model, which makes it possible
to estimate by maximum likelihood the smoothing parameter that is usually
fixed a priori. The computations use French GDP per capita since 1820 and
are done in Python, with an implementation of the filter in about thirty
lines. Solutions to the exercises are given in collapsible blocks.
#+END_QUOTE

\\
\\
\\

#+BEGIN_SRC bash :results silent :exports none :async t
  python3 -m venv /tmp/blog-kalman
  source /tmp/blog-kalman/bin/activate
  pip install numpy scipy matplotlib pandas openpyxl statsmodels
#+END_SRC

#+begin_src python :session kalman-en :exports none :results none
  import os
  import urllib.request as url
  import numpy as np
  import pandas as pd
  import matplotlib.pyplot as plt
  from scipy import sparse
  from scipy.sparse.linalg import spsolve
  from scipy.linalg import solve_discrete_lyapunov, solve_discrete_are
  from scipy.optimize import minimize, minimize_scalar
  from scipy.stats import chi2
  rng = np.random.default_rng(2014)
#+end_src

* State space models
:PROPERTIES:
:CUSTOM_ID: modele-etat-mesure
:END:

A state space model separates what is observed, a series $y_t$, from the
unobserved dynamics that generate it, summarised by a state vector
$\alpha_t$. The framework is very general: ARMA processes, unobserved
components models (trend, cycle, seasonality), regressions with
time-varying coefficients and linearised DSGE models can all be written in
this form. We follow the notation of Harvey (/Forecasting, Structural Time
Series Models and the Kalman Filter/, 1989) and of Durbin and Koopman (/Time
Series Analysis by State Space Methods/, 2012), up to a difference in the
timing of the state equation.

#+BEGIN_definition
A state space model is defined by a measurement equation and a state
equation:
\begin{equation*}
\begin{aligned}
y_t &= Z_t\alpha_t + d_t + \varepsilon_t, &\qquad \varepsilon_t &\sim \mathcal N(0, H_t),\\
\alpha_t &= T_t\alpha_{t-1} + c_t + R_t\eta_t, &\qquad \eta_t &\sim \mathcal N(0, Q_t),
\end{aligned}
\end{equation*}
for $t = 1, \dots, n$, with initial condition
\(\alpha_0 \sim \mathcal N(a_0, P_0)\). The observation $y_t$ belongs to
$\mathbb R^p$, the state $\alpha_t$ to $\mathbb R^m$ and the state
disturbance $\eta_t$ to $\mathbb R^r$. The matrices $Z_t$ of size
$p\times m$, $T_t$ of size $m\times m$, $R_t$ of size $m\times r$, $H_t$
of size $p\times p$ and $Q_t$ of size $r\times r$, as well as the vectors
$d_t$ and $c_t$, are non-random. The variables $\alpha_0$,
$\varepsilon_1, \dots, \varepsilon_n$ and $\eta_1, \dots, \eta_n$ are
mutually independent.
#+END_definition

We denote by $Y_t = \{y_1, \dots, y_t\}$ the information available at date
$t$, with $Y_0 = \emptyset$, and by $n$ the sample size, the letter $T$
being reserved for the transition matrix. The state vector summarises the
whole past: the sequence of states is a Markov chain and, given the states,
the observations are independent, $y_t$ depending only on $\alpha_t$. The
selection matrix $R_t$ allows for fewer shocks than states, which is the
rule as soon as the state contains lagged values.

** Examples
:PROPERTIES:
:CUSTOM_ID: exemples
:END:

#+BEGIN_exemple
The local level model is a random walk observed with noise:
\begin{equation*}
y_t = \mu_t + \varepsilon_t,\qquad \mu_t = \mu_{t-1} + \eta_t,
\end{equation*}
with $\varepsilon_t \sim \mathcal N(0, \sigma_\varepsilon^2)$ and
$\eta_t \sim \mathcal N(0, \sigma_\eta^2)$. The state is the level,
$\alpha_t = \mu_t$, and the system matrices are scalars:
\(Z = T = R = 1\), \(H = \sigma_\varepsilon^2\), \(Q = \sigma_\eta^2\) and
$d = c = 0$.
#+END_exemple

#+BEGIN_exemple
The local linear trend model adds a stochastic slope $\beta_t$ to the
level:
\begin{equation*}
y_t = \mu_t + \varepsilon_t,\qquad
\mu_t = \mu_{t-1} + \beta_{t-1} + \eta_t,\qquad
\beta_t = \beta_{t-1} + \zeta_t,
\end{equation*}
where the three disturbances are independent Gaussian white noises, with
variances $\sigma_\varepsilon^2$, $\sigma_\eta^2$ and $\sigma_\zeta^2$. The
state is $\alpha_t = (\mu_t, \beta_t)'$ and
\begin{equation*}
Z = \begin{pmatrix} 1 & 0 \end{pmatrix},\quad
T = \begin{pmatrix} 1 & 1\\ 0 & 1 \end{pmatrix},\quad
R = I_2,\quad
H = \sigma_\varepsilon^2,\quad
Q = \begin{pmatrix} \sigma_\eta^2 & 0\\ 0 & \sigma_\zeta^2 \end{pmatrix}.
\end{equation*}
The special case $\sigma_\eta^2 = 0$ is the /smooth trend model/: the level
is no longer disturbed directly, only the slope is, and the trend is a
smooth curve. Only the second column of $R$ remains, \(R = (0, 1)'\), and
\(Q = \sigma_\zeta^2\). This model plays the leading role in the
application.
#+END_exemple

#+BEGIN_exemple
An autoregressive process of order $p$,
\(y_t = \phi_1 y_{t-1} + \dots + \phi_p y_{t-p} + \eta_t\), is written in
state space form by stacking the lagged values in the state,
\(\alpha_t = (y_t, y_{t-1}, \dots, y_{t-p+1})'\):
\begin{equation*}
Z = \begin{pmatrix} 1 & 0 & \cdots & 0 \end{pmatrix},\quad
T = \begin{pmatrix}
\phi_1 & \phi_2 & \cdots & \phi_{p-1} & \phi_p\\
1 & 0 & \cdots & 0 & 0\\
0 & 1 & \cdots & 0 & 0\\
\vdots & & \ddots & & \vdots\\
0 & 0 & \cdots & 1 & 0
\end{pmatrix},\quad
R = \begin{pmatrix} 1\\ 0\\ \vdots\\ 0 \end{pmatrix},
\end{equation*}
with $H = 0$ and $Q = \sigma^2$. The observation is perfect, without
measurement error: all the uncertainty lies in the state.
#+END_exemple

#+BEGIN_exemple
For an ARMA($p$, $q$) process,
\(y_t = \sum_{i=1}^p \phi_i y_{t-i} + \eta_t + \sum_{j=1}^q \theta_j\eta_{t-j}\),
Harvey (1989) sets $s = \max(p, q+1)$, pads the coefficients with zeros,
$\phi_i = 0$ for $i > p$ and $\theta_j = 0$ for $j > q$, and takes a state
of dimension $s$ with
\begin{equation*}
Z = \begin{pmatrix} 1 & 0 & \cdots & 0 \end{pmatrix},\quad
T = \begin{pmatrix}
\phi_1 & 1 & 0 & \cdots & 0\\
\phi_2 & 0 & 1 & \cdots & 0\\
\vdots & & & \ddots & \\
\phi_{s-1} & 0 & 0 & \cdots & 1\\
\phi_s & 0 & 0 & \cdots & 0
\end{pmatrix},\quad
R = \begin{pmatrix} 1\\ \theta_1\\ \vdots\\ \theta_{s-1} \end{pmatrix},
\end{equation*}
$H = 0$ and $Q = \sigma^2$. For an ARMA(1,1), the state
\(\alpha_t = (\alpha_{1,t}, \alpha_{2,t})'\) satisfies
\(\alpha_{2,t} = \theta\eta_t\) and
\begin{equation*}
\alpha_{1,t} = \phi\alpha_{1,t-1} + \alpha_{2,t-1} + \eta_t
= \phi\alpha_{1,t-1} + \eta_t + \theta\eta_{t-1},
\end{equation*}
and since $y_t = \alpha_{1,t}$, we indeed recover the ARMA(1,1) process.
The case $p = 0$ gives the representation of an MA($q$) process.
#+END_exemple

#+BEGIN_remarque
The state space representation of a process is not unique. For any
invertible $m\times m$ matrix $M$, the state \(\alpha_t^* = M\alpha_t\)
defines a model with the same observations, whose matrices are
\begin{equation*}
\begin{gathered}
Z^* = ZM^{-1},\quad T^* = MTM^{-1},\quad R^* = MR,\\
c^* = Mc,\quad a_0^* = Ma_0,\quad P_0^* = MP_0M'.
\end{gathered}
\end{equation*}
The filter gives \(a_{t|t}^* = Ma_{t|t}\) and \(P_{t|t}^* = MP_{t|t}M'\),
but the innovations, their variances, the likelihood and the forecasts of
$y_t$ are unchanged. For the MA(1) process
\(y_t = \eta_t + \theta\eta_{t-1}\), Harvey's form above, with state
$(y_t, \theta\eta_t)'$, and that of Hamilton (/Time Series Analysis/, 1994),
with state $(\eta_t, \eta_{t-1})'$,
\(Z = \begin{pmatrix} 1 & \theta \end{pmatrix}\),
\(T = \begin{pmatrix} 0 & 0\\ 1 & 0 \end{pmatrix}\) and $R = (1, 0)'$, are
related by
\begin{equation*}
M = \begin{pmatrix} 1 & \theta\\ \theta & 0 \end{pmatrix},
\end{equation*}
which is invertible as soon as $\theta \neq 0$. The dimension of the state
can also be increased by components that the shocks do not reach or that
the observations do not see. A representation is minimal when it is
controllable and observable, and two minimal representations of the same
relation between shocks and observations are related by a change of basis
(Hannan and Deistler, /The Statistical Theory of Linear Systems/, 1988, chapter 2).
Finally, different shock structures can generate the same distribution of
the observations: the local level model, with two independent shocks, has
the distribution of an ARIMA(0,1,1) process with a single shock, as shown
in the section on the [[#forme-reduite][reduced form]]. The choice of
representation is therefore a matter of convenience for filtering and the
likelihood. On the other hand, the parameters of a representation whose
matrices were all free would not be identified: structural models remove
this indeterminacy by fixing $Z$, $T$ and $R$.
#+END_remarque

A stochastic seasonal component, a cycle or explanatory variables are added
in the same way, by augmenting the state vector: this is what makes the
framework useful in practice. In what follows, the system matrices are
constant and there are no deterministic terms ($d = c = 0$). The filter
code handles vector observations, but the likelihood functions and the
applications assume a scalar observation ($p = 1$). A model is represented
by a function that returns its five matrices:

#+begin_src python :session kalman-en :exports code :results none
  def local_level(s2_eps, s2_eta):
      Z, T, R = np.eye(1), np.eye(1), np.eye(1)
      return Z, T, R, np.array([[s2_eps]]), np.array([[s2_eta]])

  def local_linear_trend(s2_eps, s2_eta, s2_zeta):
      Z = np.array([[1.0, 0.0]])
      T = np.array([[1.0, 1.0],
                    [0.0, 1.0]])
      R = np.eye(2)
      return Z, T, R, np.array([[s2_eps]]), np.diag([s2_eta, s2_zeta])

  def smooth_trend(s2_eps, s2_zeta):
      Z, T, _, H, _ = local_linear_trend(s2_eps, 0.0, s2_zeta)
      R = np.array([[0.0],
                    [1.0]])
      return Z, T, R, H, np.array([[s2_zeta]])
#+end_src

* The Kalman filter
:PROPERTIES:
:CUSTOM_ID: filtre
:END:

** Objective
:PROPERTIES:
:CUSTOM_ID: objectif
:END:

We seek the distribution of the state $\alpha_t$ given information $Y_s$.
We denote by
\begin{equation*}
a_{t|s} = \mathbb E\left[\alpha_t \mid Y_s\right],\qquad
P_{t|s} = \mathbb V\left[\alpha_t \mid Y_s\right]
\end{equation*}
its conditional mean and variance matrix. Three cases arise: /prediction/
if $s = t-1$, /filtering/ if $s = t$ and /smoothing/ if $s = n$. The Kalman
filter computes $a_{t|t-1}$, $P_{t|t-1}$, $a_{t|t}$ and $P_{t|t}$ for
$t = 1, \dots, n$, recursively, incorporating one observation at each step.
In a Gaussian model, the conditional distribution of the state is normal,
hence entirely described by these two moments, and the conditional
expectation is the best predictor of $\alpha_t$ in the mean squared error
sense, among all functions of the observations.

#+BEGIN_remarque
Without the normality assumption, but with disturbances uncorrelated with
each other and over time, the recursions derived below remain valid if the
conditional expectation is replaced by the linear projection: $a_{t|t}$ is
then the best /linear/ predictor of $\alpha_t$ as a function of
$y_1, \dots, y_t$, and $P_{t|t}$ the variance matrix of its error (Harvey,
1989, chapter 3).
#+END_remarque

** A lemma on the normal distribution
:PROPERTIES:
:CUSTOM_ID: lemme-gaussien
:END:

The whole derivation rests on an elementary result, the conditional
distribution of a Gaussian vector.

#+BEGIN_lemma
Let
\begin{equation*}
\begin{pmatrix} x\\ y \end{pmatrix} \sim
\mathcal N\left(\begin{pmatrix} \mu_x\\ \mu_y \end{pmatrix},
\begin{pmatrix} \Sigma_{xx} & \Sigma_{xy}\\ \Sigma_{yx} & \Sigma_{yy} \end{pmatrix}\right)
\end{equation*}
be a Gaussian vector, where $\Sigma_{yy}$ is invertible. The distribution
of $x$ given $y$ is normal:
\begin{equation*}
x \mid y \sim \mathcal N\left(\mu_x + \Sigma_{xy}\Sigma_{yy}^{-1}(y - \mu_y),\;
\Sigma_{xx} - \Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}\right).
\end{equation*}
#+END_lemma

#+BEGIN_proof
Let \(e = x - \mu_x - \Sigma_{xy}\Sigma_{yy}^{-1}(y - \mu_y)\). The vector
$(e, y)$ is an affine transformation of $(x, y)$, hence Gaussian. The
covariance between $e$ and $y$ is zero:
\begin{equation*}
\mathrm{Cov}(e, y) = \Sigma_{xy} - \Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yy} = 0.
\end{equation*}
For a Gaussian vector, zero correlation is equivalent to independence: $e$
is independent of $y$. In the decomposition
\begin{equation*}
x = \mu_x + \Sigma_{xy}\Sigma_{yy}^{-1}(y - \mu_y) + e,
\end{equation*}
the first term is a function of $y$ and the second is independent of $y$.
Given $y$, the vector $x$ is therefore the sum of a constant and of a
Gaussian vector with distribution $\mathcal N(0, \mathbb V[e])$, which
establishes the normality of the conditional distribution and gives its
mean. Its variance is
\begin{equation*}
\mathbb V[e] = \mathrm{Cov}(e, x) = \Sigma_{xx} - \Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx},
\end{equation*}
the first equality following from the fact that $e$ differs from $x$ by an
affine function of $y$, with which $e$ is uncorrelated.
#+END_proof

** The filter equations
:PROPERTIES:
:CUSTOM_ID: equations-filtre
:END:

#+BEGIN_proposition
Starting from $a_{0|0} = a_0$ and $P_{0|0} = P_0$, the conditional moments
of the state are computed for $t = 1, \dots, n$ by a /prediction/ step,
\begin{equation*}
\begin{aligned}
a_{t|t-1} &= T_t a_{t-1|t-1} + c_t,\\
P_{t|t-1} &= T_t P_{t-1|t-1} T_t' + R_t Q_t R_t',
\end{aligned}
\end{equation*}
followed by an /updating/ step,
\begin{equation*}
\begin{aligned}
v_t &= y_t - Z_t a_{t|t-1} - d_t, & F_t &= Z_t P_{t|t-1} Z_t' + H_t,\\
K_t &= P_{t|t-1} Z_t' F_t^{-1}, & &\\
a_{t|t} &= a_{t|t-1} + K_t v_t, & P_{t|t} &= (I_m - K_t Z_t) P_{t|t-1}.
\end{aligned}
\end{equation*}
The distribution of $\alpha_t$ given $Y_t$ is normal, with mean $a_{t|t}$
and variance $P_{t|t}$. The vector $v_t$ is the /innovation/, $F_t$ its
variance matrix and $K_t$ the /gain/ of the filter.
#+END_proposition

#+BEGIN_proof
We proceed by induction. Suppose that the distribution of $\alpha_{t-1}$
given $Y_{t-1}$ is $\mathcal N(a_{t-1|t-1}, P_{t-1|t-1})$, which holds for
$t = 1$ since $Y_0$ is empty.

/Prediction./ The observations $y_1, \dots, y_{t-1}$ are functions of
$\alpha_0$ and of the disturbances $\varepsilon_s$ and $\eta_s$ for
$s \leq t-1$: like $\alpha_{t-1}$, they are independent of $\eta_t$. Given
$Y_{t-1}$, the state \(\alpha_t = T_t\alpha_{t-1} + c_t + R_t\eta_t\) is
therefore an affine transformation of two independent Gaussian vectors. Its
distribution is normal, with mean $T_t a_{t-1|t-1} + c_t$ and variance
$T_t P_{t-1|t-1} T_t' + R_t Q_t R_t'$.

/Joint distribution of the state and the observation./ The disturbance
$\varepsilon_t$ is independent of $\alpha_t$ and $Y_{t-1}$. Given
$Y_{t-1}$, the pair formed by $\alpha_t$ and
\(y_t = Z_t\alpha_t + d_t + \varepsilon_t\) is therefore Gaussian:
\begin{equation*}
\begin{pmatrix} \alpha_t\\ y_t \end{pmatrix} \Bigm| Y_{t-1} \sim
\mathcal N\left(\begin{pmatrix} a_{t|t-1}\\ Z_t a_{t|t-1} + d_t \end{pmatrix},
\begin{pmatrix} P_{t|t-1} & P_{t|t-1} Z_t'\\ Z_t P_{t|t-1} & F_t \end{pmatrix}\right),
\end{equation*}
the conditional covariance coming from
\(\mathrm{Cov}(\alpha_t, Z_t\alpha_t + \varepsilon_t \mid Y_{t-1}) = P_{t|t-1}Z_t'\).

/Updating./ Conditioning on $Y_t$ amounts to conditioning on $Y_{t-1}$ and
then on $y_t$. We therefore apply the lemma to the previous joint
distribution, with $x = \alpha_t$ and $y = y_t$. The distribution of
$\alpha_t$ given $Y_t$ is normal, with mean
\begin{equation*}
a_{t|t-1} + P_{t|t-1}Z_t'F_t^{-1}\left(y_t - Z_t a_{t|t-1} - d_t\right) = a_{t|t-1} + K_t v_t
\end{equation*}
and variance
\(P_{t|t-1} - P_{t|t-1}Z_t'F_t^{-1}Z_t P_{t|t-1} = (I_m - K_t Z_t)P_{t|t-1}\),
which completes the induction.
#+END_proof

#+BEGIN_remarque
The gain determines how much the observation corrects the prediction. In
the scalar case, \(K_t = P_{t|t-1}Z_t/(Z_t^2P_{t|t-1} + H_t)\). If the
measurement is perfect ($H_t = 0$), the gain is $1/Z_t$ and the filtered
state follows exactly from the observation. If the measurement is very
noisy ($H_t \to \infty$), or if the state is already known precisely
($P_{t|t-1} \to 0$), the gain tends to zero and the observation is ignored.
Note also that the variances $P_{t|t-1}$, $F_t$ and $P_{t|t}$ do not depend
on the observations: they could be computed before the data are available.
#+END_remarque

#+BEGIN_property
The filtered variance can also be written in /Joseph form/:
\begin{equation*}
P_{t|t} = (I_m - K_t Z_t) P_{t|t-1} (I_m - K_t Z_t)' + K_t H_t K_t'.
\end{equation*}
#+END_property

#+BEGIN_proof
Expanding the right-hand side gives
\begin{equation*}
P_{t|t-1} - K_t Z_t P_{t|t-1} - P_{t|t-1} Z_t' K_t' + K_t (Z_t P_{t|t-1} Z_t' + H_t) K_t'.
\end{equation*}
The last term equals \(K_t F_t K_t' = P_{t|t-1}Z_t'K_t'\), since
\(K_t F_t = P_{t|t-1}Z_t'\). It cancels the third one, leaving
$(I_m - K_tZ_t)P_{t|t-1}$.
#+END_proof

The two expressions of $P_{t|t}$ are equal in exact arithmetic, but not in
floating-point arithmetic. The Joseph form is a sum of two symmetric
positive semi-definite matrices: it remains so whatever the rounding errors
on the gain, whereas the short form subtracts two nearly equal terms when
$P_{t|t-1}$ is large and may produce a non-symmetric matrix, or even
negative variances. We shall see that this happens at initialisation. The
implementation follows the equations to the letter, with the Joseph form. A
missing observation, coded as =np.nan=, is handled by skipping the updating
step: without an observation, the filtered distribution is the predicted
one.

#+begin_src python :session kalman-en :exports code :results none
  def kalman_filter(y, Z, T, R, H, Q, a0, P0):
      y = np.asarray(y, dtype=float).reshape(len(y), -1)
      n, p = y.shape
      m = len(a0)
      RQR = R @ Q @ R.T
      res = {'a_pred': np.zeros((n, m)), 'P_pred': np.zeros((n, m, m)),
             'a_filt': np.zeros((n, m)), 'P_filt': np.zeros((n, m, m)),
             'v': np.full((n, p), np.nan), 'F': np.full((n, p, p), np.nan),
             'K': np.zeros((n, m, p))}
      a, P = np.asarray(a0, dtype=float), np.asarray(P0, dtype=float)
      for t in range(n):
          # Prediction
          a = T @ a
          P = T @ P @ T.T + RQR
          res['a_pred'][t], res['P_pred'][t] = a, P
          # Updating, skipped if the observation is missing
          if not np.isnan(y[t]).any():
              v = y[t] - Z @ a
              F = Z @ P @ Z.T + H
              K = np.linalg.solve(F, Z @ P).T          # P Z' F^{-1}
              a = a + K @ v
              IKZ = np.eye(m) - K @ Z
              P = IKZ @ P @ IKZ.T + K @ H @ K.T       # Joseph form
              res['v'][t], res['F'][t], res['K'][t] = v, F, K
          res['a_filt'][t], res['P_filt'][t] = a, P
      return res
#+end_src

The function returns a dictionary of arrays indexed by date, from $0$ to
$n-1$ in Python for dates $1$ to $n$ in the note.

** The innovations
:PROPERTIES:
:CUSTOM_ID: innovations
:END:

#+BEGIN_property
The innovation $v_t$ is independent of $Y_{t-1}$ and has distribution
$\mathcal N(0, F_t)$. The innovations $v_1, \dots, v_n$ are independent:
they form a Gaussian white noise, heteroskedastic if $F_t$ varies with $t$.
#+END_property

#+BEGIN_proof
By the joint distribution established in the proof of the filter, the
distribution of $y_t$ given $Y_{t-1}$ is
$\mathcal N(Z_t a_{t|t-1} + d_t, F_t)$, hence that of $v_t$ given $Y_{t-1}$
is $\mathcal N(0, F_t)$. This conditional distribution does not depend on
$Y_{t-1}$, since $F_t$ is not random. The innovation $v_t$ is therefore
independent of $Y_{t-1}$, and its distribution is $\mathcal N(0, F_t)$. The
past innovations $v_1, \dots, v_{t-1}$ are functions of $Y_{t-1}$, of which
$v_t$ is independent. Applying this argument for $t = n, n-1, \dots, 2$,
the joint density of $(v_1, \dots, v_n)$ is the product of the marginal
densities.
#+END_proof

The innovation $v_t$ is the part of $y_t$ that was not predictable at date
$t-1$. Since $v_t$ equals $y_t$ minus a function of $Y_{t-1}$, the sequence
of innovations up to date $t$ and that of the observations can be deduced
from each other: they carry the same information. If the model is correctly
specified, the /standardised innovations/ \(e_t = F_t^{-1/2}v_t\) are
independent standard normal. Their sample autocorrelations therefore provide
a simple specification diagnostic, which the application uses extensively.

** Initialisation
:PROPERTIES:
:CUSTOM_ID: initialisation
:END:

It remains to choose $a_0$ and $P_0$. If the model is stationary, with
constant matrices and eigenvalues of $T$ of modulus strictly less than one,
$\alpha_0$ is drawn from the stationary distribution of the state:
\(a_0 = (I_m - T)^{-1}c\) and $P_0$ is the solution of the Lyapunov
equation
\begin{equation*}
P_0 = T P_0 T' + R Q R',\qquad\text{that is}\qquad
\mathrm{vec}(P_0) = (I_{m^2} - T\otimes T)^{-1}\,\mathrm{vec}(RQR'),
\end{equation*}
computed by SciPy's =solve_discrete_lyapunov= function.

If the state has non-stationary components, such as the level of a random
walk, this stationary distribution does not exist. The ignorance about the
initial state is then expressed by a very dispersed distribution, $a_0 = 0$
and $P_0 = \kappa I_m$ with $\kappa$ large: this is the /approximate diffuse
initialisation/. In the local level model, the first step of the filter
gives \(P_{1|0} = \kappa + \sigma_\eta^2\), then
\begin{equation*}
K_1 = \frac{P_{1|0}}{P_{1|0} + \sigma_\varepsilon^2},\qquad
a_{1|1} = K_1 y_1,\qquad
P_{1|1} = \frac{\sigma_\varepsilon^2 P_{1|0}}{P_{1|0} + \sigma_\varepsilon^2}.
\end{equation*}
As $\kappa$ tends to infinity, $a_{1|1}$ tends to $y_1$ and $P_{1|1}$ to
$\sigma_\varepsilon^2$: the first observation pins down the level, with the
precision of the measurement, and the initial distribution is forgotten.
More generally, if $d$ components of the state are diffuse, the first $d$
observations serve to determine them, after which the effect of $\kappa$ is
of order $1/\kappa$. Three precautions are needed.

- The approximation is all the better as $\kappa$ is large relative to the
  variance of the data and to the square of their distance from $a_0$. With
  data of the order of $1\,000$, as in the application, $\kappa$ must be
  well above $10^6$.
- The variance $P_{t|t}$ of the first steps is the difference of two terms
  of order $\kappa$, and the computation loses about $\log_{10}\kappa$ of
  the sixteen significant digits of floating-point arithmetic. The Joseph
  form preserves symmetry and positivity, but an excessive value of
  $\kappa$ degrades precision. We use $\kappa = 10^7$ for the simulations
  and $\kappa = 10^9$ for the application, after checking.
- The variances $F_t$ of the first $d$ innovations are of order $\kappa$:
  these innovations carry no information about the parameters, and their
  contribution to the likelihood depends on the arbitrary value of
  $\kappa$. They must be discarded, as we shall see in the
  [[#concentree][section on estimation]].

#+BEGIN_remarque
The /exact diffuse initialisation/ handles the limit $\kappa \to \infty$
analytically, by separating in $P_{t|t-1}$ a term proportional to $\kappa$
and a finite term (De Jong, /Annals of Statistics/, 1991; Koopman, /Journal
of the American Statistical Association/, 1997; Durbin and Koopman, 2012,
chapter 5). It removes the arbitrariness of the choice of $\kappa$ at the
cost of heavier recursions during the first periods. The =statsmodels=
library implements it; it will serve as an external check in the
application.
#+END_remarque

** An illustration
:PROPERTIES:
:CUSTOM_ID: illustration
:END:

We simulate $n = 100$ observations of a local level model, with
$\sigma_\varepsilon = 1$ and $\sigma_\eta = 0.5$, and apply the filter with
an approximate diffuse initialisation.

#+begin_src python :session kalman-en :exports code :results none
  n = 100
  sigma_eps, sigma_eta = 1.0, 0.5
  mu = np.cumsum(sigma_eta*rng.standard_normal(n))
  y_sim = mu + sigma_eps*rng.standard_normal(n)
  kappa = 1e7
  Z, T, R, H, Q = local_level(sigma_eps**2, sigma_eta**2)
  res_sim = kalman_filter(y_sim, Z, T, R, H, Q, np.zeros(1), kappa*np.eye(1))
  rmse = lambda x: np.sqrt(np.mean((x - mu)**2))
  print(rmse(y_sim), rmse(res_sim['a_filt'][:, 0]))
#+end_src

The filtered state estimates the level with a root mean squared error of
$0.61$, against $1.01$ for the observation itself taken as a naive
estimator. The figure [[fig:niveau-local-filtre][below]] shows that the
filtered state follows the level with a slight lag: using only the past, it
can tell a lasting change in the level from measurement noise only after a
few observations.

#+begin_src python :session kalman-en :exports none :results none
  t_sim = np.arange(1, n+1)
  a_f, s_f = res_sim['a_filt'][:, 0], np.sqrt(res_sim['P_filt'][:, 0, 0])
  fig, ax = plt.subplots(figsize=(6, 4))
  ax.plot(t_sim, y_sim, 'k.', markersize=3, label=r'$y_t$')
  ax.plot(t_sim, mu, 'k', linewidth=0.8, label=r'$\mu_t$')
  ax.plot(t_sim, a_f, 'b', linewidth=1, label=r'$a_{t|t}$')
  ax.fill_between(t_sim, a_f-2*s_f, a_f+2*s_f, color='b', alpha=0.15, linewidth=0)
  ax.set_xlabel(r'$t$')
  ax.legend()
  fig.tight_layout()
  fig.savefig("kalman-niveau-local-filtre.svg", transparent=True)
#+end_src

#+CAPTION: *Simulated local level model ($\sigma_\varepsilon = 1$, $\sigma_\eta = 0.5$): observations, level and filtered state, with a band of plus or minus two conditional standard deviations.*
#+LABEL: fig:niveau-local-filtre
[[file:kalman-niveau-local-filtre.svg]]

* Convergence of the filter
:PROPERTIES:
:CUSTOM_ID: convergence
:END:

When the system matrices are constant, the variance recursion depends
neither on the observations nor on the date. Substituting the updating step
into the prediction step gives a recursion on the prediction variance
alone.

#+BEGIN_corollary
With constant matrices, the prediction variance satisfies the /Riccati
equation/
\begin{equation*}
P_{t+1|t} = T P_{t|t-1} T' - T P_{t|t-1} Z' F_t^{-1} Z P_{t|t-1} T' + RQR',
\qquad F_t = Z P_{t|t-1} Z' + H,
\end{equation*}
and the predicted state satisfies \(a_{t+1|t} = T a_{t|t-1} + c + T K_t v_t\).
#+END_corollary

#+BEGIN_proof
It suffices to substitute $a_{t|t} = a_{t|t-1} + K_t v_t$ and
\(P_{t|t} = P_{t|t-1} - P_{t|t-1}Z'F_t^{-1}ZP_{t|t-1}\) into the prediction
equations of date $t+1$.
#+END_proof

The sequence of prediction variances converges under very general
conditions. We take the following result for granted (Anderson and Moore,
/Optimal Filtering/, 1979, chapter 4, section 4.4).

#+BEGIN_proposition
Suppose that $H$ is positive definite, that the pair $(T, Z)$ is detectable
and that the pair $(T, RQ^{1/2})$ is stabilisable. Then, for any positive
semi-definite initial variance, $P_{t+1|t}$ converges to the unique
positive semi-definite solution $\bar P$ of the algebraic Riccati equation
\begin{equation*}
\bar P = T\bar PT' - T\bar PZ'(Z\bar PZ' + H)^{-1}Z\bar PT' + RQR',
\end{equation*}
and the gain $K_t$ converges to \(\bar K = \bar PZ'\bar F^{-1}\), with
\(\bar F = Z\bar PZ' + H\).
#+END_proposition

Detectability requires that the unstable components of the state show up in
the observations, stabilisability that the shocks reach all these
components. Once the steady state is reached, the filter is a time-invariant
linear filter, \(a_{t+1|t} = T(I_m - \bar KZ)a_{t|t-1} + T\bar Ky_t\), and
computing the variances becomes unnecessary. The matrix $\bar P$ is given by
SciPy's =solve_discrete_are= function, whose convention transposes that of
the filter: =solve_discrete_are(T.T, Z.T, R @ Q @ R.T, H)=.

** The local level case
:PROPERTIES:
:CUSTOM_ID: niveau-local-riccati
:END:

In the local level model, everything can be computed explicitly. Let
\(q = \sigma_\eta^2/\sigma_\varepsilon^2\) denote the signal-to-noise ratio
and \(p_t = P_{t|t-1}/\sigma_\varepsilon^2\) the prediction variance
relative to that of the measurement noise.

#+BEGIN_property
The sequence $p_t$ satisfies $p_{t+1} = f(p_t)$, with
\(f(p) = p/(1+p) + q\). For any $p_1 \geq 0$, it converges monotonically to
the unique positive fixed point
\begin{equation*}
\bar p = \frac{q + \sqrt{q^2 + 4q}}{2},
\end{equation*}
and the gain \(K_t = p_t/(1+p_t)\) converges to
\(\bar K = \bar p/(1 + \bar p)\).
#+END_property

#+BEGIN_proof
With $Z = T = R = 1$, the Riccati equation reads
\(P_{t+1|t} = P_{t|t-1} - P_{t|t-1}^2/(P_{t|t-1} + \sigma_\varepsilon^2) + \sigma_\eta^2\),
that is, dividing by $\sigma_\varepsilon^2$,
\(p_{t+1} = p_t - p_t^2/(1 + p_t) + q = p_t/(1+p_t) + q\). The fixed points
are the roots of the polynomial $p^2 - qp - q$, of which $\bar p$ is the
only positive one. The function $f$ is increasing on $[0, \infty)$, and the
function \(g(p) = f(p) - p = q - p^2/(1+p)\) is decreasing and vanishes at
$\bar p$. If $p_t < \bar p$, we therefore have
\(p_t < f(p_t) < f(\bar p) = \bar p\): the sequence is increasing and
bounded above by $\bar p$. Symmetrically, it is decreasing and bounded below
by $\bar p$ if $p_t > \bar p$. In both cases it converges, and its limit is
a fixed point of $f$, that is $\bar p$.
#+END_proof

Near the fixed point, the distance to $\bar p$ is multiplied at each period
by \(f'(\bar p) = (1 + \bar p)^{-2}\): convergence is geometric, and all the
faster as the signal-to-noise ratio is high. The figure [[fig:gain][below]]
checks this for three values of $q$, with the approximate diffuse
initialisation. The gain starts at one, since the first observation pins
down the level, and reaches its steady-state value, $0.200$, $0.390$ and
$0.618$ for $q$ equal to $0.05$, $0.25$ and $1$, within $10^{-4}$ after
$19$, $10$ and $6$ observations respectively. For $q = 1$, $\bar p$ is the
golden ratio and $\bar K$ its inverse.

#+begin_src python :session kalman-en :exports code :results none
  def local_level_fixed_point(q):
      p_bar = (q + np.sqrt(q**2 + 4*q))/2
      return p_bar, p_bar/(1 + p_bar)

  for q in (0.05, 0.25, 1):
      r = kalman_filter(np.zeros(30), *local_level(1.0, q), np.zeros(1), kappa*np.eye(1))
      K_bar = local_level_fixed_point(q)[1]
      print(q, K_bar, 1 + np.argmax(np.abs(r['K'][:, 0, 0] - K_bar) < 1e-4))
#+end_src

#+begin_src python :session kalman-en :exports none :results none
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  dates = np.arange(1, 31)
  for q, c in ((0.05, 'b'), (0.25, 'g'), (1, 'r')):
      r = kalman_filter(np.zeros(30), *local_level(1.0, q), np.zeros(1), kappa*np.eye(1))
      p_bar, K_bar = local_level_fixed_point(q)
      axes[0].plot(dates, r['P_pred'][:, 0, 0], c+'o-', markersize=3, linewidth=1, label=fr'$q = {q}$')
      axes[0].axhline(p_bar, color=c, linestyle='--', linewidth=0.8)
      axes[1].plot(dates, r['K'][:, 0, 0], c+'o-', markersize=3, linewidth=1, label=fr'$q = {q}$')
      axes[1].axhline(K_bar, color=c, linestyle='--', linewidth=0.8)
  axes[0].set_yscale('log')
  axes[0].set_ylim(0.1, 10)
  axes[0].set_xlabel(r'$t$')
  axes[0].set_ylabel(r'$P_{t|t-1}/\sigma_\varepsilon^2$')
  axes[0].legend()
  axes[1].set_xlabel(r'$t$')
  axes[1].set_ylabel(r'$K_t$')
  axes[1].legend()
  fig.tight_layout()
  fig.savefig("kalman-gain.svg", transparent=True)
#+end_src

#+CAPTION: *Convergence of the prediction variance (left, logarithmic scale) and of the gain (right) in the local level model, for three values of the signal-to-noise ratio. Dotted lines show the steady-state values $\bar p$ and $\bar K$.*
#+LABEL: fig:gain
[[file:kalman-gain.svg]]

** Reduced form and exponential smoothing
:PROPERTIES:
:CUSTOM_ID: forme-reduite
:END:

The steady state of the filter has a remarkable interpretation in the local
level model, through the reduced form of this model.

#+BEGIN_property
In the local level model, the first difference
\(\Delta y_t = \eta_t + \varepsilon_t - \varepsilon_{t-1}\) is an invertible
MA(1) process, \(\Delta y_t = \xi_t + \theta\xi_{t-1}\), where $\xi_t$ is a
white noise with variance \(\sigma_\xi^2 = -\sigma_\varepsilon^2/\theta\)
and
\begin{equation*}
\theta = \frac{\sqrt{q^2 + 4q} - q - 2}{2} \in (-1, 0).
\end{equation*}
The process $y_t$ is therefore an ARIMA(0,1,1).
#+END_property

#+BEGIN_proof
The autocovariances of $\Delta y_t$ are
\(\gamma_0 = \sigma_\eta^2 + 2\sigma_\varepsilon^2\),
\(\gamma_1 = -\sigma_\varepsilon^2\) and $\gamma_k = 0$ for $k \geq 2$:
they are those of an MA(1) process. Those of $\xi_t + \theta\xi_{t-1}$ are
\(\gamma_0 = (1 + \theta^2)\sigma_\xi^2\) and
\(\gamma_1 = \theta\sigma_\xi^2\). Equating the first-order
autocorrelations, \(\theta/(1 + \theta^2) = -1/(q + 2)\), that is
\(\theta^2 + (q + 2)\theta + 1 = 0\). The discriminant, $q^2 + 4q$, is
positive, both roots are negative and their product is one: only one lies
between $-1$ and $0$, the invertible root. Finally,
\(\sigma_\xi^2 = \gamma_1/\theta = -\sigma_\varepsilon^2/\theta\).
#+END_proof

#+BEGIN_property
The steady-state gain is $\bar K = 1 + \theta$, and the steady-state
innovation variance, \(\bar F = \sigma_\varepsilon^2(1 + \bar p)\), equals
$\sigma_\xi^2$. In the steady state, the prediction follows the
/exponential smoothing/ recursion,
\begin{equation*}
a_{t+1|t} = (1 - \bar K)\,a_{t|t-1} + \bar K y_t,
\end{equation*}
and the innovations of the filter are those of the reduced form:
\(\Delta y_t = v_t + \theta v_{t-1}\).
#+END_property

#+BEGIN_proof
Let $r = \sqrt{q^2 + 4q}$, so that $\bar p = (q + r)/2$ and
$1 + \theta = (r - q)/2$. Multiplying numerator and denominator by
$q + 2 + r$, and since $r^2 = q^2 + 4q$,
\begin{equation*}
-\frac{1 + \theta}{\theta} = \frac{r - q}{q + 2 - r}
= \frac{(r - q)(q + 2 + r)}{(q + 2)^2 - r^2}
= \frac{2q + 2r}{4} = \bar p.
\end{equation*}
Hence $1 + \bar p = -1/\theta$, then
\(\bar K = \bar p/(1 + \bar p) = 1 + \theta\) and
\(\bar F = \sigma_\varepsilon^2(1 + \bar p) = -\sigma_\varepsilon^2/\theta = \sigma_\xi^2\).
In the steady state, with $Z = T = 1$, the predicted state is
\(a_{t+1|t} = a_{t|t} = a_{t|t-1} + \bar K(y_t - a_{t|t-1})\), which is the
stated recursion. Finally, $y_t = a_{t|t-1} + v_t$ and
\(a_{t|t-1} - a_{t-1|t-2} = \bar K v_{t-1}\) give
\begin{equation*}
\Delta y_t = \bar K v_{t-1} + v_t - v_{t-1} = v_t - (1 - \bar K)v_{t-1} = v_t + \theta v_{t-1}.
\end{equation*}
#+END_proof

The Kalman filter in steady state is therefore the exponential smoothing
forecasting method, and exponential smoothing is the optimal forecast of an
ARIMA(0,1,1) process (Muth, /Journal of the American Statistical
Association/, 1960). The filter adds two things: the smoothing constant
follows from the parameters of the model, and the gain varies during the
transition phase to account for the initial uncertainty. We check this on
the simulated series, with $q = 0.25$, by comparing the prediction of the
filter with exponential smoothing with constant $\bar K$ initialised at the
first observation:

#+begin_src python :session kalman-en :exports code :results none
  q = sigma_eta**2/sigma_eps**2
  p_bar, K_bar = local_level_fixed_point(q)
  theta = (np.sqrt(q**2 + 4*q) - q - 2)/2
  print(K_bar, 1 + theta)
  exp_smoothing = np.zeros(n)
  exp_smoothing[1] = y_sim[0]
  for t in range(1, n-1):
      exp_smoothing[t+1] = (1 - K_bar)*exp_smoothing[t] + K_bar*y_sim[t]
  gap = np.abs(exp_smoothing[1:] - res_sim['a_pred'][1:, 0])
  print(gap[[3, 8, 18]])
#+end_src

The gap between the two predictions, of the order of $2\times 10^{-2}$ at
date $5$, is only $3\times 10^{-3}$ at date $10$ and $2\times 10^{-5}$ at
date $20$: it decreases geometrically, at rate
$1 - \bar K = -\theta \approx 0.61$ per period. The right panel of the
figure [[fig:lissage-exponentiel][below]] also checks that the standardised
innovations $e_t$ are not autocorrelated, as they should not be when the
model is correctly specified. The sample autocorrelation function will be
used again in the application.

#+begin_src python :session kalman-en :exports code :results none
  def acf(x, kmax):
      x = x - x.mean()
      return np.array([np.sum(x[k:]*x[:-k]) for k in range(1, kmax+1)])/np.sum(x**2)

  e_sim = res_sim['v'][1:, 0]/np.sqrt(res_sim['F'][1:, 0, 0])
  print(acf(e_sim, 10))
#+end_src

#+begin_src python :session kalman-en :exports none :results none
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  axes[0].semilogy(np.arange(2, 41), gap[:39], 'b', linewidth=1)
  axes[0].set_xlabel(r'$t$')
  axes[0].set_ylabel('gap between the two predictions')
  k = np.arange(1, 11)
  band = 1.96/np.sqrt(len(e_sim))
  axes[1].bar(k, acf(e_sim, 10), color='b', width=0.5)
  for s in (-1, 1):
      axes[1].axhline(s*band, color='k', linestyle='--', linewidth=0.8)
  axes[1].axhline(0, color='k', linewidth=0.5)
  axes[1].set_xticks(k)
  axes[1].set_xlabel('lag')
  axes[1].set_ylabel(r'autocorrelation of $e_t$')
  fig.tight_layout()
  fig.savefig("kalman-lissage-exponentiel.svg", transparent=True)
#+end_src

#+CAPTION: *Left, absolute gap between the filter prediction $a_{t|t-1}$ and exponential smoothing with constant $\bar K$ (logarithmic scale). Right, autocorrelations of the standardised innovations, with the band $\pm 1.96/\sqrt{n}$.*
#+LABEL: fig:lissage-exponentiel
[[file:kalman-lissage-exponentiel.svg]]

* The smoother
:PROPERTIES:
:CUSTOM_ID: lisseur
:END:

The filter estimates the state at date $t$ with the information available
at that date. When the whole sample is observed, it is natural to also use
the later observations: this is the purpose of /smoothing/, which computes
$a_{t|n}$ and $P_{t|n}$ for $t = 1, \dots, n$. The filter is one-sided, the
smoother two-sided. The starting point is a Markov property: once the state
at date $t+1$ is known, the later observations bring nothing more about the
state at date $t$.

#+BEGIN_lemma
For $t < n$, the distribution of $\alpha_t$ given $\alpha_{t+1}$ and $Y_n$
is equal to the distribution of $\alpha_t$ given $\alpha_{t+1}$ and $Y_t$.
#+END_lemma

#+BEGIN_proof
Let $Y_{t+1:n} = (y_{t+1}, \dots, y_n)$ denote the observations after date
$t$. By the model equations, they can be written as a function of
$\alpha_{t+1}$ and of the disturbances
\(\varepsilon_{t+1}, \dots, \varepsilon_n, \eta_{t+2}, \dots, \eta_n\),
gathered in a vector $\omega$. The triple formed by $\alpha_t$,
$\alpha_{t+1}$ and $Y_t$ is a function of $\alpha_0$, of
\(\varepsilon_1, \dots, \varepsilon_t\) and of
\(\eta_1, \dots, \eta_{t+1}\), of which $\omega$ is independent. Given
$\alpha_t$, $\alpha_{t+1}$ and $Y_t$, the distribution of $Y_{t+1:n}$ is
therefore that of a function of $\alpha_{t+1}$ and $\omega$: it depends
only on $\alpha_{t+1}$. By Bayes' formula, denoting densities by $p$,
\begin{equation*}
p(\alpha_t \mid \alpha_{t+1}, Y_n)
\propto p(Y_{t+1:n} \mid \alpha_t, \alpha_{t+1}, Y_t)\, p(\alpha_t \mid \alpha_{t+1}, Y_t)
= p(Y_{t+1:n} \mid \alpha_{t+1})\, p(\alpha_t \mid \alpha_{t+1}, Y_t),
\end{equation*}
where the symbol $\propto$ denotes proportionality as functions of
$\alpha_t$. The first factor does not depend on $\alpha_t$: the density of
$\alpha_t$ given $\alpha_{t+1}$ and $Y_n$ is proportional, hence equal, to
that of $\alpha_t$ given $\alpha_{t+1}$ and $Y_t$.
#+END_proof

#+BEGIN_proposition
Starting from $a_{n|n}$ and $P_{n|n}$, computed by the filter, the smoothed
moments are obtained by the backward recursion, for
$t = n-1, n-2, \dots, 1$,
\begin{equation*}
\begin{aligned}
J_t &= P_{t|t} T_{t+1}' P_{t+1|t}^{-1},\\
a_{t|n} &= a_{t|t} + J_t\left(a_{t+1|n} - a_{t+1|t}\right),\\
P_{t|n} &= P_{t|t} + J_t\left(P_{t+1|n} - P_{t+1|t}\right)J_t',
\end{aligned}
\end{equation*}
assuming $P_{t+1|t}$ is invertible. The distribution of $\alpha_t$ given
$Y_n$ is normal, with mean $a_{t|n}$ and variance $P_{t|n}$.
#+END_proposition

This recursion is due to Rauch, Tung and Striebel (/AIAA Journal/, 1965).
It requires only one backward pass over the output of the filter.

#+BEGIN_proof
/Joint distribution./ Given $Y_t$, the state $\alpha_t$ has distribution
$\mathcal N(a_{t|t}, P_{t|t})$ and
\(\alpha_{t+1} = T_{t+1}\alpha_t + c_{t+1} + R_{t+1}\eta_{t+1}\), where
$\eta_{t+1}$ is independent of $\alpha_t$ and $Y_t$. Given $Y_t$, the pair
formed by $\alpha_t$ and $\alpha_{t+1}$ is therefore Gaussian:
\begin{equation*}
\begin{pmatrix} \alpha_t\\ \alpha_{t+1} \end{pmatrix} \Bigm| Y_t \sim
\mathcal N\left(\begin{pmatrix} a_{t|t}\\ a_{t+1|t} \end{pmatrix},
\begin{pmatrix} P_{t|t} & P_{t|t}T_{t+1}'\\ T_{t+1}P_{t|t} & P_{t+1|t} \end{pmatrix}\right).
\end{equation*}

/Conditioning on the next state./ The lemma on the normal distribution
gives
\begin{equation*}
\alpha_t \mid \alpha_{t+1}, Y_t \sim
\mathcal N\left(a_{t|t} + J_t(\alpha_{t+1} - a_{t+1|t}),\; P_{t|t} - J_tP_{t+1|t}J_t'\right),
\end{equation*}
and, by the previous lemma, this is also the distribution of $\alpha_t$
given $\alpha_{t+1}$ and $Y_n$.

/Smoothed mean./ By the law of iterated expectations,
\begin{equation*}
a_{t|n} = \mathbb E\left[\mathbb E[\alpha_t \mid \alpha_{t+1}, Y_n] \mid Y_n\right]
= a_{t|t} + J_t\left(a_{t+1|n} - a_{t+1|t}\right).
\end{equation*}

/Smoothed variance./ By the law of total variance,
\begin{equation*}
\begin{aligned}
P_{t|n} &= \mathbb E\left[\mathbb V[\alpha_t \mid \alpha_{t+1}, Y_n] \mid Y_n\right]
+ \mathbb V\left[\mathbb E[\alpha_t \mid \alpha_{t+1}, Y_n] \mid Y_n\right]\\
&= P_{t|t} - J_tP_{t+1|t}J_t' + J_tP_{t+1|n}J_t'.
\end{aligned}
\end{equation*}
Finally, the distribution of $\alpha_t$ given $Y_n$ is normal because
$\alpha_t$ and $Y_n$ form a Gaussian vector.
#+END_proof

#+BEGIN_remarque
The same law of total variance, applied by conditioning on $Y_t$ rather
than on $\alpha_{t+1}$, gives
\(P_{t|t} = P_{t|n} + \mathbb V\left[a_{t|n} \mid Y_t\right]\): the smoothed
variance is smaller than the filtered variance, in the sense of positive
semi-definite matrices. When $P_{t+1|t}$ is singular, which happens when the
state contains components observed without error, the inverse can be
replaced by a generalised inverse, or one can use the recursions of De Jong
(/Journal of the American Statistical Association/, 1989), which avoid any
inversion.
#+END_remarque

#+begin_src python :session kalman-en :exports code :results none
  def kalman_smoother(res, T):
      a_smooth, P_smooth = res['a_filt'].copy(), res['P_filt'].copy()
      for t in range(len(a_smooth) - 2, -1, -1):
          J = np.linalg.solve(res['P_pred'][t+1], T @ res['P_filt'][t]).T   # P_{t|t} T' P_{t+1|t}^{-1}
          a_smooth[t] = res['a_filt'][t] + J @ (a_smooth[t+1] - res['a_pred'][t+1])
          P_smooth[t] = res['P_filt'][t] + J @ (P_smooth[t+1] - res['P_pred'][t+1]) @ J.T
      return a_smooth, P_smooth

  a_smooth, P_smooth = kalman_smoother(res_sim, T)
  print(rmse(a_smooth[:, 0]))
  print(res_sim['P_filt'][50, 0, 0], P_smooth[50, 0, 0])
#+end_src

On the simulated series, the smoother reduces the root mean squared error
from $0.61$ to $0.48$. In the middle of the sample, the smoothed variance
is $0.243$ against $0.390$ for the filtered variance: the confidence band
is about a fifth narrower, except at the ends of the sample, where the
smoother does not have observations on both sides and joins the filter at
date $n$. The figure [[fig:niveau-local-lisseur][below]] also shows that
the smoother corrects the lag of the filter: it anticipates changes in the
level, since it sees the following observations.

#+begin_src python :session kalman-en :exports none :results none
  s_s = np.sqrt(P_smooth[:, 0, 0])
  fig, ax = plt.subplots(figsize=(6, 4))
  ax.plot(t_sim, mu, 'k', linewidth=0.8, label=r'$\mu_t$')
  ax.plot(t_sim, a_f, 'b', linewidth=1, label=r'$a_{t|t}$')
  ax.fill_between(t_sim, a_f-2*s_f, a_f+2*s_f, color='b', alpha=0.12, linewidth=0)
  ax.plot(t_sim, a_smooth[:, 0], 'r', linewidth=1, label=r'$a_{t|n}$')
  ax.fill_between(t_sim, a_smooth[:, 0]-2*s_s, a_smooth[:, 0]+2*s_s, color='r', alpha=0.2, linewidth=0)
  ax.set_xlabel(r'$t$')
  ax.legend()
  fig.tight_layout()
  fig.savefig("kalman-niveau-local-lisseur.svg", transparent=True)
#+end_src

#+CAPTION: *Simulated level, filtered state and smoothed state, with their bands of plus or minus two conditional standard deviations.*
#+LABEL: fig:niveau-local-lisseur
[[file:kalman-niveau-local-lisseur.svg]]

* Maximum likelihood estimation
:PROPERTIES:
:CUSTOM_ID: estimation
:END:

The system matrices generally depend on unknown parameters, gathered in a
vector $\psi$: the variances of the local level model, the coefficients and
the variance of an ARMA process. The Kalman filter delivers the likelihood
of these parameters at almost no extra cost.

** Prediction error decomposition
:PROPERTIES:
:CUSTOM_ID: decomposition
:END:

#+BEGIN_proposition
The log-likelihood of the observations $y_1, \dots, y_n$ is
\begin{equation*}
\log L(\psi) = -\frac{np}{2}\log 2\pi - \frac12\sum_{t=1}^n \log\det F_t - \frac12\sum_{t=1}^n v_t'F_t^{-1}v_t,
\end{equation*}
where the innovations $v_t$ and their variances $F_t$ are computed by the
filter with the value $\psi$ of the parameters.
#+END_proposition

#+BEGIN_proof
The joint density factorises as a product of conditional densities,
\begin{equation*}
f(y_1, \dots, y_n) = \prod_{t=1}^n f(y_t \mid Y_{t-1}),
\end{equation*}
with the convention $f(y_1 \mid Y_0) = f(y_1)$. By the proof of the filter,
the distribution of $y_t$ given $Y_{t-1}$ is normal, with mean
$Z_t a_{t|t-1} + d_t$ and variance $F_t$, and its density at $y_t$ is
\begin{equation*}
f(y_t \mid Y_{t-1}) = (2\pi)^{-p/2}(\det F_t)^{-1/2}\exp\left(-\frac12 v_t'F_t^{-1}v_t\right).
\end{equation*}
It remains to take the logarithm of the product.
#+END_proof

Computing the likelihood directly would require inverting the variance
matrix of the sample, of size $np \times np$, at a cost of order $n^3p^3$.
The filter implicitly computes a triangular factorisation of it, at a cost
proportional to $n$. For an ARMA process initialised with its stationary
distribution, this is the /exact/ likelihood, which accounts for the
distribution of the first observations, and not the likelihood conditional
on the first observations that least squares maximise (see the note on the
[[https://stephane-adjemian.fr/en/posts/ols-bias-in-an-ar1-process/][bias of the least squares estimator of an AR(1)]]
and the [[ex-h-nul][exercise]] on the autoregressive case). The code allows
the first $d$ innovations to be discarded, for the reason given below:

#+begin_src python :session kalman-en :exports code :results none
  def loglik(res, d=0):
      v, F = res['v'][d:, 0], res['F'][d:, 0, 0]
      ok = ~np.isnan(v)
      return -0.5*np.sum(np.log(2*np.pi*F[ok]) + v[ok]**2/F[ok])
#+end_src

** Concentrated likelihood
:PROPERTIES:
:CUSTOM_ID: concentree
:END:

In unobserved components models, the parameters are often variances, and
one of them can be eliminated analytically. We assume $p = 1$ and that the
parameters split into a scale factor $\sigma^2$ and a vector $\psi$, with
$H = \sigma^2H_*(\psi)$, $Q = \sigma^2Q_*(\psi)$ and
$P_0 = \sigma^2P_{0*}(\psi)$. We denote by $v_t^*$ and $F_t^*$ the
quantities computed by the filter with $\sigma^2 = 1$, and by $N$ the
number of innovations included in the likelihood.

#+BEGIN_property
The innovations do not depend on $\sigma^2$, $v_t = v_t^*$, and
$F_t = \sigma^2F_t^*$. The log-likelihood is maximised in $\sigma^2$ at
\begin{equation*}
\hat\sigma^2(\psi) = \frac1N\sum_t \frac{(v_t^*)^2}{F_t^*},
\end{equation*}
and the /concentrated log-likelihood/, obtained by replacing $\sigma^2$
with $\hat\sigma^2(\psi)$, is
\begin{equation*}
\log L_c(\psi) = -\frac N2\left(\log 2\pi + 1\right) - \frac12\sum_t\log F_t^* - \frac N2\log\hat\sigma^2(\psi).
\end{equation*}
#+END_property

#+BEGIN_proof
By induction, if $P_{t-1|t-1} = \sigma^2P_{t-1|t-1}^*$, the filter
equations give $P_{t|t-1} = \sigma^2P_{t|t-1}^*$, $F_t = \sigma^2F_t^*$ and
\(K_t = P_{t|t-1}Z'F_t^{-1} = K_t^*\): the gain, and therefore the predicted
states and the innovations, do not depend on $\sigma^2$, while
$P_{t|t} = \sigma^2P_{t|t}^*$. The log-likelihood then reads
\begin{equation*}
\log L = -\frac N2\log 2\pi - \frac N2\log\sigma^2 - \frac12\sum_t\log F_t^* - \frac{1}{2\sigma^2}\sum_t\frac{(v_t^*)^2}{F_t^*}.
\end{equation*}
As a function of $1/\sigma^2$, it is concave, and its derivative vanishes
at $\hat\sigma^2(\psi)$. Substituting, the last term equals $-N/2$.
#+END_proof

With the approximate diffuse initialisation, $P_0 = \kappa I_m$ is not
proportional to $\sigma^2$, but the discrepancy is of order $1/\kappa$.
More importantly, the first $d$ innovations, whose variance is of order
$\kappa$, contribute to the likelihood through terms $-\frac12\log F_t$
that depend on the arbitrary value of $\kappa$: they are discarded, and the
sums run over dates $t = d+1, \dots, n$, with $N = n - d$. This is the
/diffuse likelihood/ (Harvey, 1989, section 3.4.3). It differs from the
likelihood computed with the exact diffuse initialisation by a constant
that does not depend on the parameters, provided the diffuse part of the
model, in $Z$ and $T$, does not depend on them either, as in the trend
models of this note (Francke, Koopman and De Vos, /Journal of Time Series
Analysis/, 2010): both then give the same estimates, as we shall check in
the application.

#+begin_src python :session kalman-en :exports code :results none
  def concentrated_loglik(res, d=0):
      v, F = res['v'][d:, 0], res['F'][d:, 0, 0]
      ok = ~np.isnan(v)
      N = ok.sum()
      s2 = np.mean(v[ok]**2/F[ok])
      return -0.5*N*(np.log(2*np.pi) + 1) - 0.5*np.sum(np.log(F[ok])) - 0.5*N*np.log(s2), s2
#+end_src

** Precision of the estimators and tests
:PROPERTIES:
:CUSTOM_ID: erreurs-standard
:END:

Under regularity conditions, which require in particular that the true
value of the parameters lies in the interior of the parameter space, the
maximum likelihood estimator is asymptotically normal, and its variance
matrix is estimated by the inverse of minus the Hessian of the
log-likelihood at the maximum. The Hessian is computed by central finite
differences:
\begin{equation*}
\begin{aligned}
\frac{\partial^2 \log L}{\partial\psi_i\partial\psi_j} \approx \frac{1}{4h^2}\Big[
&\log L(\psi + he_i + he_j) - \log L(\psi + he_i - he_j)\\
&- \log L(\psi - he_i + he_j) + \log L(\psi - he_i - he_j)\Big],
\end{aligned}
\end{equation*}
where $e_i$ is the \(i\)-th vector of the canonical basis.

#+begin_src python :session kalman-en :exports code :results none
  def hessian(f, x, h=1e-3):
      x = np.atleast_1d(np.asarray(x, dtype=float))
      E = h*np.eye(len(x))
      return np.array([[(f(x+ei+ej) - f(x+ei-ej) - f(x-ei+ej) + f(x-ei-ej))/(4*h**2)
                        for ej in E] for ei in E])
#+end_src

The variance parameters are estimated in logarithms, which removes the
positivity constraint. Confidence intervals are built for the logarithm and
then transformed, which makes them asymmetric. The curvature of the
concentrated log-likelihood directly gives the precision of $\hat\psi$: the
inverse of minus its Hessian equals the corresponding block of the inverse
of minus the Hessian of the full
log-likelihood (Patefield, /Biometrika/, 1985). Finally, two
nested models
are compared by the likelihood ratio test,
\(\mathrm{LR} = 2(\log L_1 - \log L_0)\), whose asymptotic distribution is a
$\chi^2$ with as many degrees of freedom as restrictions, provided the
restricted model does not put a parameter on the boundary of the parameter
space, such as a zero variance.

* Forecasting
:PROPERTIES:
:CUSTOM_ID: prevision
:END:

#+BEGIN_proposition
For $h \geq 1$, the distribution of $\alpha_{n+h}$ given $Y_n$ is
$\mathcal N(a_{n+h|n}, P_{n+h|n})$, where
\begin{equation*}
a_{n+h|n} = T a_{n+h-1|n} + c,\qquad P_{n+h|n} = T P_{n+h-1|n} T' + RQR',
\end{equation*}
starting from $a_{n|n}$ and $P_{n|n}$. The forecast of the observation is
\(\hat y_{n+h|n} = Z a_{n+h|n} + d\), with variance
\(F_{n+h|n} = Z P_{n+h|n} Z' + H\), and, in the scalar case, a forecast
interval with level $1 - \delta$ is
\(\hat y_{n+h|n} \pm z_{1-\delta/2}\sqrt{F_{n+h|n}}\), where
$z_{1-\delta/2}$ is a quantile of the standard normal distribution.
#+END_proposition

#+BEGIN_proof
No observation is available after date $n$: the distribution of
$\alpha_{n+h}$ given $Y_n$ is obtained by applying $h$ times the prediction
step of the filter, without updating, and the distribution of $y_{n+h}$
follows from the measurement equation.
#+END_proof

In other words, forecasting amounts to filtering a series extended with
missing observations at dates $n+1, \dots, n+h$. With constant matrices, we
get \(a_{n+h|n} = T^h a_{n|n} + \sum_{j=0}^{h-1}T^jc\) and
\begin{equation*}
P_{n+h|n} = T^h P_{n|n} (T')^h + \sum_{j=0}^{h-1} T^j RQR' (T')^j.
\end{equation*}

#+BEGIN_exemple
In the smooth trend model,
\(T^h = \begin{pmatrix} 1 & h\\ 0 & 1 \end{pmatrix}\) and \(ZT^jR = j\).
The forecast extends the filtered trend linearly,
\(\hat y_{n+h|n} = a_{\mu, n|n} + h\,a_{\beta, n|n}\), where
$a_{\mu, n|n}$ and $a_{\beta, n|n}$ are the components of $a_{n|n}$. Since
\(\mu_{n+h} = \mu_n + h\beta_n + \sum_{j=1}^{h-1} j\,\zeta_{n+h-j}\), its
variance is
\begin{equation*}
F_{n+h|n} = \sigma_\varepsilon^2 + \begin{pmatrix} 1 & h \end{pmatrix} P_{n|n} \begin{pmatrix} 1\\ h \end{pmatrix}
+ \sigma_\zeta^2\sum_{j=1}^{h-1} j^2.
\end{equation*}
The last term, equal to $\sigma_\zeta^2(h-1)h(2h-1)/6$, grows like
$\sigma_\zeta^2h^3/3$: the forecast intervals of a trend integrated of
order two widen like $h^{3/2}$.
#+END_exemple

#+begin_src python :session kalman-en :exports code :results none
  def forecast(a, P, Z, T, R, H, Q, h):
      RQR = R @ Q @ R.T
      y_fc, F_fc = np.zeros(h), np.zeros(h)
      for j in range(h):
          a = T @ a
          P = T @ P @ T.T + RQR
          y_fc[j], F_fc[j] = (Z @ a)[0], (Z @ P @ Z.T + H)[0, 0]
      return y_fc, F_fc
#+end_src

* The Hodrick-Prescott filter as a Kalman smoother
:PROPERTIES:
:CUSTOM_ID: application
:END:

The Hodrick and Prescott filter (/Journal of Money, Credit and Banking/,
1997, from a 1980 working paper) is the most widely used tool to separate
the trend and the cycle of a macroeconomic series. We show that it
coincides with the Kalman smoother of a smooth trend model, which makes it
possible to estimate its parameter by maximum likelihood, and then apply
the result to French GDP per capita.

** The Hodrick-Prescott filter
:PROPERTIES:
:CUSTOM_ID: filtre-hp
:END:

#+BEGIN_definition
For a series $y_1, \dots, y_n$ and a parameter $\lambda > 0$, the
Hodrick-Prescott trend is the sequence $\hat\mu_1, \dots, \hat\mu_n$ that
minimises
\begin{equation*}
S(\mu) = \sum_{t=1}^n (y_t - \mu_t)^2 + \lambda\sum_{t=3}^n \left(\Delta^2\mu_t\right)^2,
\qquad \Delta^2\mu_t = \mu_t - 2\mu_{t-1} + \mu_{t-2}.
\end{equation*}
The cycle is the deviation $y_t - \hat\mu_t$.
#+END_definition

The criterion trades off the fit to the data against the smoothness of the
trend, measured by its second differences. For $\lambda = 0$, the trend is
the series itself; as $\lambda$ tends to infinity, it tends to a straight
line (see the [[ex-limites-hp][exercise]]).

#+BEGIN_property
The Hodrick-Prescott trend is the unique solution of the linear system
\((I_n + \lambda D'D)\hat\mu = y\), where $D$ is the $(n-2)\times n$ matrix
whose row $t$ contains the coefficients $1$, $-2$ and $1$ in columns $t$,
$t+1$ and $t+2$, and zeros elsewhere.
#+END_property

#+BEGIN_proof
The criterion reads \(S(\mu) = (y - \mu)'(y - \mu) + \lambda\mu'D'D\mu\). It
is a quadratic function with Hessian $2(I_n + \lambda D'D)$, positive
definite since $D'D$ is positive semi-definite: the criterion is strictly
convex, and its unique minimum sets the gradient
\(-2(y - \mu) + 2\lambda D'D\mu\) to zero.
#+END_proof

The matrix $I_n + \lambda D'D$ is pentadiagonal: the system is solved with a
sparse solver, in a number of operations proportional to $n$.

#+begin_src python :session kalman-en :exports code :results none
  def hp_filter(y, lam):
      n = len(y)
      D = sparse.diags([1.0, -2.0, 1.0], [0, 1, 2], shape=(n-2, n))
      return spsolve((sparse.eye(n) + lam*D.T @ D).tocsc(), y)
#+end_src

#+BEGIN_remarque
Applied to an infinite series, the filter is a symmetric moving average
whose gain for the trend is \(1/(1 + 16\lambda\sin^4(\omega/2))\) at
frequency $\omega$ (King and Rebelo, /Journal of Economic Dynamics and
Control/, 1993). This gain equals one half for cycles of period
$2\pi/\omega$ such that \(16\lambda\sin^4(\omega/2) = 1\): shorter
fluctuations are attributed to the cycle, longer ones to the trend. The
usual value $\lambda = 1\,600$ for quarterly data puts this period at about
forty quarters, that is about ten years. For annual data, $\lambda = 100$,
used by Backus and Kehoe (/American Economic Review/, 1992), puts it at
about twenty years, and $\lambda = 6.25$, recommended by Ravn and Uhlig
(/Review of Economics and Statistics/, 2002) to match the usual quarterly
filter, at about ten years.
#+END_remarque

** The smooth trend model
:PROPERTIES:
:CUSTOM_ID: tendance-lisse
:END:

#+BEGIN_theorem
Consider the smooth trend model
\begin{equation*}
y_t = \mu_t + \varepsilon_t,\qquad
\mu_t = \mu_{t-1} + \beta_{t-1},\qquad
\beta_t = \beta_{t-1} + \zeta_t,
\end{equation*}
with $\varepsilon_t \sim \mathcal N(0, \sigma_\varepsilon^2)$,
$\zeta_t \sim \mathcal N(0, \sigma_\zeta^2)$ and a diffuse prior on the
initial level and slope. The smoothed trend $\mathbb E[\mu_t \mid Y_n]$ is
the Hodrick-Prescott trend with parameter
\(\lambda = \sigma_\varepsilon^2/\sigma_\zeta^2\).
#+END_theorem

#+BEGIN_proof
The model equations give, for $t \geq 3$,
\begin{equation*}
\Delta^2\mu_t = (\mu_t - \mu_{t-1}) - (\mu_{t-1} - \mu_{t-2}) = \beta_{t-1} - \beta_{t-2} = \zeta_{t-1}.
\end{equation*}
The trend $(\mu_1, \dots, \mu_n)$ is therefore the image of
$(\mu_1, \beta_1, \zeta_2, \dots, \zeta_{n-1})$ under the linear map
defined by $\mu_2 = \mu_1 + \beta_1$ and
\(\mu_t = 2\mu_{t-1} - \mu_{t-2} + \zeta_{t-1}\) for $t \geq 3$. This map is
bijective, with a unit triangular matrix, and its Jacobian is one. With a
constant prior density for $(\mu_1, \beta_1)$, the limit of a normal
distribution whose variance tends to infinity, and independent $\zeta_t$
with distribution $\mathcal N(0, \sigma_\zeta^2)$, the prior density of the
trend is
\begin{equation*}
p(\mu_1, \dots, \mu_n) \propto \exp\left(-\frac{1}{2\sigma_\zeta^2}\sum_{t=3}^n\left(\Delta^2\mu_t\right)^2\right).
\end{equation*}
Given the trend, the observations are independent, with distribution
$\mathcal N(\mu_t, \sigma_\varepsilon^2)$. By Bayes' formula, the posterior
density of the trend is proportional to the product of the likelihood and
the prior density:
\begin{equation*}
p(\mu_1, \dots, \mu_n \mid Y_n) \propto
\exp\left(-\frac{1}{2\sigma_\varepsilon^2}\left[\sum_{t=1}^n(y_t - \mu_t)^2 + \lambda\sum_{t=3}^n\left(\Delta^2\mu_t\right)^2\right]\right).
\end{equation*}
This is the density of a normal distribution. Its mean equals its mode,
which minimises the bracket, that is the Hodrick-Prescott trend. Now the
Kalman smoother, with the diffuse initialisation, computes precisely the
posterior mean $\mathbb E[\mu_t \mid Y_n]$. The slope $\beta_n$, which does
not enter the observations, plays no role.
#+END_proof

This result extends the Bayesian interpretation of smoothing splines
(Wahba, /Journal of the Royal Statistical Society B/, 1978). Harvey and
Jaeger (/Journal of Applied Econometrics/, 1993) drew its consequences for
business cycle analysis: the Hodrick-Prescott trend depends only on the
ratio of the two variances, and the filter is optimal if the series is
generated by the smooth trend model with
\(\sigma_\varepsilon^2/\sigma_\zeta^2 = \lambda\). In return, the model
provides what the filter does not: confidence bands around the trend, a
one-sided trend, computed by the filter, that uses only the past, forecasts
and above all a likelihood, which makes it possible to estimate $\lambda$
instead of fixing it.

#+BEGIN_remarque
The smoothed trend is not a property of the distribution of the
observations, but of the chosen decomposition. The smooth trend model
assumes that the shocks to the trend and those to the deviation from trend
are independent; other assumptions are compatible with the same
distribution of $y_t$ and give other trends. At the other extreme, the
Beveridge and Nelson decomposition (/Journal of Monetary Economics/, 1981)
defines the trend as the long-run forecast of the series, net of trend
growth, which amounts to assuming that the shocks to the trend and to the
cycle are perfectly correlated: the trend is then a random walk. Morley,
Nelson and Zivot (/Review of Economics and Statistics/, 2003) note that,
for US GDP, this trend accounts for most of the fluctuations, and show that
Clark's trend and cycle model,
estimated without imposing independent shocks, gives filtered estimates of
the trend and the cycle identical to those of Beveridge and Nelson. For US
quarterly GDP from 1947 to 1998, the estimated correlation between the
shocks is $-0.9$, and independence is rejected at the 5% level, narrowly.
For a random walk with drift, which the data favour
[[#tendance-lineaire-locale][below]] for French GDP per capita over the
long run, the Beveridge-Nelson trend is the series itself and the cycle is
zero. The Hodrick-Prescott trend is therefore only one decomposition among
others, set by two assumptions: independent components and the shape of
the trend, whose smoothness is governed by $\lambda$.
#+END_remarque

** The data
:PROPERTIES:
:CUSTOM_ID: donnees
:END:

We use French GDP per capita from the Maddison Project Database, 2020
version, downloaded as in the note on
[[https://stephane-adjemian.fr/en/posts/plotting-a-gdp-time-series/][plotting a GDP time series]].
The annual series is continuous from 1820 to 2018, that is $n = 199$
observations. We work with \(y_t = 100\log(\text{GDP per capita})\), so that
first differences are growth rates in percent and standard deviations are
in percentage points.

#+begin_src python :session kalman-en :exports code :results none
  MADDISON = 'mpd2020'
  if not os.path.exists(MADDISON+'.xlsx'):
      url.urlretrieve('https://www.rug.nl/ggdc/historicaldevelopment/maddison/data/'+MADDISON+'.xlsx',
                      MADDISON+'.xlsx')
  data = pd.read_excel(MADDISON+'.xlsx', 'Full data')
  france = data[(data.countrycode == 'FRA') & (data.year >= 1820)]
  years = france.year.values.astype(int)
  y = 100*np.log(france.gdppc.values)
  dy = np.diff(y)
  for start, end in ((1820, 2018), (1820, 1913), (1950, 1973), (1974, 2018)):
      print(start, end, (y[years == end] - y[years == start])/(end - start))
  print(years[1:][np.argsort(dy)[:2]], np.sort(dy)[:2], years[1:][np.argmax(dy)], dy.max())
#+end_src

Average annual growth is $1.54$ points over the whole period: $1.21$ points
per year between 1820 and 1913, $3.94$ between 1950 and 1973 and $1.39$
since 1974. The two world wars leave considerable falls, of $21.8$ points
in 1918 and $20.0$ points in 1941, followed by large rebounds, up to
$40.4$ points in 1946. These exceptional years weigh on all the results
that follow.

** Numerical check
:PROPERTIES:
:CUSTOM_ID: verification
:END:

We compare the Hodrick-Prescott trend, with $\lambda = 100$, and the
smoothed trend of the model with $\sigma_\varepsilon^2 = 1$ and
$\sigma_\zeta^2 = 1/100$, for several values of $\kappa$:

#+begin_src python :session kalman-en :exports code :results none
  lam = 100
  hp_trend = hp_filter(y, lam)

  def smoothed_trend(y, s2_eps, s2_zeta, kappa):
      Z, T, R, H, Q = smooth_trend(s2_eps, s2_zeta)
      res = kalman_filter(y, Z, T, R, H, Q, np.zeros(2), kappa*np.eye(2))
      return res, kalman_smoother(res, T)

  for kap in (1e3, 1e5, 1e7, 1e9, 1e11):
      _, (a_s, _) = smoothed_trend(y, 1.0, 1/lam, kap)
      print(kap, np.max(np.abs(hp_trend - a_s[:, 0])))
  kappa = 1e9
#+end_src

| $\kappa$  | Maximum difference   |
|-----------+----------------------|
| $10^3$    | $3.3\times 10^{-1}$  |
| $10^5$    | $3.3\times 10^{-3}$  |
| $10^7$    | $3.3\times 10^{-5}$  |
| $10^9$    | $3.2\times 10^{-7}$  |
| $10^{11}$ | $2.2\times 10^{-5}$  |

The difference decreases like $1/\kappa$ as long as the diffuse
approximation dominates, then rises again when rounding errors take over:
with $\kappa = 10^{11}$, computing the first filtered variances loses about
ten significant digits. It is largest at the beginning of the sample, where
the initialisation matters, as shown in the figure [[fig:pib-france][below]].
From now on we use $\kappa = 10^9$.

#+begin_src python :session kalman-en :exports none :results none
  res_hp, (a_hp, P_hp) = smoothed_trend(y, 1.0, 1/lam, kappa)
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  axes[0].plot(years, y, 'k', linewidth=0.8, label=r'$y_t$')
  axes[0].plot(years, hp_trend, 'b', linewidth=1.5, label='Hodrick-Prescott')
  axes[0].plot(years, a_hp[:, 0], 'r--', linewidth=1.5, label='Kalman smoother')
  axes[0].set_xlabel('year')
  axes[0].legend()
  axes[1].plot(years, 1e7*(hp_trend - a_hp[:, 0]), 'b', linewidth=1)
  axes[1].axhline(0, color='k', linewidth=0.5)
  axes[1].set_xlabel('year')
  axes[1].set_ylabel(r'difference ($\times 10^{-7}$)')
  fig.tight_layout()
  fig.savefig("kalman-pib-france.svg", transparent=True)
#+end_src

#+CAPTION: *Left, logarithm of French GDP per capita (multiplied by $100$), Hodrick-Prescott trend with $\lambda = 100$ and smoothed trend of the smooth trend model with $\kappa = 10^9$. Right, difference between the two trends.*
#+LABEL: fig:pib-france
[[file:kalman-pib-france.svg]]

The external check uses the =UnobservedComponents= class of =statsmodels=,
which implements the exact diffuse initialisation. We compare the smoothed
trends and the log-likelihoods for three pairs of variances:

#+begin_src python :session kalman-en :exports code :results none
  from statsmodels.tsa.statespace.structural import UnobservedComponents
  model_sm = UnobservedComponents(y, level='smooth trend', use_exact_diffuse=True)
  for s2e, s2z in ((1.0, 0.01), (10.0, 5.0), (40.0, 0.4)):
      res, (a_s, _) = smoothed_trend(y, s2e, s2z, kappa)
      trend_sm = model_sm.smooth([s2e, s2z]).smoothed_state[0]
      print(np.max(np.abs(trend_sm - a_s[:, 0])),
            model_sm.loglike([s2e, s2z]) - loglik(res, d=2))
#+end_src

The trends differ by less than $1.4\times 10^{-5}$, and the difference
between the log-likelihoods is $-1.83787$ within $2\times 10^{-5}$ for the
three pairs, that is $-\log 2\pi$: as announced, the diffuse likelihood and
the exact diffuse likelihood differ only by a constant, which affects
neither the estimates nor the tests.

** Estimating the smoothing parameter
:PROPERTIES:
:CUSTOM_ID: estimer-lambda
:END:

The parameter $\lambda$ is the inverse of the signal-to-noise ratio
\(q = \sigma_\zeta^2/\sigma_\varepsilon^2\). We concentrate the likelihood
with respect to $\sigma_\varepsilon^2$ and maximise it in $\log q$,
discarding the first $d = 2$ innovations, since the initial level and slope
are diffuse. We also compute the likelihood with $\lambda = 100$, and the
estimate over two subperiods.

#+begin_src python :session kalman-en :exports code :results none
  def logLc_smooth_trend(log_q, y):
      Z, T, R, H, Q = smooth_trend(1.0, np.exp(log_q))
      res = kalman_filter(y, Z, T, R, H, Q, np.zeros(2), kappa*np.eye(2))
      return concentrated_loglik(res, d=2)

  def estimate_lambda(y):
      opt = minimize_scalar(lambda x: -logLc_smooth_trend(x, y)[0],
                            bounds=(np.log(1e-6), np.log(1e3)), method='bounded',
                            options={'xatol': 1e-8})
      return opt.x

  log_q = estimate_lambda(y)
  logL_max, s2_eps = logLc_smooth_trend(log_q, y)
  lam_hat = np.exp(-log_q)
  se = 1/np.sqrt(-hessian(lambda x: logLc_smooth_trend(x[0], y)[0], [log_q])[0, 0])
  ci = np.exp(-log_q + np.array([-1.96, 1.96])*se)
  logL_100, s2_eps_100 = logLc_smooth_trend(np.log(1/100), y)
  LR = 2*(logL_max - logL_100)
  print(lam_hat, se, ci, LR, chi2.sf(LR, 1))
  print(2*np.pi/(2*np.arcsin((16*lam_hat)**-0.25)))
  for start in (1870, 1950):
      print(start, np.exp(-estimate_lambda(y[years >= start])))
#+end_src

The estimate is $\hat\lambda = 1.24$, with a 95% confidence interval of
$[0.65, 2.36]$, built from the standard error of $\log\hat q$, which is
$0.33$. This is very far from the usual values: the likelihood ratio test
rejects $\lambda = 100$ with a statistic of $93.8$, for a critical value of
$3.84$ at the 5% level. By the remark on the frequency response, the gain
of the trend equals one half for a period of $6.4$ years with
$\hat\lambda$, against about twenty years with $\lambda = 100$: the
estimated trend absorbs most medium-term fluctuations. The result is not
due to the wars alone: the estimate is $0.86$ from 1870 and $0.53$
from 1950. The figure [[fig:vraisemblance][below]] shows the concentrated
log-likelihood as a function of $\lambda$: $\lambda = 6.25$ is already
clearly less likely than $\hat\lambda$, and $\lambda = 100$ or
$\lambda = 1\,600$ much less so.

To understand this rejection, we compare the two models, with the estimated
variances, using the standardised innovations and the cycle
$y_t - a_{\mu, t|n}$:

#+begin_src python :session kalman-en :exports code :results none
  def fit(y, s2_eps, lam):
      res, (a_s, P_s) = smoothed_trend(y, s2_eps, s2_eps/lam, kappa)
      e = res['v'][2:, 0]/np.sqrt(res['F'][2:, 0, 0])
      return res, a_s, P_s, e, y - a_s[:, 0]

  results = {'100': fit(y, s2_eps_100, 100), 'hat': fit(y, s2_eps, lam_hat)}
  for name, (res, a_s, P_s, e, cycle) in results.items():
      print(name, acf(e, 3), cycle.std(), acf(cycle, 1))
#+end_src

|                                     | $\lambda = 100$ | $\hat\lambda = 1.24$ |
|-------------------------------------+-----------------+----------------------|
| $\hat\sigma_\varepsilon$            | $6.85$          | $3.40$               |
| $\hat\sigma_\zeta$                  | $0.68$          | $3.06$               |
| Log-likelihood                      | $-704.4$        | $-657.5$             |
| Autocorrelation of innovations (1)  | $0.66$          | $0.08$               |
| Autocorrelation of innovations (2)  | $0.34$          | $-0.26$              |
| Standard deviation of the cycle     | $5.94$          | $2.70$               |
| Autocorrelation of the cycle (1)    | $0.57$          | $-0.22$              |

#+begin_src python :session kalman-en :exports none :results none
  lam_grid = np.logspace(-2, 4, 121)
  profile = [logLc_smooth_trend(-np.log(l), y)[0] for l in lam_grid]
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  axes[0].plot(lam_grid, profile, 'b', linewidth=1)
  for l, c, label in ((lam_hat, 'r', r'$\hat\lambda$'), (6.25, 'g', r'$\lambda = 6.25$'),
                      (100, 'k', r'$\lambda = 100$')):
      axes[0].axvline(l, color=c, linestyle='--', linewidth=0.8, label=label)
  axes[0].set_xscale('log')
  axes[0].set_xlabel(r'$\lambda$')
  axes[0].set_ylabel('concentrated log-likelihood')
  axes[0].legend()
  k = np.arange(1, 11)
  axes[1].bar(k - 0.17, acf(results['100'][3], 10), width=0.34, color='k', label=r'$\lambda = 100$')
  axes[1].bar(k + 0.17, acf(results['hat'][3], 10), width=0.34, color='r', label=r'$\hat\lambda$')
  band = 1.96/np.sqrt(len(results['hat'][3]))
  for s in (-1, 1):
      axes[1].axhline(s*band, color='k', linestyle='--', linewidth=0.8)
  axes[1].axhline(0, color='k', linewidth=0.5)
  axes[1].set_xticks(k)
  axes[1].set_xlabel('lag')
  axes[1].set_ylabel('autocorrelation of innovations')
  axes[1].legend()
  fig.tight_layout()
  fig.savefig("kalman-vraisemblance.svg", transparent=True)
#+end_src

#+CAPTION: *Left, concentrated log-likelihood of the smooth trend model as a function of $\lambda$ (logarithmic scale). Right, autocorrelations of the standardised innovations for $\lambda = 100$ and for $\hat\lambda$, with the band $\pm 1.96/\sqrt{n-2}$.*
#+LABEL: fig:vraisemblance
[[file:kalman-vraisemblance.svg]]

The two values of $\lambda$ tell very different stories (figure
[[fig:cycles][below]]). With $\lambda = 100$, the trend is smooth and the
cycle, with a standard deviation of $5.9$ points, is very persistent: it
records the wars as recessions lasting several years. But the model
assumes that the deviation from trend, $\varepsilon_t$, is white noise, and
its standardised innovations have a first-order autocorrelation of $0.66$:
the specification is rejected by its own diagnostics. Maximum likelihood
resolves the contradiction by making the trend flexible. With
$\hat\lambda$, the trend absorbs the persistent fluctuations and the
residual cycle, with a standard deviation of $2.7$ points, is no longer
persistent. The smooth trend model has no room for a persistent cycle: it
knows only the trend and the noise. Fixing $\lambda$ therefore amounts to
choosing a definition of the cycle, not to estimating it, a conclusion
already reached by Harvey and Jaeger (1993).

#+begin_src python :session kalman-en :exports none :results none
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  axes[0].plot(years, y, 'k', linewidth=0.8, label=r'$y_t$')
  axes[0].plot(years, results['100'][1][:, 0], 'b', linewidth=1.2, label=r'trend, $\lambda = 100$')
  axes[0].plot(years, results['hat'][1][:, 0], 'r', linewidth=1.2, label=r'trend, $\hat\lambda$')
  axes[0].set_xlabel('year')
  axes[0].legend()
  axes[1].plot(years, results['100'][4], 'b', linewidth=1, label=r'$\lambda = 100$')
  axes[1].plot(years, results['hat'][4], 'r', linewidth=1, label=r'$\hat\lambda$')
  axes[1].axhline(0, color='k', linewidth=0.5)
  axes[1].set_xlabel('year')
  axes[1].set_ylabel(r'$y_t - a_{\mu, t|n}$')
  axes[1].legend()
  fig.tight_layout()
  fig.savefig("kalman-cycles.svg", transparent=True)
#+end_src

#+CAPTION: *Smoothed trends (left) and cycles (right) for $\lambda = 100$ and for $\hat\lambda$.*
#+LABEL: fig:cycles
[[file:kalman-cycles.svg]]

*** Freeing the level variance
:PROPERTIES:
:CUSTOM_ID: tendance-lineaire-locale
:END:

The smooth trend model imposes $\sigma_\eta^2 = 0$. The local linear trend
model, which estimates all three variances, decides otherwise. We
concentrate the likelihood with respect to $\sigma_\eta^2$, the ratios of
the other two variances to $\sigma_\eta^2$ being estimated in logarithms and
bounded below by $e^{-20}$.

#+begin_src python :session kalman-en :exports code :results none
  def logLc_llt(x, y):
      Z, T, R, H, Q = local_linear_trend(np.exp(x[0]), 1.0, np.exp(x[1]))
      res = kalman_filter(y, Z, T, R, H, Q, np.zeros(2), kappa*np.eye(2))
      return concentrated_loglik(res, d=2)

  opt_llt = minimize(lambda x: -logLc_llt(x, y)[0], [0.0, -2.0], method='L-BFGS-B',
                     bounds=[(-20.0, 5.0)]*2)
  logL_llt, s2_eta = logLc_llt(opt_llt.x, y)
  print(opt_llt.x, s2_eta, logL_llt, 2*(logL_llt - logL_max))
#+end_src

The optimum lies on the boundary: both ratios are at, or very close to,
their lower bound, and $\hat\sigma_\eta = 5.99$. The model reduces to a
random walk with drift, $\Delta y_t = \beta + \eta_t$, whose constant slope
is the average growth rate, $1.54$ points per year, and whose shocks have
the standard deviation of annual growth. The log-likelihood, $-635.0$,
exceeds that of the smooth trend model by $22.6$. The likelihood ratio
statistic, $45.1$, does not follow a $\chi^2$ distribution, since the
hypothesis $\sigma_\eta^2 = 0$ puts a parameter on the boundary of the
parameter space. In the setting of independent observations studied by
Self and Liang (/Journal of the American Statistical Association/, 1987),
its asymptotic distribution would be an equal mixture of a point mass at
zero and a $\chi^2(1)$ distribution; for a state variance in a time series
model, it may differ. With a statistic of this size, however, the
conclusion is beyond doubt. The
random walk with drift, which has only one variance parameter, is far more
likely than the smooth trend model, which has two: each annual change in
GDP per capita is a permanent shock to the level, including the falls and
rebounds of the wars.

#+BEGIN_remarque
We find again the pile-up phenomenon at the boundary mentioned in the note
on the
[[https://stephane-adjemian.fr/en/posts/ols-bias-in-an-ar1-process/#etat-mesure][bias of the least squares estimator of an AR(1)]]:
in unobserved components models, estimated variances are exactly zero with
positive probability. More importantly, if the series is a random walk, the
Hodrick-Prescott filter creates a persistent cycle that the series does not
contain (Cogley and Nason, /Journal of Economic Dynamics and Control/,
1995): the cycle obtained with $\lambda = 100$ is partly such an artefact.
#+END_remarque

** Filtered and smoothed trends
:PROPERTIES:
:CUSTOM_ID: filtre-vs-lisseur
:END:

The filtered trend $a_{\mu, t|t}$ uses only past observations: it is the
one-sided, real-time version of the Hodrick-Prescott filter (Stock and
Watson, /Journal of Monetary Economics/, 1999). Its gain quickly reaches its
steady-state value, computed from the algebraic Riccati equation:

#+begin_src python :session kalman-en :exports code :results none
  for name, l in (('100', 100), ('hat', lam_hat)):
      Z, T, R, H, Q = smooth_trend(1.0, 1/l)
      P_bar = solve_discrete_are(T.T, Z.T, R @ Q @ R.T, H)
      K_bar = P_bar @ Z.T/(Z @ P_bar @ Z.T + H)
      K_t = results[name][0]['K'][:, :, 0]
      print(name, K_bar.ravel(), years[np.argmax(np.all(np.abs(K_t - K_bar.ravel()) < 1e-3, axis=1))])

  res, a_s, P_s, e, cycle = results['100']
  revision = res['a_filt'][:, 0] - a_s[:, 0]
  print(years[np.argmax(np.abs(revision))], revision.min())
  print(np.sqrt(res['P_filt'][years == 1930, 0, 0]), np.sqrt(P_s[years == 1930, 0, 0]))
  print(acf(y - res['a_filt'][:, 0], 1), acf(y - a_s[:, 0], 1))
#+end_src

With $\lambda = 100$, the steady-state gain is $(0.362, 0.080)'$, reached
within $10^{-3}$ from 1834: each year, the filtered trend corrects its level
by 36% of the innovation and its slope by 8%. With $\hat\lambda$, it is
$(0.750, 0.450)'$, reached from 1825. The figure
[[fig:filtre-vs-lisseur][below]] compares the two trends for
$\lambda = 100$ around the two wars. In real time, the collapses are partly
attributed to the trend, which drops; with hindsight, the smoother, which
sees the reconstruction, attributes them to the cycle. In 1945, the
filtered trend is $31.5$ points below the smoothed trend. The conditional
standard deviation of the trend is $4.1$ points for the filter and $2.3$
points for the smoother, away from the ends of the sample. The one-sided
cycle is more persistent than the two-sided cycle, with a first-order
autocorrelation of $0.66$ against $0.57$.

#+begin_src python :session kalman-en :exports none :results none
  window = (years >= 1900) & (years <= 1960)
  af, sf = res['a_filt'][:, 0], np.sqrt(res['P_filt'][:, 0, 0])
  asm, ssm = a_s[:, 0], np.sqrt(P_s[:, 0, 0])
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  axes[0].plot(years[window], y[window], 'k', linewidth=0.8, label=r'$y_t$')
  axes[0].plot(years[window], af[window], 'b', linewidth=1.2, label=r'filtered trend $a_{\mu, t|t}$')
  axes[0].fill_between(years[window], (af-2*sf)[window], (af+2*sf)[window], color='b', alpha=0.12, linewidth=0)
  axes[0].plot(years[window], asm[window], 'r', linewidth=1.2, label=r'smoothed trend $a_{\mu, t|n}$')
  axes[0].fill_between(years[window], (asm-2*ssm)[window], (asm+2*ssm)[window], color='r', alpha=0.2, linewidth=0)
  axes[0].set_xlabel('year')
  axes[0].legend(loc='upper left')
  axes[1].plot(years, y - af, 'b', linewidth=1, label='one-sided cycle')
  axes[1].plot(years, y - asm, 'r', linewidth=1, label='two-sided cycle')
  axes[1].axhline(0, color='k', linewidth=0.5)
  axes[1].set_xlabel('year')
  axes[1].legend()
  fig.tight_layout()
  fig.savefig("kalman-filtre-vs-lisseur.svg", transparent=True)
#+end_src

#+CAPTION: *Smooth trend model with $\lambda = 100$. Left, filtered and smoothed trends from 1900 to 1960, with bands of plus or minus two conditional standard deviations. Right, one-sided cycle, $y_t - a_{\mu, t|t}$, and two-sided cycle, $y_t - a_{\mu, t|n}$, which is the Hodrick-Prescott cycle.*
#+LABEL: fig:filtre-vs-lisseur
[[file:kalman-filtre-vs-lisseur.svg]]

#+BEGIN_remarque
The handling of missing observations makes it possible to interpolate the
trend. If the years 1939 to 1945 are removed, the smoother extends the
trend through the war: it estimates it at $880.8$ points in 1942, with a
standard deviation of $3.7$ points, against $846.6$ observed and a trend of
$861.4$ with all the data.
#+END_remarque

#+begin_src python :session kalman-en :exports code :results none
  y_missing = y.copy()
  y_missing[(years >= 1939) & (years <= 1945)] = np.nan
  _, (a_m, P_m) = smoothed_trend(y_missing, s2_eps_100, s2_eps_100/100, kappa)
  print(a_m[years == 1942, 0], np.sqrt(P_m[years == 1942, 0, 0]))
  print(y[years == 1942], results['100'][1][years == 1942, 0])
#+end_src

** A trend and a cycle
:PROPERTIES:
:CUSTOM_ID: tendance-cycle
:END:

To obtain a persistent cycle that is not imposed by the choice of
$\lambda$, it has to be modelled. We replace the noise $\varepsilon_t$ by a
stationary second-order autoregressive process. This is a variant of the
model of Clark (/Quarterly Journal of Economics/, 1987), whose trend also
has a shock on the level:
\begin{equation*}
y_t = \mu_t + \psi_t,\qquad
\mu_t = \mu_{t-1} + \beta_{t-1},\qquad
\beta_t = \beta_{t-1} + \zeta_t,\qquad
\psi_t = \phi_1\psi_{t-1} + \phi_2\psi_{t-2} + u_t,
\end{equation*}
with \(u_t \sim \mathcal N(0, \sigma_u^2)\). The state is
\(\alpha_t = (\mu_t, \beta_t, \psi_t, \psi_{t-1})'\), the measurement is
perfect ($H = 0$) and
\begin{equation*}
Z = \begin{pmatrix} 1 & 0 & 1 & 0 \end{pmatrix},\quad
T = \begin{pmatrix}
1 & 1 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & \phi_1 & \phi_2\\
0 & 0 & 1 & 0
\end{pmatrix},\quad
R = \begin{pmatrix} 0 & 0\\ 1 & 0\\ 0 & 1\\ 0 & 0 \end{pmatrix},\quad
Q = \begin{pmatrix} \sigma_\zeta^2 & 0\\ 0 & \sigma_u^2 \end{pmatrix}.
\end{equation*}
The initialisation is mixed: diffuse for the trend, stationary, through the
Lyapunov equation, for the cycle. We concentrate the likelihood with respect
to $\sigma_u^2$. To guarantee the stationarity of the cycle,
$(\phi_1, \phi_2)$ are parameterised by the partial autocorrelations
$r_1, r_2 \in (-1, 1)$, with $\phi_2 = r_2$ and \(\phi_1 = r_1(1 - r_2)\)
(Barndorff-Nielsen and Schou, /Journal of Multivariate Analysis/, 1973),
themselves written \(r_i = \tanh(x_i)\) with $x_i$ real.

#+begin_src python :session kalman-en :exports code :results none
  def trend_cycle(s2_zeta, phi1, phi2, s2_u):
      Z = np.array([[1.0, 0.0, 1.0, 0.0]])
      T = np.array([[1.0, 1.0, 0.0, 0.0],
                    [0.0, 1.0, 0.0, 0.0],
                    [0.0, 0.0, phi1, phi2],
                    [0.0, 0.0, 1.0, 0.0]])
      R = np.array([[0.0, 0.0],
                    [1.0, 0.0],
                    [0.0, 1.0],
                    [0.0, 0.0]])
      return Z, T, R, np.zeros((1, 1)), np.diag([s2_zeta, s2_u])

  def init_trend_cycle(T, R, Q):
      P0 = np.zeros((4, 4))
      P0[:2, :2] = kappa*np.eye(2)
      P0[2:, 2:] = solve_discrete_lyapunov(T[2:, 2:], (R @ Q @ R.T)[2:, 2:])
      return np.zeros(4), P0

  def parameters(x):
      r1, r2 = np.tanh(x[1]), np.tanh(x[2])
      return np.exp(x[0]), r1*(1 - r2), r2

  def logLc_trend_cycle(x, y):
      Z, T, R, H, Q = trend_cycle(*parameters(x), 1.0)
      res = kalman_filter(y, Z, T, R, H, Q, *init_trend_cycle(T, R, Q))
      return concentrated_loglik(res, d=2)
#+end_src

The wars make it impossible to estimate this model on the full sample: the
cycle absorbs the falls and rebounds, its dominant autoregressive root is
$0.95$ and the variance of the slope becomes almost zero, so that the model
gets close to a random walk with drift. We therefore estimate it over
1950--2018, from four starting points, because the likelihood of unobserved
components models often has several local maxima.

#+begin_src python :session kalman-en :exports code :results none
  y50, years50 = y[years >= 1950], years[years >= 1950]
  best = None
  for x0 in ([-2.0, 0.5, 0.0], [-5.0, 1.0, -0.5], [0.0, 0.2, 0.2], [-1.0, 1.5, -0.8]):
      opt = minimize(lambda x: -logLc_trend_cycle(x, y50)[0], x0, method='Nelder-Mead',
                     options={'xatol': 1e-8, 'fatol': 1e-10, 'maxiter': 5000})
      print(x0, -opt.fun)
      if best is None or opt.fun < best.fun:
          best = opt
  logL_tc, s2_u = logLc_trend_cycle(best.x, y50)
  q_zeta, phi1, phi2 = parameters(best.x)
  roots = np.roots([1, -phi1, -phi2])
  print(q_zeta*s2_u, phi1, phi2, s2_u, np.abs(roots), 2*np.pi/np.abs(np.angle(roots[0])))
  logL_st50 = logLc_smooth_trend(estimate_lambda(y50), y50)[0]
  print(logL_tc, logL_st50, chi2.sf(2*(logL_tc - logL_st50), 2))
#+end_src

One of the four optimisations stops at a local maximum where the first
partial autocorrelation of the cycle tends to one; the other three converge
to the same point. The estimates are $\hat\phi_1 = 0.807$ and
$\hat\phi_2 = -0.306$: the roots of the autoregressive polynomial are
complex, with modulus $0.55$, and the cycle has a period of $8.3$ years.
The standard deviation of the cycle shocks is $1.03$ points, that of the
slope of the trend $0.40$ points. The model clearly improves on the
log-likelihood of the smooth trend model estimated over the same period,
which goes from $-126.4$ to $-119.1$: the likelihood ratio statistic is
$14.6$ for the two restrictions $\phi_1 = \phi_2 = 0$, a p-value of
$0.0007$.

#+begin_src python :session kalman-en :exports code :results none
  Z_tc, T_tc, R_tc, H_tc, Q_tc = trend_cycle(q_zeta*s2_u, phi1, phi2, s2_u)
  res_tc = kalman_filter(y50, Z_tc, T_tc, R_tc, H_tc, Q_tc, *init_trend_cycle(T_tc, R_tc, Q_tc))
  a_tc, P_tc = kalman_smoother(res_tc, T_tc)
  hp_distance = lambda l: np.sqrt(np.mean((hp_filter(y50, l) - a_tc[:, 0])**2))
  opt = minimize_scalar(lambda x: hp_distance(np.exp(x)), bounds=(np.log(0.1), np.log(1e4)), method='bounded')
  print(np.exp(opt.x), hp_distance(np.exp(opt.x)), hp_distance(100))
  print(acf(res_tc['v'][2:, 0]/np.sqrt(res_tc['F'][2:, 0, 0]), 3))
  print(a_tc[-1, 2], (y50 - hp_filter(y50, 100))[-1])
#+end_src

The standardised innovations no longer show any autocorrelation. The trend
of the model is very close to a Hodrick-Prescott trend, but not with
parameter $100$: the value of $\lambda$ that best reproduces it is $27$,
with a root mean squared difference of $0.10$ points, against $0.39$ for
$\lambda = 100$. The smoothed cycle (figure [[fig:tendance-cycle][below]])
picks up the recessions of 1975, 1993 and 2009, and puts the economy
slightly above its trend in 2018, by $0.8$ points against $1.6$ for the
Hodrick-Prescott cycle.

#+begin_src python :session kalman-en :exports none :results none
  hp50 = hp_filter(y50, 100)
  c, sc = a_tc[:, 2], np.sqrt(P_tc[:, 2, 2])
  fig, axes = plt.subplots(1, 2, figsize=(12, 4))
  axes[0].plot(years50, y50, 'k', linewidth=0.8, label=r'$y_t$')
  axes[0].plot(years50, a_tc[:, 0], 'r', linewidth=1.2, label='model trend')
  axes[0].plot(years50, hp50, 'b--', linewidth=1.2, label=r'Hodrick-Prescott, $\lambda = 100$')
  axes[0].set_xlabel('year')
  axes[0].legend()
  axes[1].plot(years50, c, 'r', linewidth=1.2, label=r'cycle $a_{\psi, t|n}$')
  axes[1].fill_between(years50, c-2*sc, c+2*sc, color='r', alpha=0.2, linewidth=0)
  axes[1].plot(years50, y50 - hp50, 'b--', linewidth=1, label=r'HP cycle, $\lambda = 100$')
  axes[1].axhline(0, color='k', linewidth=0.5)
  axes[1].set_xlabel('year')
  axes[1].legend()
  fig.tight_layout()
  fig.savefig("kalman-tendance-cycle.svg", transparent=True)
#+end_src

#+CAPTION: *Trend and cycle model estimated over 1950--2018. Left, smoothed trend and Hodrick-Prescott trend with $\lambda = 100$. Right, smoothed cycle with a band of plus or minus two conditional standard deviations, and Hodrick-Prescott cycle.*
#+LABEL: fig:tendance-cycle
[[file:kalman-tendance-cycle.svg]]

** Forecasts
:PROPERTIES:
:CUSTOM_ID: previsions-pib
:END:

Finally, we compare the forecasts from 2019 to 2028 of two models estimated
over 1950--2018: the smooth trend model with $\lambda = 100$, the variance
$\sigma_\varepsilon^2$ being estimated, and the trend and cycle model.

#+begin_src python :session kalman-en :exports code :results none
  h = 10
  s2_100_50 = logLc_smooth_trend(np.log(1/100), y50)[1]
  model_st = smooth_trend(s2_100_50, s2_100_50/100)
  res_st = kalman_filter(y50, *model_st, np.zeros(2), kappa*np.eye(2))
  fc = {'st': forecast(res_st['a_filt'][-1], res_st['P_filt'][-1], *model_st, h),
        'tc': forecast(res_tc['a_filt'][-1], res_tc['P_filt'][-1], Z_tc, T_tc, R_tc, H_tc, Q_tc, h)}
  for name, (y_fc, F_fc) in fc.items():
      print(name, y_fc[0] - y50[-1], np.sqrt(F_fc[[0, 4, 9]]))
  print(res_st['a_filt'][-1][1], res_tc['a_filt'][-1][1], np.sqrt(s2_100_50/100), np.sqrt(q_zeta*s2_u))
#+end_src

With $\lambda = 100$, the forecast for 2019 is $1.0$ point below the 2018
level: the cycle, assumed to have no memory, disappears immediately and the
forecast starts from the filtered trend, which it extends with a slope of
$0.60$ points per year. The model with a cycle forecasts instead a growth
of $0.8$ points in 2019, since the cycle fades only gradually, followed by
a slope of $1.04$ points. The forecast standard deviations are $2.2$, $3.3$
and $5.6$ points at one, five and ten years with $\lambda = 100$, against
$1.4$, $4.9$ and $10.2$ points with the cycle (figure
[[fig:previsions][below]]). The apparent long-run precision of the first
model comes from the small variance it imposes on the slope of the trend,
$\sigma_\zeta = 0.17$ points against $0.40$: fixing $\lambda$ also means
fixing the uncertainty about future growth.

#+begin_src python :session kalman-en :exports none :results none
  years_fc = np.arange(years50[-1] + 1, years50[-1] + h + 1)
  window = years50 >= 1990
  fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
  for ax, name, title in ((axes[0], 'st', r'smooth trend, $\lambda = 100$'),
                          (axes[1], 'tc', 'trend and AR(2) cycle')):
      y_fc, F_fc = fc[name]
      ax.plot(years50[window], y50[window], 'k', linewidth=1)
      for z, alpha in ((1.645, 0.15), (0.674, 0.3)):
          ax.fill_between(years_fc, y_fc - z*np.sqrt(F_fc), y_fc + z*np.sqrt(F_fc),
                          color='b', alpha=alpha, linewidth=0)
      ax.plot(years_fc, y_fc, 'b', linewidth=1.2)
      ax.set_title(title)
      ax.set_xlabel('year')
  fig.tight_layout()
  fig.savefig("kalman-previsions.svg", transparent=True)
#+end_src

#+CAPTION: *Forecasts of the logarithm of GDP per capita (multiplied by $100$) from 2019 to 2028, with 50% and 90% forecast intervals, for the smooth trend model with $\lambda = 100$ (left) and for the trend and cycle model (right), estimated over 1950--2018.*
#+LABEL: fig:previsions
[[file:kalman-previsions.svg]]

At the end of the application, the Hodrick-Prescott filter appears as a
special case of an explicit model, whose likelihood assesses the
assumptions. For French GDP per capita, it rejects the usual value of the
smoothing parameter, prefers a random walk with drift to the smooth trend
model over the long run, and, since 1950, finds an eight-year cycle once
the cycle is modelled explicitly.

* Exercises
:PROPERTIES:
:CUSTOM_ID: exercices
:END:

#+NAME: ex-ar2-mesure
#+BEGIN_exercice
A stationary second-order autoregressive process is observed with
measurement error: $y_t = x_t + \varepsilon_t$, with
\(x_t = \phi_1x_{t-1} + \phi_2x_{t-2} + \eta_t\), where $\varepsilon_t$ and
$\eta_t$ are independent Gaussian white noises, with variances
$\sigma_\varepsilon^2$ and $\sigma_\eta^2$. Write the model in state space
form and explain how to initialise the filter. Simulate $n = 200$
observations with $\phi_1 = 1.2$, $\phi_2 = -0.5$ and
$\sigma_\varepsilon = \sigma_\eta = 1$, then compare the root mean squared
errors of $y_t$, of the filtered state and of the smoothed state as
estimators of $x_t$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
The state is $\alpha_t = (x_t, x_{t-1})'$, and
\begin{equation*}
Z = \begin{pmatrix} 1 & 0 \end{pmatrix},\quad
T = \begin{pmatrix} \phi_1 & \phi_2\\ 1 & 0 \end{pmatrix},\quad
R = \begin{pmatrix} 1\\ 0 \end{pmatrix},\quad
H = \sigma_\varepsilon^2,\quad
Q = \sigma_\eta^2.
\end{equation*}
The roots of the polynomial $1 - 1.2z + 0.5z^2$ have modulus
$\sqrt 2 > 1$: the process is stationary, and the filter is initialised
with its stationary distribution, $a_0 = 0$ and $P_0$ the solution of the
Lyapunov equation.

#+begin_src python :session kalman-en :exports code :results none
  phi_1, phi_2, n_ex = 1.2, -0.5, 200
  x = np.zeros(n_ex + 100)
  for t in range(2, n_ex + 100):
      x[t] = phi_1*x[t-1] + phi_2*x[t-2] + rng.standard_normal()
  x = x[100:]
  y_ex = x + rng.standard_normal(n_ex)
  Z = np.array([[1.0, 0.0]])
  T = np.array([[phi_1, phi_2],
                [1.0, 0.0]])
  R = np.array([[1.0],
                [0.0]])
  H, Q = np.eye(1), np.eye(1)
  P0 = solve_discrete_lyapunov(T, R @ Q @ R.T)
  res_ex = kalman_filter(y_ex, Z, T, R, H, Q, np.zeros(2), P0)
  a_ex, _ = kalman_smoother(res_ex, T)
  for estimator in (y_ex, res_ex['a_filt'][:, 0], a_ex[:, 0]):
      print(np.sqrt(np.mean((estimator - x)**2)))
#+end_src

The root mean squared errors are $1.03$ for the observation, $0.78$ for the
filtered state and $0.72$ for the smoothed state. Since the stationary
variance of $x_t$ is $3.7$, the signal dominates the noise, and filtering
gains less than in the local level model of the
[[#illustration][illustration]] section. The first $100$ simulated values
are discarded so that $x_t$ is drawn from its stationary distribution.
#+END_details

#+NAME: ex-arima022
#+BEGIN_exercice
Show that in the smooth trend model, the second difference $\Delta^2y_t$ is
an MA(2) process, and compute its autocorrelations as functions of
\(q = \sigma_\zeta^2/\sigma_\varepsilon^2\). Compare them, for
$\lambda = 100$ and for $\hat\lambda$, with the sample autocorrelations of
the second difference of the logarithm of French GDP per capita. What do
you conclude?
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
Since $\Delta^2\mu_t = \zeta_{t-1}$,
\begin{equation*}
\Delta^2y_t = \zeta_{t-1} + \varepsilon_t - 2\varepsilon_{t-1} + \varepsilon_{t-2}.
\end{equation*}
The autocovariances are \(\gamma_0 = \sigma_\zeta^2 + 6\sigma_\varepsilon^2\),
\(\gamma_1 = -4\sigma_\varepsilon^2\), \(\gamma_2 = \sigma_\varepsilon^2\)
and $\gamma_k = 0$ for $k \geq 3$: $\Delta^2 y_t$ is an MA(2), and $y_t$ an
ARIMA(0,2,2), with autocorrelations
\begin{equation*}
\rho_1 = -\frac{4}{6 + q},\qquad \rho_2 = \frac{1}{6 + q},\qquad \rho_k = 0 \text{ for } k \geq 3.
\end{equation*}

#+begin_src python :session kalman-en :exports code :results none
  print(acf(np.diff(y, 2), 4))
  for l in (100, lam_hat):
      print(l, -4/(6 + 1/l), 1/(6 + 1/l))
#+end_src

The model implies $\rho_1 = -0.67$ and $\rho_2 = 0.17$ for $\lambda = 100$,
$\rho_1 = -0.59$ and $\rho_2 = 0.15$ for $\hat\lambda$. The sample
autocorrelations are $-0.41$, $-0.22$, $0.24$ and $0.01$ at the first four
lags. The negative sign of $\hat\rho_2$ and the positive value of
$\hat\rho_3$, both outside the band $\pm 1.96/\sqrt{197} \approx \pm 0.14$,
are incompatible with the model whatever $q$, since it imposes
$\rho_2 > 0$ and $\rho_3 = 0$: the smooth trend model is misspecified for
this series, as confirmed by the superiority of the random walk with drift.
#+END_details

#+NAME: ex-h-nul
#+BEGIN_exercice
Consider the state space representation of an AR($p$) process given in the
examples, with $H = 0$, and an arbitrary initialisation. Show that for
$t \geq p$, the filtered state is \(a_{t|t} = (y_t, \dots, y_{t-p+1})'\) and
that $P_{t|t} = 0$. Deduce that, for $t > p$, the innovation is the
forecast error \(y_t - \phi_1y_{t-1} - \dots - \phi_py_{t-p}\), with
variance $\sigma^2$, and that maximising the terms of the log-likelihood
corresponding to dates $t > p$ gives the ordinary least squares estimator.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
At each date, the observation $y_t$ equals, without error, the first
component of the state. The updating step therefore sets the variance of
this component to zero: if $F_t > 0$, the gain is
\(K_t = P_{t|t-1}Z'/F_t\) with \(F_t = Z P_{t|t-1} Z'\), and the first row
of $I - K_tZ$ is zero since \(ZK_t = 1\), so that $P_{t|t}$ has a zero first
row and a zero first column. The prediction step shifts the state: component
$j+1$ of $\alpha_{t+1}$ is component $j$ of $\alpha_t$, known exactly. By
induction, after $p$ observations, all the components of the state are past
observations, known exactly: \(a_{p|p} = (y_p, \dots, y_1)'\) and
$P_{p|p} = 0$.

For $t > p$, we then have \(a_{t|t-1} = Ta_{t-1|t-1}\), whose first
component is \(\phi_1y_{t-1} + \dots + \phi_py_{t-p}\), and
\(P_{t|t-1} = RQR'\), whose only non-zero element is the first one,
$\sigma^2$. The innovation is
\(v_t = y_t - \phi_1y_{t-1} - \dots - \phi_py_{t-p}\), with variance
$F_t = \sigma^2$, the gain is the first vector of the canonical basis and,
again, $a_{t|t} = (y_t, \dots, y_{t-p+1})'$ and $P_{t|t} = 0$. The terms of
the log-likelihood for $t > p$ are
\begin{equation*}
-\frac{n-p}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2}\sum_{t=p+1}^n\left(y_t - \phi_1y_{t-1} - \dots - \phi_py_{t-p}\right)^2.
\end{equation*}
For any $\sigma^2$, they are maximised in $(\phi_1, \dots, \phi_p)$ when the
sum of squares is minimal: this is the ordinary least squares estimator of
the regression of $y_t$ on its $p$ lagged values. The exact likelihood,
obtained with the stationary initialisation, adds the first $p$ terms, which
depend on the coefficients through the stationary distribution of the first
observations. This is the only difference between the two estimators,
negligible asymptotically but not in small samples (see the note on the
[[https://stephane-adjemian.fr/en/posts/ols-bias-in-an-ar1-process/][bias of the least squares estimator of an AR(1)]]).
#+END_details

#+NAME: ex-limites-hp
#+BEGIN_exercice
Find the limits of the Hodrick-Prescott trend as $\lambda$ tends to zero
and to infinity, and interpret the second one in the smooth trend model.
Check numerically on French GDP per capita, and explain why convergence is
so slow as $\lambda$ increases.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
As $\lambda$ tends to zero, $I_n + \lambda D'D$ tends to $I_n$ and the trend
to the series. As $\lambda$ tends to infinity, the penalty imposes in the
limit $D\mu = 0$, that is zero second differences: the trend is an affine
sequence $\mu_t = a + bt$, and minimising $\sum(y_t - a - bt)^2$ gives the
least squares line. In the smooth trend model, $\lambda \to \infty$
corresponds to $\sigma_\zeta^2 \to 0$: the slope is constant and the trend
deterministic, with a diffuse initial level and slope, whose estimate is
the least squares one.

#+begin_src python :session kalman-en :exports code :results none
  t_ = np.arange(len(y))
  line = np.polyval(np.polyfit(t_, y, 1), t_)
  for l in (1e2, 1e4, 1e6, 1e8, 1e10):
      print(l, np.max(np.abs(hp_filter(y, l) - line)))
  D = sparse.diags([1.0, -2.0, 1.0], [0, 1, 2], shape=(len(y)-2, len(y))).toarray()
  print(np.linalg.eigvalsh(D.T @ D)[:3])
#+end_src

The maximum distance to the line is $62.1$ points for $\lambda = 100$,
$36.5$ for $10^4$, $27.2$ for $10^6$, $1.10$ for $10^8$ and $0.012$ for
$10^{10}$. To understand this, decompose $y$ on an orthonormal basis of
eigenvectors of $D'D$, with eigenvalues $\nu_k$: the trend multiplies the
component of $y$ on the \(k\)-th vector by $1/(1 + \lambda\nu_k)$. The
kernel of $D$, of dimension two, is the set of affine sequences, whose
components are preserved. The others vanish only when $\lambda\nu_k$ is
large, and the smallest non-zero eigenvalue of $D'D$, which decreases like
$1/n^4$, is only $3.2\times 10^{-7}$ for $n = 199$: the corresponding
component, a very slow oscillation, is still multiplied by $0.76$ for
$\lambda = 10^6$, and $\lambda$ must be well above $3\times 10^6$ for it to
disappear.
#+END_details

#+NAME: ex-lisseur-regime-permanent
#+BEGIN_exercice
In the local level model in steady state, show that the smoothing gain is
$\bar J = -\theta$, where $\theta$ is the coefficient of the reduced form,
that the smoother reads
\(a_{t|n} = (1 + \theta)\,a_{t|t} - \theta\,a_{t+1|n}\), and that the
smoothed variance is
\(\bar P_{t|n} = \sigma_\varepsilon^2(1 + \theta)/(1 - \theta)\). Check this
result on the simulated series of the [[#illustration][illustration]]
section.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
In the steady state, $P_{t|t} = \sigma_\varepsilon^2\bar K$ and
$P_{t+1|t} = \sigma_\varepsilon^2\bar p$, hence, with $T = 1$,
\begin{equation*}
\bar J = \frac{\bar K}{\bar p} = \frac{1}{1 + \bar p} = -\theta,
\end{equation*}
by the properties of the reduced form. Since $a_{t+1|t} = a_{t|t}$, the
smoother recursion reads
\(a_{t|n} = a_{t|t} + \bar J(a_{t+1|n} - a_{t|t}) = (1 + \theta)\,a_{t|t} - \theta\,a_{t+1|n}\):
the smoother applies backwards to the filtered series an exponential
smoothing with the same constant $\bar K = 1 + \theta$. The steady-state
smoothed variance is the fixed point of
\(P = \sigma_\varepsilon^2\bar K + \theta^2(P - \sigma_\varepsilon^2\bar p)\),
that is
\begin{equation*}
\bar P_{t|n} = \sigma_\varepsilon^2\,\frac{\bar K - \theta^2\bar p}{1 - \theta^2}
= \sigma_\varepsilon^2\,\frac{(1 + \theta) + \theta(1 + \theta)}{(1 - \theta)(1 + \theta)}
= \sigma_\varepsilon^2\,\frac{1 + \theta}{1 - \theta},
\end{equation*}
using $\bar p = -(1 + \theta)/\theta$. Smoothing therefore divides the
filtered variance $\sigma_\varepsilon^2(1 + \theta)$ by $1 - \theta$.

#+begin_src python :session kalman-en :exports code :results none
  q = sigma_eta**2/sigma_eps**2
  theta = (np.sqrt(q**2 + 4*q) - q - 2)/2
  print(sigma_eps**2*(1 + theta)/(1 - theta), P_smooth[50, 0, 0])
#+end_src

With $q = 0.25$, $\theta = -0.610$, and the formula gives $0.2425$, the
value computed by the smoother in the middle of the simulated sample,
against $0.3904$ for the filtered variance.
#+END_details