#+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: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+LANGUAGE: fr #+STARTUP: latexpreview #+TITLE: Le filtre de Kalman et le filtre de Hodrick-Prescott #+DATE: Septembre 2026 #+AUTHOR: Stéphane Adjemian #+EMAIL: stephane.adjemian@univ-lemans.fr #+PROPERTY: header-args:python :python /tmp/blog-kalman/bin/python #+BEGIN_QUOTE Le filtre de Kalman calcule récursivement la loi d'un vecteur d'état non observé sachant les observations passées, dans un modèle linéaire et gaussien écrit sous forme état-mesure. Cette note dérive en détail le filtre, le lisseur, qui exploite aussi les observations futures, et la vraisemblance, qui s'obtient comme un sous-produit du filtre. Elle étudie la convergence du filtre vers un régime permanent et, pour le modèle de niveau local, son lien avec le lissage exponentiel. L'application montre que le filtre de Hodrick-Prescott est le lisseur d'un modèle de tendance particulier, ce qui permet d'estimer par le maximum de vraisemblance le paramètre de lissage que l'on fixe d'ordinaire a priori. Les calculs portent sur le PIB par tête français depuis 1820 et sont faits en Python, avec une implémentation du filtre en une trentaine de lignes. Les corrigés des exercices sont donnés dans des blocs dépliables. #+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 :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) def virgule(x): return f'{x}'.replace('.', '{,}') #+end_src * Modèles état-mesure :PROPERTIES: :CUSTOM_ID: modele-etat-mesure :END: Un modèle état-mesure distingue ce que l'on observe, une série $y_t$, de la dynamique non observée qui l'engendre, résumée par un vecteur d'état $\alpha_t$. Le cadre est très général\nbsp{}: les processus ARMA, les modèles à composantes inobservées (tendance, cycle, saisonnalité), les régressions à coefficients variables ou les modèles DSGE linéarisés s'y écrivent tous. On reprend la notation de Harvey (/Forecasting, Structural Time Series Models and the Kalman Filter/, 1989) et de Durbin et Koopman (/Time Series Analysis by State Space Methods/, 2012), à une différence de datation près dans l'équation d'état. #+BEGIN_definition Un modèle état-mesure est défini par une équation de mesure et une équation d'état\nbsp{}: \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*} pour $t = 1, \dots, n$, avec la condition initiale $\alpha_0 \sim \mathcal N(a_0, P_0)$. L'observation $y_t$ appartient à $\mathbb R^p$, l'état $\alpha_t$ à $\mathbb R^m$ et la perturbation d'état $\eta_t$ à $\mathbb R^r$. Les matrices $Z_t$ de taille $p\times m$, $T_t$ de taille $m\times m$, $R_t$ de taille $m\times r$, $H_t$ de taille $p\times p$ et $Q_t$ de taille $r\times r$, ainsi que les vecteurs $d_t$ et $c_t$, sont non aléatoires. Les variables $\alpha_0$, $\varepsilon_1, \dots, \varepsilon_n$ et $\eta_1, \dots, \eta_n$ sont mutuellement indépendantes. #+END_definition On note $Y_t = \{y_1, \dots, y_t\}$ l'information disponible à la date $t$, avec $Y_0 = \emptyset$, et $n$ la taille de l'échantillon, la lettre $T$ étant réservée à la matrice de transition. Le vecteur d'état résume tout le passé\nbsp{}: la suite des états est une chaîne de Markov et, sachant les états, les observations sont indépendantes, $y_t$ ne dépendant que de $\alpha_t$. La matrice de sélection $R_t$ permet d'avoir moins de chocs que d'états, ce qui est la règle dès que l'état contient des valeurs retardées. ** Exemples :PROPERTIES: :CUSTOM_ID: exemples :END: #+BEGIN_exemple Le modèle de niveau local est une marche aléatoire observée avec un bruit\nbsp{}: \begin{equation*} y_t = \mu_t + \varepsilon_t,\qquad \mu_t = \mu_{t-1} + \eta_t, \end{equation*} avec $\varepsilon_t \sim \mathcal N(0, \sigma_\varepsilon^2)$ et $\eta_t \sim \mathcal N(0, \sigma_\eta^2)$. L'état est le niveau, $\alpha_t = \mu_t$, et les matrices du système sont des scalaires\nbsp{}: $Z = T = R = 1$, $H = \sigma_\varepsilon^2$, $Q = \sigma_\eta^2$ et $d = c = 0$. #+END_exemple #+BEGIN_exemple Le modèle de tendance linéaire locale ajoute une pente stochastique $\beta_t$ au niveau\nbsp{}: \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*} où les trois perturbations sont des bruits blancs gaussiens indépendants, de variances $\sigma_\varepsilon^2$, $\sigma_\eta^2$ et $\sigma_\zeta^2$. L'état est $\alpha_t = (\mu_t, \beta_t)'$ et \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*} Le cas particulier $\sigma_\eta^2 = 0$ est le /modèle de tendance lisse/\nbsp{}: le niveau n'est plus perturbé directement, seule la pente l'est, et la tendance est une courbe régulière. Seule la seconde colonne de $R$ subsiste, $R = (0, 1)'$, et $Q = \sigma_\zeta^2$. Ce modèle joue le premier rôle dans l'application. #+END_exemple #+BEGIN_exemple Un processus autorégressif d'ordre $p$, $y_t = \phi_1 y_{t-1} + \dots + \phi_p y_{t-p} + \eta_t$, s'écrit sous forme état-mesure en empilant les valeurs retardées dans l'état, $\alpha_t = (y_t, y_{t-1}, \dots, y_{t-p+1})'$\nbsp{}: \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*} avec $H = 0$ et $Q = \sigma^2$. L'observation est parfaite, sans erreur de mesure\nbsp{}: toute l'incertitude porte sur l'état. #+END_exemple #+BEGIN_exemple Pour un processus ARMA($p$, $q$), $y_t = \sum_{i=1}^p \phi_i y_{t-i} + \eta_t + \sum_{j=1}^q \theta_j\eta_{t-j}$, Harvey (1989) pose $s = \max(p, q+1)$, complète les coefficients par des zéros, $\phi_i = 0$ pour $i > p$ et $\theta_j = 0$ pour $j > q$, et prend un état de dimension $s$ avec \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$ et $Q = \sigma^2$. Pour un ARMA(1,1), l'état $\alpha_t = (\alpha_{1,t}, \alpha_{2,t})'$ vérifie $\alpha_{2,t} = \theta\eta_t$ et \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*} et comme $y_t = \alpha_{1,t}$, on retrouve bien le processus ARMA(1,1). Le cas $p = 0$ donne la représentation d'un processus MA($q$). #+END_exemple #+BEGIN_remarque La représentation état-mesure d'un processus n'est pas unique. Pour toute matrice inversible $M$ de taille $m\times m$, l'état $\alpha_t^* = M\alpha_t$ définit un modèle aux mêmes observations, dont les matrices sont \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*} Le filtre donne $a_{t|t}^* = Ma_{t|t}$ et $P_{t|t}^* = MP_{t|t}M'$, mais les innovations, leurs variances, la vraisemblance et les prévisions de $y_t$ sont inchangées. Pour le processus MA(1) $y_t = \eta_t + \theta\eta_{t-1}$, la forme de Harvey ci-dessus, d'état $(y_t, \theta\eta_t)'$, et celle de Hamilton (/Time Series Analysis/, 1994), d'état $(\eta_t, \eta_{t-1})'$, avec $Z = \begin{pmatrix} 1 & \theta \end{pmatrix}$, $T = \begin{pmatrix} 0 & 0\\ 1 & 0 \end{pmatrix}$ et $R = (1, 0)'$, se déduisent l'une de l'autre par \begin{equation*} M = \begin{pmatrix} 1 & \theta\\ \theta & 0 \end{pmatrix}, \end{equation*} inversible dès que $\theta \neq 0$. On peut aussi augmenter la dimension de l'état par des composantes que les chocs n'atteignent pas ou que les observations ne voient pas. Une représentation est minimale lorsqu'elle est commandable et observable, et deux représentations minimales de la même relation entre les chocs et les observations se déduisent l'une de l'autre par un changement de base (Hannan et Deistler, /The Statistical Theory of Linear Systems/, 1988, chapitre 2). Enfin, des structures de chocs différentes peuvent engendrer la même loi des observations\nbsp{}: le modèle de niveau local, à deux chocs indépendants, a la loi d'un processus ARIMA(0,1,1) à un seul choc, comme on le verra dans la section sur la [[#forme-reduite][forme réduite]]. Le choix de la représentation est donc affaire de commodité pour le filtrage et la vraisemblance. En revanche, les paramètres d'une représentation dont toutes les matrices seraient libres ne seraient pas identifiés\nbsp{}: les modèles structurels lèvent cette indétermination en fixant $Z$, $T$ et $R$. #+END_remarque Une composante saisonnière stochastique, un cycle ou des variables explicatives s'ajoutent de la même façon, en augmentant le vecteur d'état\nbsp{}: c'est ce qui fait l'intérêt pratique du cadre. Dans la suite, les matrices du système sont constantes et il n'y a pas de termes déterministes ($d = c = 0$). Le code du filtre traite des observations vectorielles, mais les fonctions de vraisemblance et les applications supposent une observation scalaire ($p = 1$). Un modèle est représenté par une fonction qui renvoie ses cinq matrices\nbsp{}: #+begin_src python :session kalman :exports code :results none def niveau_local(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 tendance_lineaire_locale(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 tendance_lisse(s2_eps, s2_zeta): Z, T, _, H, _ = tendance_lineaire_locale(s2_eps, 0.0, s2_zeta) R = np.array([[0.0], [1.0]]) return Z, T, R, H, np.array([[s2_zeta]]) #+end_src * Le filtre de Kalman :PROPERTIES: :CUSTOM_ID: filtre :END: ** Objectif :PROPERTIES: :CUSTOM_ID: objectif :END: On cherche la loi de l'état $\alpha_t$ sachant une information $Y_s$. On note \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*} son espérance et sa matrice de variance conditionnelles. Trois cas se présentent\nbsp{}: la /prédiction/ si $s = t-1$, le /filtrage/ si $s = t$ et le /lissage/ si $s = n$. Le filtre de Kalman calcule $a_{t|t-1}$, $P_{t|t-1}$, $a_{t|t}$ et $P_{t|t}$ pour $t = 1, \dots, n$, récursivement, en intégrant une observation à chaque étape. Dans un modèle gaussien, la loi conditionnelle de l'état est normale, donc entièrement décrite par ces deux moments, et l'espérance conditionnelle est le meilleur prédicteur de $\alpha_t$ au sens de l'erreur quadratique moyenne, parmi toutes les fonctions des observations. #+BEGIN_remarque Sans l'hypothèse de normalité, mais avec des perturbations non corrélées entre elles et d'une date à l'autre, les récurrences établies plus bas restent valides si l'on remplace l'espérance conditionnelle par la projection linéaire\nbsp{}: $a_{t|t}$ est alors le meilleur prédicteur /linéaire/ de $\alpha_t$ en fonction de $y_1, \dots, y_t$, et $P_{t|t}$ la matrice de variance de son erreur (Harvey, 1989, chapitre 3). #+END_remarque ** Un lemme sur la loi normale :PROPERTIES: :CUSTOM_ID: lemme-gaussien :END: Toute la dérivation repose sur un résultat élémentaire, la loi conditionnelle d'un vecteur gaussien. #+BEGIN_lemma Soit un vecteur gaussien \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*} où $\Sigma_{yy}$ est inversible. La loi de $x$ sachant $y$ est normale\nbsp{}: \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 Posons $e = x - \mu_x - \Sigma_{xy}\Sigma_{yy}^{-1}(y - \mu_y)$. Le vecteur $(e, y)$ est une transformation affine de $(x, y)$, il est donc gaussien. La covariance entre $e$ et $y$ est nulle\nbsp{}: \begin{equation*} \mathrm{Cov}(e, y) = \Sigma_{xy} - \Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yy} = 0. \end{equation*} Pour un vecteur gaussien, l'absence de corrélation équivaut à l'indépendance\nbsp{}: $e$ est indépendant de $y$. Dans la décomposition \begin{equation*} x = \mu_x + \Sigma_{xy}\Sigma_{yy}^{-1}(y - \mu_y) + e, \end{equation*} le premier terme est une fonction de $y$ et le second est indépendant de $y$. Sachant $y$, le vecteur $x$ est donc la somme d'une constante et d'un vecteur gaussien de loi $\mathcal N(0, \mathbb V[e])$, ce qui établit la normalité de la loi conditionnelle et donne son espérance. Sa variance est \begin{equation*} \mathbb V[e] = \mathrm{Cov}(e, x) = \Sigma_{xx} - \Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}, \end{equation*} la première égalité venant de ce que $e$ diffère de $x$ par une fonction affine de $y$, avec laquelle $e$ n'est pas corrélé. #+END_proof ** Les équations du filtre :PROPERTIES: :CUSTOM_ID: equations-filtre :END: #+BEGIN_proposition Partant de $a_{0|0} = a_0$ et $P_{0|0} = P_0$, les moments conditionnels de l'état se calculent pour $t = 1, \dots, n$ par une étape de /prédiction/, \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*} suivie d'une étape de /mise à jour/, \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*} La loi de $\alpha_t$ sachant $Y_t$ est normale, d'espérance $a_{t|t}$ et de variance $P_{t|t}$. Le vecteur $v_t$ est l'/innovation/, $F_t$ sa matrice de variance et $K_t$ le /gain/ du filtre. #+END_proposition #+BEGIN_proof On raisonne par récurrence. Supposons que la loi de $\alpha_{t-1}$ sachant $Y_{t-1}$ soit $\mathcal N(a_{t-1|t-1}, P_{t-1|t-1})$, ce qui est vrai pour $t = 1$ puisque $Y_0$ est vide. /Prédiction./ Les observations $y_1, \dots, y_{t-1}$ sont des fonctions de $\alpha_0$ et des perturbations $\varepsilon_s$ et $\eta_s$ pour $s \leq t-1$\nbsp{}: comme $\alpha_{t-1}$, elles sont indépendantes de $\eta_t$. Sachant $Y_{t-1}$, l'état $\alpha_t = T_t\alpha_{t-1} + c_t + R_t\eta_t$ est donc une transformation affine de deux vecteurs gaussiens indépendants. Sa loi est normale, d'espérance $T_t a_{t-1|t-1} + c_t$ et de variance $T_t P_{t-1|t-1} T_t' + R_t Q_t R_t'$. /Loi jointe de l'état et de l'observation./ La perturbation $\varepsilon_t$ est indépendante de $\alpha_t$ et de $Y_{t-1}$. Sachant $Y_{t-1}$, le couple formé par $\alpha_t$ et $y_t = Z_t\alpha_t + d_t + \varepsilon_t$ est donc gaussien\nbsp{}: \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*} la covariance conditionnelle venant de $\mathrm{Cov}(\alpha_t, Z_t\alpha_t + \varepsilon_t \mid Y_{t-1}) = P_{t|t-1}Z_t'$. /Mise à jour./ Conditionner par $Y_t$, c'est conditionner par $Y_{t-1}$ puis par $y_t$. On applique donc le lemme à la loi jointe précédente, avec $x = \alpha_t$ et $y = y_t$. La loi de $\alpha_t$ sachant $Y_t$ est normale, d'espérance \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*} et de 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}$, ce qui achève la récurrence. #+END_proof #+BEGIN_remarque Le gain dose la correction apportée par l'observation. Dans le cas scalaire, $K_t = P_{t|t-1}Z_t/(Z_t^2P_{t|t-1} + H_t)$. Si la mesure est parfaite ($H_t = 0$), le gain vaut $1/Z_t$ et l'état filtré se déduit exactement de l'observation. Si la mesure est très bruitée ($H_t \to \infty$), ou si l'état est déjà connu avec précision ($P_{t|t-1} \to 0$), le gain tend vers zéro et l'observation est ignorée. On remarque aussi que les variances $P_{t|t-1}$, $F_t$ et $P_{t|t}$ ne dépendent pas des observations\nbsp{}: on pourrait les calculer avant même de disposer des données. #+END_remarque #+BEGIN_property La variance filtrée s'écrit aussi sous la /forme de Joseph/\nbsp{}: \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 En développant le membre de droite, on obtient \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*} Le dernier terme vaut $K_t F_t K_t' = P_{t|t-1}Z_t'K_t'$, puisque $K_t F_t = P_{t|t-1}Z_t'$. Il compense le troisième, et il reste $(I_m - K_tZ_t)P_{t|t-1}$. #+END_proof Les deux expressions de $P_{t|t}$ sont égales en arithmétique exacte, mais pas en arithmétique flottante. La forme de Joseph est une somme de deux matrices symétriques semi-définies positives\nbsp{}: elle le reste quelles que soient les erreurs d'arrondi sur le gain, alors que la forme courte retranche deux termes presque égaux lorsque $P_{t|t-1}$ est grand et peut produire une matrice non symétrique, voire des variances négatives. On verra que ce cas se présente à l'initialisation. L'implémentation suit les équations à la lettre, avec la forme de Joseph. Une observation manquante, codée =np.nan=, est traitée en sautant l'étape de mise à jour\nbsp{}: sans observation, la loi filtrée est la loi prédite. #+begin_src python :session kalman :exports code :results none def filtre_kalman(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): # Prédiction a = T @ a P = T @ P @ T.T + RQR res['a_pred'][t], res['P_pred'][t] = a, P # Mise à jour, sautée si l'observation est manquante 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 # forme de Joseph 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 La fonction renvoie un dictionnaire de tableaux indexés par la date, de $0$ à $n-1$ en Python pour les dates $1$ à $n$ de la note. ** Les innovations :PROPERTIES: :CUSTOM_ID: innovations :END: #+BEGIN_property L'innovation $v_t$ est indépendante de $Y_{t-1}$ et suit la loi $\mathcal N(0, F_t)$. Les innovations $v_1, \dots, v_n$ sont indépendantes\nbsp{}: elles forment un bruit blanc gaussien, hétéroscédastique si $F_t$ varie avec $t$. #+END_property #+BEGIN_proof D'après la loi jointe établie dans la preuve du filtre, la loi de $y_t$ sachant $Y_{t-1}$ est $\mathcal N(Z_t a_{t|t-1} + d_t, F_t)$, donc celle de $v_t$ sachant $Y_{t-1}$ est $\mathcal N(0, F_t)$. Cette loi conditionnelle ne dépend pas de $Y_{t-1}$, puisque $F_t$ n'est pas aléatoire. L'innovation $v_t$ est donc indépendante de $Y_{t-1}$, et sa loi est $\mathcal N(0, F_t)$. Les innovations passées $v_1, \dots, v_{t-1}$ sont des fonctions de $Y_{t-1}$, dont $v_t$ est indépendante. En appliquant cet argument pour $t = n, n-1, \dots, 2$, la densité jointe de $(v_1, \dots, v_n)$ est le produit des densités marginales. #+END_proof L'innovation $v_t$ est la partie de $y_t$ qui n'était pas prévisible à la date $t-1$. Comme $v_t$ est égale à $y_t$ moins une fonction de $Y_{t-1}$, la suite des innovations jusqu'à la date $t$ et celle des observations se déduisent l'une de l'autre\nbsp{}: elles portent la même information. Si le modèle est bien spécifié, les /innovations standardisées/ $e_t = F_t^{-1/2}v_t$ sont indépendantes et de loi normale centrée réduite. Leurs autocorrélations empiriques fournissent donc un diagnostic simple de la spécification, dont l'application fera un large usage. ** Initialisation :PROPERTIES: :CUSTOM_ID: initialisation :END: Il reste à choisir $a_0$ et $P_0$. Si le modèle est stationnaire, avec des matrices constantes et des valeurs propres de $T$ de module strictement inférieur à un, on tire $\alpha_0$ dans la loi stationnaire de l'état\nbsp{}: $a_0 = (I_m - T)^{-1}c$ et $P_0$ est la solution de l'équation de Lyapunov \begin{equation*} P_0 = T P_0 T' + R Q R',\qquad\text{soit}\qquad \mathrm{vec}(P_0) = (I_{m^2} - T\otimes T)^{-1}\,\mathrm{vec}(RQR'), \end{equation*} que calcule la fonction =solve_discrete_lyapunov= de SciPy. Si l'état a des composantes non stationnaires, comme le niveau d'une marche aléatoire, cette loi stationnaire n'existe pas. On traduit alors l'ignorance sur l'état initial par une loi très dispersée, $a_0 = 0$ et $P_0 = \kappa I_m$ avec $\kappa$ grand\nbsp{}: c'est l'/initialisation diffuse approchée/. Dans le modèle de niveau local, la première étape du filtre donne $P_{1|0} = \kappa + \sigma_\eta^2$, puis \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*} Lorsque $\kappa$ tend vers l'infini, $a_{1|1}$ tend vers $y_1$ et $P_{1|1}$ vers $\sigma_\varepsilon^2$\nbsp{}: la première observation fixe le niveau, avec la précision de la mesure, et la loi initiale est oubliée. Plus généralement, si $d$ composantes de l'état sont diffuses, les $d$ premières observations servent à les déterminer, après quoi l'effet de $\kappa$ est d'ordre $1/\kappa$. Trois précautions s'imposent. - L'approximation est d'autant meilleure que $\kappa$ est grand devant la variance des données et devant le carré de leur écart à $a_0$. Avec des données de l'ordre de $1\,000$, comme dans l'application, $\kappa$ doit être bien supérieur à $10^6$. - La variance $P_{t|t}$ des premières étapes est la différence de deux termes d'ordre $\kappa$, et le calcul perd environ $\log_{10}\kappa$ chiffres significatifs sur les seize de l'arithmétique flottante. La forme de Joseph préserve la symétrie et la positivité, mais une valeur excessive de $\kappa$ dégrade la précision. On retient $\kappa = 10^7$ pour les simulations et $\kappa = 10^9$ pour l'application, après vérification. - Les variances $F_t$ des $d$ premières innovations sont d'ordre $\kappa$\nbsp{}: ces innovations ne portent aucune information sur les paramètres, et leur contribution à la vraisemblance dépend de la valeur arbitraire de $\kappa$. Il faut les écarter, comme on le verra dans la [[#concentree][section sur l'estimation]]. #+BEGIN_remarque L'/initialisation diffuse exacte/ traite la limite $\kappa \to \infty$ analytiquement, en séparant dans $P_{t|t-1}$ un terme proportionnel à $\kappa$ et un terme fini (De Jong, /Annals of Statistics/, 1991\nbsp{}; Koopman, /Journal of the American Statistical Association/, 1997\nbsp{}; Durbin et Koopman, 2012, chapitre 5). Elle supprime l'arbitraire du choix de $\kappa$ au prix de récurrences plus lourdes pendant les premières périodes. La bibliothèque =statsmodels= la met en œuvre\nbsp{}; elle servira de contrôle externe dans l'application. #+END_remarque ** Une illustration :PROPERTIES: :CUSTOM_ID: illustration :END: On simule $n = 100$ observations d'un modèle de niveau local, avec $\sigma_\varepsilon = 1$ et $\sigma_\eta = 0{,}5$, et on applique le filtre avec une initialisation diffuse approchée. #+begin_src python :session kalman :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 = niveau_local(sigma_eps**2, sigma_eta**2) res_sim = filtre_kalman(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 L'état filtré estime le niveau avec une erreur quadratique moyenne de $0{,}61$, contre $1{,}01$ pour l'observation elle-même prise comme estimateur naïf. La figure [[fig:niveau-local-filtre][ci-dessous]] montre que l'état filtré suit le niveau avec un léger retard\nbsp{}: n'utilisant que le passé, il ne peut distinguer une variation durable du niveau d'un bruit de mesure qu'après quelques observations. #+begin_src python :session kalman :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: *Modèle de niveau local simulé ($\sigma_\varepsilon = 1$, $\sigma_\eta = 0{,}5$)\nbsp{}: observations, niveau et état filtré, avec une bande de plus ou moins deux écarts-types conditionnels.* #+LABEL: fig:niveau-local-filtre [[file:kalman-niveau-local-filtre.svg]] * Convergence du filtre :PROPERTIES: :CUSTOM_ID: convergence :END: Lorsque les matrices du système sont constantes, la récurrence sur les variances ne dépend ni des observations ni de la date. En reportant l'étape de mise à jour dans l'étape de prédiction, on obtient une récurrence sur la seule variance de prédiction. #+BEGIN_corollary Avec des matrices constantes, la variance de prédiction vérifie l'/équation de Riccati/ \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*} et l'état prédit vérifie $a_{t+1|t} = T a_{t|t-1} + c + T K_t v_t$. #+END_corollary #+BEGIN_proof Il suffit de reporter $a_{t|t} = a_{t|t-1} + K_t v_t$ et $P_{t|t} = P_{t|t-1} - P_{t|t-1}Z'F_t^{-1}ZP_{t|t-1}$ dans les équations de prédiction de la date $t+1$. #+END_proof La suite des variances de prédiction converge sous des conditions très générales. On admet le résultat suivant (Anderson et Moore, /Optimal Filtering/, 1979, chapitre 4, section 4.4). #+BEGIN_proposition Supposons $H$ définie positive, la paire $(T, Z)$ détectable et la paire $(T, RQ^{1/2})$ stabilisable. Alors, pour toute variance initiale semi-définie positive, $P_{t+1|t}$ converge vers l'unique solution semi-définie positive $\bar P$ de l'équation de Riccati algébrique \begin{equation*} \bar P = T\bar PT' - T\bar PZ'(Z\bar PZ' + H)^{-1}Z\bar PT' + RQR', \end{equation*} et le gain $K_t$ converge vers $\bar K = \bar PZ'\bar F^{-1}$, avec $\bar F = Z\bar PZ' + H$. #+END_proposition La détectabilité exige que les composantes instables de l'état se reflètent dans les observations, la stabilisabilité que les chocs atteignent toutes ces composantes. Une fois le régime permanent atteint, le filtre est un filtre linéaire invariant, $a_{t+1|t} = T(I_m - \bar KZ)a_{t|t-1} + T\bar Ky_t$, et le calcul des variances devient inutile. La matrice $\bar P$ est donnée par la fonction =solve_discrete_are= de SciPy, dont la convention transpose celle du filtre\nbsp{}: =solve_discrete_are(T.T, Z.T, R @ Q @ R.T, H)=. ** Le cas du niveau local :PROPERTIES: :CUSTOM_ID: niveau-local-riccati :END: Dans le modèle de niveau local, tout se calcule explicitement. On note $q = \sigma_\eta^2/\sigma_\varepsilon^2$ le rapport signal sur bruit et $p_t = P_{t|t-1}/\sigma_\varepsilon^2$ la variance de prédiction rapportée à celle du bruit de mesure. #+BEGIN_property La suite $p_t$ vérifie $p_{t+1} = f(p_t)$, avec $f(p) = p/(1+p) + q$. Pour tout $p_1 \geq 0$, elle converge de façon monotone vers l'unique point fixe positif \begin{equation*} \bar p = \frac{q + \sqrt{q^2 + 4q}}{2}, \end{equation*} et le gain $K_t = p_t/(1+p_t)$ converge vers $\bar K = \bar p/(1 + \bar p)$. #+END_property #+BEGIN_proof Avec $Z = T = R = 1$, l'équation de Riccati s'écrit $P_{t+1|t} = P_{t|t-1} - P_{t|t-1}^2/(P_{t|t-1} + \sigma_\varepsilon^2) + \sigma_\eta^2$, soit, en divisant par $\sigma_\varepsilon^2$, $p_{t+1} = p_t - p_t^2/(1 + p_t) + q = p_t/(1+p_t) + q$. Les points fixes annulent le polynôme $p^2 - qp - q$, dont $\bar p$ est la seule racine positive. La fonction $f$ est croissante sur $[0, \infty[$, et la fonction $g(p) = f(p) - p = q - p^2/(1+p)$ est décroissante et s'annule en $\bar p$. Si $p_t < \bar p$, on a donc $p_t < f(p_t) < f(\bar p) = \bar p$\nbsp{}: la suite est croissante et majorée par $\bar p$. Symétriquement, elle est décroissante et minorée par $\bar p$ si $p_t > \bar p$. Dans les deux cas elle converge, et sa limite est un point fixe de $f$, c'est-à-dire $\bar p$. #+END_proof Au voisinage du point fixe, l'écart à $\bar p$ est multiplié à chaque période par $f'(\bar p) = (1 + \bar p)^{-2}$\nbsp{}: la convergence est géométrique, d'autant plus rapide que le rapport signal sur bruit est élevé. La figure [[fig:gain][ci-dessous]] le vérifie pour trois valeurs de $q$, avec l'initialisation diffuse approchée. Le gain part de un, puisque la première observation fixe le niveau, et rejoint sa valeur stationnaire, $0{,}200$, $0{,}390$ et $0{,}618$ pour $q$ égal à $0{,}05$, $0{,}25$ et $1$, à $10^{-4}$ près après respectivement $19$, $10$ et $6$ observations. Pour $q = 1$, $\bar p$ est le nombre d'or et $\bar K$ son inverse. #+begin_src python :session kalman :exports code :results none def point_fixe_niveau_local(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 = filtre_kalman(np.zeros(30), *niveau_local(1.0, q), np.zeros(1), kappa*np.eye(1)) K_bar = point_fixe_niveau_local(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 :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 = filtre_kalman(np.zeros(30), *niveau_local(1.0, q), np.zeros(1), kappa*np.eye(1)) p_bar, K_bar = point_fixe_niveau_local(q) axes[0].plot(dates, r['P_pred'][:, 0, 0], c+'o-', markersize=3, linewidth=1, label=fr'$q = {virgule(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 = {virgule(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 de la variance de prédiction (à gauche, échelle logarithmique) et du gain (à droite) dans le modèle de niveau local, pour trois valeurs du rapport signal sur bruit. Les pointillés indiquent les valeurs stationnaires $\bar p$ et $\bar K$.* #+LABEL: fig:gain [[file:kalman-gain.svg]] ** Forme réduite et lissage exponentiel :PROPERTIES: :CUSTOM_ID: forme-reduite :END: Le régime permanent du filtre a une interprétation remarquable dans le modèle de niveau local, qui passe par la forme réduite de ce modèle. #+BEGIN_property Dans le modèle de niveau local, la différence première $\Delta y_t = \eta_t + \varepsilon_t - \varepsilon_{t-1}$ est un processus MA(1) inversible, $\Delta y_t = \xi_t + \theta\xi_{t-1}$, où $\xi_t$ est un bruit blanc de variance $\sigma_\xi^2 = -\sigma_\varepsilon^2/\theta$ et \begin{equation*} \theta = \frac{\sqrt{q^2 + 4q} - q - 2}{2} \in \left]-1, 0\right[. \end{equation*} Le processus $y_t$ est donc un ARIMA(0,1,1). #+END_property #+BEGIN_proof Les autocovariances de $\Delta y_t$ sont $\gamma_0 = \sigma_\eta^2 + 2\sigma_\varepsilon^2$, $\gamma_1 = -\sigma_\varepsilon^2$ et $\gamma_k = 0$ pour $k \geq 2$\nbsp{}: ce sont celles d'un processus MA(1). Celles de $\xi_t + \theta\xi_{t-1}$ sont $\gamma_0 = (1 + \theta^2)\sigma_\xi^2$ et $\gamma_1 = \theta\sigma_\xi^2$. En égalant les autocorrélations d'ordre un, $\theta/(1 + \theta^2) = -1/(q + 2)$, soit $\theta^2 + (q + 2)\theta + 1 = 0$. Le discriminant, $q^2 + 4q$, est positif, les deux racines sont négatives et leur produit vaut un\nbsp{}: une seule est comprise entre $-1$ et $0$, c'est la racine inversible. Enfin, $\sigma_\xi^2 = \gamma_1/\theta = -\sigma_\varepsilon^2/\theta$. #+END_proof #+BEGIN_property Le gain stationnaire vaut $\bar K = 1 + \theta$, et la variance stationnaire de l'innovation, $\bar F = \sigma_\varepsilon^2(1 + \bar p)$, est égale à $\sigma_\xi^2$. En régime permanent, la prévision suit la récurrence du /lissage exponentiel/, \begin{equation*} a_{t+1|t} = (1 - \bar K)\,a_{t|t-1} + \bar K y_t, \end{equation*} et les innovations du filtre sont celles de la forme réduite\nbsp{}: $\Delta y_t = v_t + \theta v_{t-1}$. #+END_property #+BEGIN_proof Notons $r = \sqrt{q^2 + 4q}$, de sorte que $\bar p = (q + r)/2$ et $1 + \theta = (r - q)/2$. En multipliant numérateur et dénominateur par $q + 2 + r$, et puisque $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*} Il en résulte $1 + \bar p = -1/\theta$, puis $\bar K = \bar p/(1 + \bar p) = 1 + \theta$ et $\bar F = \sigma_\varepsilon^2(1 + \bar p) = -\sigma_\varepsilon^2/\theta = \sigma_\xi^2$. En régime permanent, avec $Z = T = 1$, l'état prédit est $a_{t+1|t} = a_{t|t} = a_{t|t-1} + \bar K(y_t - a_{t|t-1})$, qui est la récurrence annoncée. Enfin, $y_t = a_{t|t-1} + v_t$ et $a_{t|t-1} - a_{t-1|t-2} = \bar K v_{t-1}$ donnent \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 Le filtre de Kalman en régime permanent est donc la méthode de prévision par lissage exponentiel, et ce lissage est la prévision optimale d'un processus ARIMA(0,1,1) (Muth, /Journal of the American Statistical Association/, 1960). Le filtre apporte deux choses de plus\nbsp{}: la constante de lissage se déduit des paramètres du modèle, et le gain varie pendant la phase transitoire pour tenir compte de l'incertitude initiale. On le vérifie sur la série simulée, avec $q = 0{,}25$, en comparant la prédiction du filtre avec un lissage exponentiel de constante $\bar K$ initialisé à la première observation\nbsp{}: #+begin_src python :session kalman :exports code :results none q = sigma_eta**2/sigma_eps**2 p_bar, K_bar = point_fixe_niveau_local(q) theta = (np.sqrt(q**2 + 4*q) - q - 2)/2 print(K_bar, 1 + theta) lissage = np.zeros(n) lissage[1] = y_sim[0] for t in range(1, n-1): lissage[t+1] = (1 - K_bar)*lissage[t] + K_bar*y_sim[t] ecart = np.abs(lissage[1:] - res_sim['a_pred'][1:, 0]) print(ecart[[3, 8, 18]]) #+end_src L'écart entre les deux prédictions, de l'ordre de $2\times 10^{-2}$ à la date $5$, n'est plus que de $3\times 10^{-3}$ à la date $10$ et $2\times 10^{-5}$ à la date $20$\nbsp{}: il décroît géométriquement, au taux $1 - \bar K = -\theta \approx 0{,}61$ par période. Le panneau droit de la figure [[fig:lissage-exponentiel][ci-dessous]] vérifie par ailleurs que les innovations standardisées $e_t$ ne sont pas autocorrélées, comme elles doivent l'être lorsque le modèle est bien spécifié. La fonction d'autocorrélation empirique servira encore dans l'application. #+begin_src python :session kalman :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 :exports none :results none fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].semilogy(np.arange(2, 41), ecart[:39], 'b', linewidth=1) axes[0].set_xlabel(r'$t$') axes[0].set_ylabel('écart entre les deux prédictions') k = np.arange(1, 11) bande = 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*bande, color='k', linestyle='--', linewidth=0.8) axes[1].axhline(0, color='k', linewidth=0.5) axes[1].set_xticks(k) axes[1].set_xlabel('retard') axes[1].set_ylabel(r'autocorrélation de $e_t$') fig.tight_layout() fig.savefig("kalman-lissage-exponentiel.svg", transparent=True) #+end_src #+CAPTION: *À gauche, écart en valeur absolue entre la prédiction du filtre $a_{t|t-1}$ et le lissage exponentiel de constante $\bar K$ (échelle logarithmique). À droite, autocorrélations des innovations standardisées, avec la bande $\pm 1{,}96/\sqrt{n}$.* #+LABEL: fig:lissage-exponentiel [[file:kalman-lissage-exponentiel.svg]] * Le lisseur :PROPERTIES: :CUSTOM_ID: lisseur :END: Le filtre estime l'état à la date $t$ avec l'information disponible à cette date. Lorsque tout l'échantillon est observé, il est naturel d'utiliser aussi les observations postérieures\nbsp{}: c'est l'objet du /lissage/, qui calcule $a_{t|n}$ et $P_{t|n}$ pour $t = 1, \dots, n$. Le filtre est unilatéral, le lisseur bilatéral. Le point de départ est une propriété de Markov\nbsp{}: une fois connu l'état à la date $t+1$, les observations postérieures n'apportent plus rien sur l'état à la date $t$. #+BEGIN_lemma Pour $t < n$, la loi de $\alpha_t$ sachant $\alpha_{t+1}$ et $Y_n$ est égale à la loi de $\alpha_t$ sachant $\alpha_{t+1}$ et $Y_t$. #+END_lemma #+BEGIN_proof Notons $Y_{t+1:n} = (y_{t+1}, \dots, y_n)$ les observations postérieures à la date $t$. D'après les équations du modèle, elles s'écrivent comme une fonction de $\alpha_{t+1}$ et des perturbations $\varepsilon_{t+1}, \dots, \varepsilon_n, \eta_{t+2}, \dots, \eta_n$, que l'on regroupe dans un vecteur $\omega$. Le triplet formé par $\alpha_t$, $\alpha_{t+1}$ et $Y_t$ est une fonction de $\alpha_0$, de $\varepsilon_1, \dots, \varepsilon_t$ et de $\eta_1, \dots, \eta_{t+1}$, dont $\omega$ est indépendant. Sachant $\alpha_t$, $\alpha_{t+1}$ et $Y_t$, la loi de $Y_{t+1:n}$ est donc celle d'une fonction de $\alpha_{t+1}$ et de $\omega$\nbsp{}: elle ne dépend que de $\alpha_{t+1}$. Par la formule de Bayes, en notant $p$ les densités, \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*} où le symbole $\propto$ indique la proportionnalité en tant que fonctions de $\alpha_t$. Le premier facteur ne dépend pas de $\alpha_t$\nbsp{}: la densité de $\alpha_t$ sachant $\alpha_{t+1}$ et $Y_n$ est proportionnelle, donc égale, à celle de $\alpha_t$ sachant $\alpha_{t+1}$ et $Y_t$. #+END_proof #+BEGIN_proposition Partant de $a_{n|n}$ et $P_{n|n}$, calculés par le filtre, les moments lissés s'obtiennent par la récurrence rétrograde, pour $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*} en supposant $P_{t+1|t}$ inversible. La loi de $\alpha_t$ sachant $Y_n$ est normale, d'espérance $a_{t|n}$ et de variance $P_{t|n}$. #+END_proposition Cette récurrence est due à Rauch, Tung et Striebel (/AIAA Journal/, 1965). Elle ne demande qu'une passe arrière sur les résultats du filtre. #+BEGIN_proof /Loi jointe./ Sachant $Y_t$, l'état $\alpha_t$ suit la loi $\mathcal N(a_{t|t}, P_{t|t})$ et $\alpha_{t+1} = T_{t+1}\alpha_t + c_{t+1} + R_{t+1}\eta_{t+1}$, où $\eta_{t+1}$ est indépendant de $\alpha_t$ et de $Y_t$. Sachant $Y_t$, le couple formé par $\alpha_t$ et $\alpha_{t+1}$ est donc gaussien\nbsp{}: \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*} /Conditionnement par l'état suivant./ Le lemme sur la loi normale donne \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*} et, d'après le lemme précédent, c'est aussi la loi de $\alpha_t$ sachant $\alpha_{t+1}$ et $Y_n$. /Espérance lissée./ Par la loi des espérances itérées, \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*} /Variance lissée./ Par la formule de la variance totale, \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*} Enfin, la loi de $\alpha_t$ sachant $Y_n$ est normale parce que $\alpha_t$ et $Y_n$ forment un vecteur gaussien. #+END_proof #+BEGIN_remarque La même formule de la variance totale, appliquée en conditionnant par $Y_t$ plutôt que par $\alpha_{t+1}$, donne $P_{t|t} = P_{t|n} + \mathbb V\left[a_{t|n} \mid Y_t\right]$\nbsp{}: la variance lissée est inférieure à la variance filtrée, au sens des matrices semi-définies positives. Lorsque $P_{t+1|t}$ est singulière, ce qui arrive quand l'état contient des composantes observées sans erreur, on peut remplacer l'inverse par un inverse généralisé ou utiliser les récurrences de De Jong (/Journal of the American Statistical Association/, 1989), qui évitent toute inversion. #+END_remarque #+begin_src python :session kalman :exports code :results none def lisseur_kalman(res, T): a_liss, P_liss = res['a_filt'].copy(), res['P_filt'].copy() for t in range(len(a_liss) - 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_liss[t] = res['a_filt'][t] + J @ (a_liss[t+1] - res['a_pred'][t+1]) P_liss[t] = res['P_filt'][t] + J @ (P_liss[t+1] - res['P_pred'][t+1]) @ J.T return a_liss, P_liss a_liss, P_liss = lisseur_kalman(res_sim, T) print(rmse(a_liss[:, 0])) print(res_sim['P_filt'][50, 0, 0], P_liss[50, 0, 0]) #+end_src Sur la série simulée, le lisseur réduit l'erreur quadratique moyenne de $0{,}61$ à $0{,}48$. Au milieu de l'échantillon, la variance lissée vaut $0{,}243$ contre $0{,}390$ pour la variance filtrée\nbsp{}: la bande de confiance est plus étroite d'un cinquième environ, sauf aux extrémités de l'échantillon, où le lisseur ne dispose pas d'observations de part et d'autre et rejoint le filtre à la date $n$. La figure [[fig:niveau-local-lisseur][ci-dessous]] montre aussi que le lisseur corrige le retard du filtre\nbsp{}: il anticipe les variations du niveau, puisqu'il voit les observations suivantes. #+begin_src python :session kalman :exports none :results none s_s = np.sqrt(P_liss[:, 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_liss[:, 0], 'r', linewidth=1, label=r'$a_{t|n}$') ax.fill_between(t_sim, a_liss[:, 0]-2*s_s, a_liss[:, 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: *Niveau simulé, état filtré et état lissé, avec leurs bandes de plus ou moins deux écarts-types conditionnels.* #+LABEL: fig:niveau-local-lisseur [[file:kalman-niveau-local-lisseur.svg]] * Estimation par le maximum de vraisemblance :PROPERTIES: :CUSTOM_ID: estimation :END: Les matrices du système dépendent en général de paramètres inconnus, regroupés dans un vecteur $\psi$\nbsp{}: les variances du modèle de niveau local, les coefficients et la variance d'un processus ARMA. Le filtre de Kalman donne la vraisemblance de ces paramètres presque sans calcul supplémentaire. ** Décomposition de l'erreur de prévision :PROPERTIES: :CUSTOM_ID: decomposition :END: #+BEGIN_proposition La log-vraisemblance des observations $y_1, \dots, y_n$ s'écrit \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*} où les innovations $v_t$ et leurs variances $F_t$ sont calculées par le filtre avec la valeur $\psi$ des paramètres. #+END_proposition #+BEGIN_proof La densité jointe se factorise en un produit de densités conditionnelles, \begin{equation*} f(y_1, \dots, y_n) = \prod_{t=1}^n f(y_t \mid Y_{t-1}), \end{equation*} avec la convention $f(y_1 \mid Y_0) = f(y_1)$. D'après la preuve du filtre, la loi de $y_t$ sachant $Y_{t-1}$ est normale, d'espérance $Z_t a_{t|t-1} + d_t$ et de variance $F_t$, et sa densité en $y_t$ vaut \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*} Il reste à prendre le logarithme du produit. #+END_proof Le calcul direct de la vraisemblance demanderait d'inverser la matrice de variance de l'échantillon, de taille $np \times np$, pour un coût d'ordre $n^3p^3$. Le filtre en calcule implicitement une factorisation triangulaire, pour un coût proportionnel à $n$. Pour un processus ARMA initialisé par sa loi stationnaire, c'est la vraisemblance /exacte/, qui tient compte de la loi des premières observations, et non la vraisemblance conditionnelle aux premières observations que maximisent les moindres carrés (voir la note sur le [[https://stephane-adjemian.fr/posts/biais-de-l-estimateur-mco-d-un-processus-ar1/][biais de l'estimateur des moindres carrés d'un AR(1)]] et l'[[ex-h-nul][exercice]] consacré au cas autorégressif). Le code prévoit d'écarter les $d$ premières innovations, pour la raison indiquée plus bas\nbsp{}: #+begin_src python :session kalman :exports code :results none def log_vraisemblance(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 ** Vraisemblance concentrée :PROPERTIES: :CUSTOM_ID: concentree :END: Dans les modèles à composantes inobservées, les paramètres sont souvent des variances, et l'une d'elles peut être éliminée analytiquement. On suppose $p = 1$ et que les paramètres se décomposent en un facteur d'échelle $\sigma^2$ et un vecteur $\psi$, avec $H = \sigma^2H_*(\psi)$, $Q = \sigma^2Q_*(\psi)$ et $P_0 = \sigma^2P_{0*}(\psi)$. On note $v_t^*$ et $F_t^*$ les quantités calculées par le filtre avec $\sigma^2 = 1$, et $N$ le nombre d'innovations prises en compte dans la vraisemblance. #+BEGIN_property Les innovations ne dépendent pas de $\sigma^2$, $v_t = v_t^*$, et $F_t = \sigma^2F_t^*$. La log-vraisemblance est maximale en $\sigma^2$ pour \begin{equation*} \hat\sigma^2(\psi) = \frac1N\sum_t \frac{(v_t^*)^2}{F_t^*}, \end{equation*} et la /log-vraisemblance concentrée/, obtenue en remplaçant $\sigma^2$ par $\hat\sigma^2(\psi)$, vaut \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 Par récurrence, si $P_{t-1|t-1} = \sigma^2P_{t-1|t-1}^*$, les équations du filtre donnent $P_{t|t-1} = \sigma^2P_{t|t-1}^*$, $F_t = \sigma^2F_t^*$ et $K_t = P_{t|t-1}Z'F_t^{-1} = K_t^*$\nbsp{}: le gain, et donc les états prédits et les innovations, ne dépendent pas de $\sigma^2$, tandis que $P_{t|t} = \sigma^2P_{t|t}^*$. La log-vraisemblance s'écrit alors \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*} Vue comme une fonction de $1/\sigma^2$, elle est concave, et sa dérivée s'annule en $\hat\sigma^2(\psi)$. En reportant, le dernier terme vaut $-N/2$. #+END_proof Avec l'initialisation diffuse approchée, $P_0 = \kappa I_m$ n'est pas proportionnelle à $\sigma^2$, mais l'écart est d'ordre $1/\kappa$. Surtout, les $d$ premières innovations, dont la variance est d'ordre $\kappa$, contribuent à la vraisemblance par des termes $-\frac12\log F_t$ qui dépendent de la valeur arbitraire de $\kappa$\nbsp{}: on les écarte, et les sommes portent sur les dates $t = d+1, \dots, n$, avec $N = n - d$. C'est la /vraisemblance diffuse/ (Harvey, 1989, section 3.4.3). Elle diffère de la vraisemblance calculée avec l'initialisation diffuse exacte par une constante qui ne dépend pas des paramètres, pourvu que la partie diffuse du modèle, dans $Z$ et $T$, n'en dépende pas non plus, comme dans les modèles de tendance de cette note (Francke, Koopman et De Vos, /Journal of Time Series Analysis/, 2010)\nbsp{}: les deux donnent alors les mêmes estimations, ce que l'on vérifiera dans l'application. #+begin_src python :session kalman :exports code :results none def vraisemblance_concentree(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 ** Précision des estimateurs et tests :PROPERTIES: :CUSTOM_ID: erreurs-standard :END: Sous des conditions de régularité, qui supposent notamment que la vraie valeur des paramètres est intérieure à l'espace des paramètres, l'estimateur du maximum de vraisemblance est asymptotiquement normal, et sa matrice de variance est estimée par l'inverse de l'opposé de la hessienne de la log-vraisemblance au maximum. La hessienne se calcule par différences finies centrées\nbsp{}: \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*} où $e_i$ est le $i$-ème vecteur de la base canonique. #+begin_src python :session kalman :exports code :results none def hessienne(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 Les paramètres de variance sont estimés en logarithme, ce qui supprime la contrainte de positivité. Les intervalles de confiance sont construits pour le logarithme puis transformés, ce qui les rend asymétriques. La courbure de la log-vraisemblance concentrée donne directement la précision de $\hat\psi$\nbsp{}: l'inverse de l'opposé de sa hessienne est égal au bloc correspondant de l'inverse de l'opposé de la hessienne de la log-vraisemblance complète (Patefield, /Biometrika/, 1985). Enfin, deux modèles emboîtés se comparent par le test du rapport de vraisemblance, $\mathrm{LR} = 2(\log L_1 - \log L_0)$, dont la loi asymptotique est celle d'un $\chi^2$ à autant de degrés de liberté que de restrictions, pourvu que le modèle restreint ne place pas un paramètre au bord de l'espace des paramètres, comme une variance nulle. * Prévision :PROPERTIES: :CUSTOM_ID: prevision :END: #+BEGIN_proposition Pour $h \geq 1$, la loi de $\alpha_{n+h}$ sachant $Y_n$ est $\mathcal N(a_{n+h|n}, P_{n+h|n})$, où \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*} à partir de $a_{n|n}$ et $P_{n|n}$. La prévision de l'observation est $\hat y_{n+h|n} = Z a_{n+h|n} + d$, de variance $F_{n+h|n} = Z P_{n+h|n} Z' + H$, et, dans le cas scalaire, un intervalle de prévision de niveau $1 - \delta$ est $\hat y_{n+h|n} \pm z_{1-\delta/2}\sqrt{F_{n+h|n}}$, où $z_{1-\delta/2}$ est un quantile de la loi normale centrée réduite. #+END_proposition #+BEGIN_proof Aucune observation n'est disponible après la date $n$\nbsp{}: la loi de $\alpha_{n+h}$ sachant $Y_n$ s'obtient en appliquant $h$ fois l'étape de prédiction du filtre, sans mise à jour, et la loi de $y_{n+h}$ s'en déduit par l'équation de mesure. #+END_proof Autrement dit, prévoir revient à filtrer une série complétée par des observations manquantes aux dates $n+1, \dots, n+h$. Avec des matrices constantes, on obtient $a_{n+h|n} = T^h a_{n|n} + \sum_{j=0}^{h-1}T^jc$ et \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 Dans le modèle de tendance lisse, $T^h = \begin{pmatrix} 1 & h\\ 0 & 1 \end{pmatrix}$ et $ZT^jR = j$. La prévision prolonge linéairement la tendance filtrée, $\hat y_{n+h|n} = a_{\mu, n|n} + h\,a_{\beta, n|n}$, où $a_{\mu, n|n}$ et $a_{\beta, n|n}$ sont les composantes de $a_{n|n}$. Comme $\mu_{n+h} = \mu_n + h\beta_n + \sum_{j=1}^{h-1} j\,\zeta_{n+h-j}$, sa variance est \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*} Le dernier terme, égal à $\sigma_\zeta^2(h-1)h(2h-1)/6$, croît comme $\sigma_\zeta^2h^3/3$\nbsp{}: les intervalles de prévision d'une tendance intégrée d'ordre deux s'élargissent comme $h^{3/2}$. #+END_exemple #+begin_src python :session kalman :exports code :results none def previsions(a, P, Z, T, R, H, Q, h): RQR = R @ Q @ R.T y_prev, F_prev = np.zeros(h), np.zeros(h) for j in range(h): a = T @ a P = T @ P @ T.T + RQR y_prev[j], F_prev[j] = (Z @ a)[0], (Z @ P @ Z.T + H)[0, 0] return y_prev, F_prev #+end_src * Le filtre de Hodrick-Prescott, un lisseur de Kalman :PROPERTIES: :CUSTOM_ID: application :END: Le filtre de Hodrick et Prescott (/Journal of Money, Credit and Banking/, 1997, d'après un document de travail de 1980) est l'outil le plus répandu pour séparer la tendance et le cycle d'une série macroéconomique. On montre qu'il coïncide avec le lisseur de Kalman d'un modèle de tendance lisse, ce qui permet d'estimer son paramètre par le maximum de vraisemblance, puis on applique le résultat au PIB par tête français. ** Le filtre de Hodrick-Prescott :PROPERTIES: :CUSTOM_ID: filtre-hp :END: #+BEGIN_definition Pour une série $y_1, \dots, y_n$ et un paramètre $\lambda > 0$, la tendance de Hodrick-Prescott est la suite $\hat\mu_1, \dots, \hat\mu_n$ qui minimise \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*} Le cycle est l'écart $y_t - \hat\mu_t$. #+END_definition Le critère arbitre entre l'ajustement aux données et la régularité de la tendance, mesurée par ses différences secondes. Pour $\lambda = 0$, la tendance est la série elle-même\nbsp{}; quand $\lambda$ tend vers l'infini, elle tend vers une droite (voir l'[[ex-limites-hp][exercice]]). #+BEGIN_property La tendance de Hodrick-Prescott est l'unique solution du système linéaire $(I_n + \lambda D'D)\hat\mu = y$, où $D$ est la matrice de taille $(n-2)\times n$ dont la ligne $t$ contient les coefficients $1$, $-2$ et $1$ dans les colonnes $t$, $t+1$ et $t+2$, et des zéros ailleurs. #+END_property #+BEGIN_proof Le critère s'écrit $S(\mu) = (y - \mu)'(y - \mu) + \lambda\mu'D'D\mu$. C'est une fonction quadratique de hessienne $2(I_n + \lambda D'D)$, définie positive puisque $D'D$ est semi-définie positive\nbsp{}: le critère est strictement convexe, et son unique minimum annule le gradient $-2(y - \mu) + 2\lambda D'D\mu$. #+END_proof La matrice $I_n + \lambda D'D$ est pentadiagonale\nbsp{}: le système se résout avec un solveur pour matrices creuses, en un nombre d'opérations proportionnel à $n$. #+begin_src python :session kalman :exports code :results none def filtre_hp(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 Appliqué à une série infinie, le filtre est une moyenne mobile symétrique dont le gain vaut $1/(1 + 16\lambda\sin^4(\omega/2))$ pour la tendance, à la fréquence $\omega$ (King et Rebelo, /Journal of Economic Dynamics and Control/, 1993). Ce gain vaut un demi pour les cycles de période $2\pi/\omega$ telle que $16\lambda\sin^4(\omega/2) = 1$\nbsp{}: les fluctuations plus courtes sont attribuées au cycle, les plus longues à la tendance. La valeur usuelle $\lambda = 1\,600$ en données trimestrielles place cette période à une quarantaine de trimestres, soit une dizaine d'années. En données annuelles, $\lambda = 100$, retenu par Backus et Kehoe (/American Economic Review/, 1992), la place à une vingtaine d'années, et $\lambda = 6{,}25$, préconisé par Ravn et Uhlig (/Review of Economics and Statistics/, 2002) pour retrouver le filtre trimestriel usuel, à une dizaine d'années. #+END_remarque ** Le modèle de tendance lisse :PROPERTIES: :CUSTOM_ID: tendance-lisse :END: #+BEGIN_theorem Considérons le modèle de tendance lisse \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*} avec $\varepsilon_t \sim \mathcal N(0, \sigma_\varepsilon^2)$, $\zeta_t \sim \mathcal N(0, \sigma_\zeta^2)$ et une loi a priori diffuse sur le niveau et la pente initiaux. La tendance lissée $\mathbb E[\mu_t \mid Y_n]$ est la tendance de Hodrick-Prescott de paramètre $\lambda = \sigma_\varepsilon^2/\sigma_\zeta^2$. #+END_theorem #+BEGIN_proof Les équations du modèle donnent, pour $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*} La tendance $(\mu_1, \dots, \mu_n)$ est donc l'image de $(\mu_1, \beta_1, \zeta_2, \dots, \zeta_{n-1})$ par l'application linéaire définie par $\mu_2 = \mu_1 + \beta_1$ et $\mu_t = 2\mu_{t-1} - \mu_{t-2} + \zeta_{t-1}$ pour $t \geq 3$. Cette application est bijective, de matrice triangulaire à diagonale unité, et son jacobien vaut un. Avec une densité a priori constante pour $(\mu_1, \beta_1)$, limite d'une loi normale dont la variance tend vers l'infini, et des $\zeta_t$ indépendants de loi $\mathcal N(0, \sigma_\zeta^2)$, la densité a priori de la tendance est \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*} Sachant la tendance, les observations sont indépendantes, de loi $\mathcal N(\mu_t, \sigma_\varepsilon^2)$. Par la formule de Bayes, la densité a posteriori de la tendance est proportionnelle au produit de la vraisemblance et de la densité a priori\nbsp{}: \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*} C'est la densité d'une loi normale. Son espérance est égale à son mode, qui minimise le crochet, c'est-à-dire à la tendance de Hodrick-Prescott. Or le lisseur de Kalman, avec l'initialisation diffuse, calcule précisément l'espérance a posteriori $\mathbb E[\mu_t \mid Y_n]$. La pente $\beta_n$, qui n'intervient pas dans les observations, ne joue aucun rôle. #+END_proof Ce résultat prolonge l'interprétation bayésienne des splines de lissage (Wahba, /Journal of the Royal Statistical Society B/, 1978). Harvey et Jaeger (/Journal of Applied Econometrics/, 1993) en ont tiré les conséquences pour l'analyse du cycle\nbsp{}: la tendance de Hodrick-Prescott ne dépend que du rapport des deux variances, et le filtre est optimal si la série est engendrée par le modèle de tendance lisse avec $\sigma_\varepsilon^2/\sigma_\zeta^2 = \lambda$. Le modèle apporte en retour ce que le filtre ne donne pas\nbsp{}: des bandes de confiance autour de la tendance, une tendance unilatérale, calculée par le filtre, qui n'utilise que le passé, des prévisions et surtout une vraisemblance, qui permet d'estimer $\lambda$ au lieu de le fixer. #+BEGIN_remarque La tendance lissée n'est pas une propriété de la loi des observations, mais de la décomposition retenue. Le modèle de tendance lisse suppose indépendants les chocs de la tendance et ceux de l'écart à la tendance\nbsp{}; d'autres hypothèses sont compatibles avec la même loi de $y_t$ et donnent d'autres tendances. À l'opposé, la décomposition de Beveridge et Nelson (/Journal of Monetary Economics/, 1981) définit la tendance comme la prévision de long terme de la série, corrigée de la croissance tendancielle, ce qui revient à supposer parfaitement corrélés les chocs de la tendance et ceux du cycle\nbsp{}: la tendance est alors une marche aléatoire. Morley, Nelson et Zivot (/Review of Economics and Statistics/, 2003) rappellent que, sur le PIB américain, cette tendance absorbe l'essentiel des fluctuations, et montrent que le modèle de tendance et de cycle de Clark, estimé sans imposer l'indépendance des chocs, donne des estimations filtrées de la tendance et du cycle identiques à celles de Beveridge et Nelson. Sur le PIB trimestriel américain de 1947 à 1998, la corrélation estimée entre les chocs vaut $-0{,}9$, et l'hypothèse d'indépendance est rejetée au seuil de 5\nbsp{}%, de justesse. Pour une marche aléatoire avec dérive, qui s'impose [[#tendance-lineaire-locale][plus bas]] pour le PIB par tête français sur longue période, la tendance de Beveridge et Nelson est la série elle-même et le cycle est nul. La tendance de Hodrick-Prescott n'est donc qu'une décomposition parmi d'autres, fixée par deux hypothèses\nbsp{}: l'indépendance des composantes et la forme de la tendance, dont $\lambda$ règle la régularité. #+END_remarque ** Les données :PROPERTIES: :CUSTOM_ID: donnees :END: On utilise le PIB par tête de la France de la base du projet Maddison, dans sa version 2020, téléchargée comme dans la note sur la [[https://stephane-adjemian.fr/posts/representer-graphiquement-une-serie-temporelle-du-pib/][représentation graphique d'une série temporelle du PIB]]. La série annuelle est continue de 1820 à 2018, soit $n = 199$ observations. On travaille sur $y_t = 100\log(\text{PIB par tête})$, de sorte que les différences premières s'interprètent comme des taux de croissance en pourcentage et les écarts-types en points de pourcentage. #+begin_src python :session kalman :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)] annees = france.year.values.astype(int) y = 100*np.log(france.gdppc.values) dy = np.diff(y) for debut, fin in ((1820, 2018), (1820, 1913), (1950, 1973), (1974, 2018)): print(debut, fin, (y[annees == fin] - y[annees == debut])/(fin - debut)) print(annees[1:][np.argsort(dy)[:2]], np.sort(dy)[:2], annees[1:][np.argmax(dy)], dy.max()) #+end_src La croissance annuelle moyenne est de $1{,}54$ point sur l'ensemble de la période\nbsp{}: $1{,}21$ point par an entre 1820 et 1913, $3{,}94$ entre 1950 et 1973 et $1{,}39$ depuis 1974. Les deux guerres mondiales laissent des chutes considérables, de $21{,}8$ points en 1918 et de $20{,}0$ points en 1941, suivies de rattrapages de grande ampleur, jusqu'à $40{,}4$ points en 1946. Ces années exceptionnelles pèsent sur tous les résultats qui suivent. ** Vérification numérique :PROPERTIES: :CUSTOM_ID: verification :END: On compare la tendance de Hodrick-Prescott, avec $\lambda = 100$, et la tendance lissée du modèle avec $\sigma_\varepsilon^2 = 1$ et $\sigma_\zeta^2 = 1/100$, pour plusieurs valeurs de $\kappa$\nbsp{}: #+begin_src python :session kalman :exports code :results none lam = 100 tendance_hp = filtre_hp(y, lam) def lisse_tendance(y, s2_eps, s2_zeta, kappa): Z, T, R, H, Q = tendance_lisse(s2_eps, s2_zeta) res = filtre_kalman(y, Z, T, R, H, Q, np.zeros(2), kappa*np.eye(2)) return res, lisseur_kalman(res, T) for kap in (1e3, 1e5, 1e7, 1e9, 1e11): _, (a_l, _) = lisse_tendance(y, 1.0, 1/lam, kap) print(kap, np.max(np.abs(tendance_hp - a_l[:, 0]))) kappa = 1e9 #+end_src | $\kappa$ | Écart maximal | |-----------+-----------------------| | $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}$ | L'écart décroît comme $1/\kappa$ tant que l'approximation diffuse domine, puis remonte lorsque les erreurs d'arrondi l'emportent\nbsp{}: avec $\kappa = 10^{11}$, le calcul des premières variances filtrées perd une dizaine de chiffres significatifs. Il est maximal en début d'échantillon, là où l'initialisation compte, comme le montre la figure [[fig:pib-france][ci-dessous]]. On retient désormais $\kappa = 10^9$. #+begin_src python :session kalman :exports none :results none res_hp, (a_hp, P_hp) = lisse_tendance(y, 1.0, 1/lam, kappa) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(annees, y, 'k', linewidth=0.8, label=r'$y_t$') axes[0].plot(annees, tendance_hp, 'b', linewidth=1.5, label='Hodrick-Prescott') axes[0].plot(annees, a_hp[:, 0], 'r--', linewidth=1.5, label='lisseur de Kalman') axes[0].set_xlabel('année') axes[0].legend() axes[1].plot(annees, 1e7*(tendance_hp - a_hp[:, 0]), 'b', linewidth=1) axes[1].axhline(0, color='k', linewidth=0.5) axes[1].set_xlabel('année') axes[1].set_ylabel(r'écart ($\times 10^{-7}$)') fig.tight_layout() fig.savefig("kalman-pib-france.svg", transparent=True) #+end_src #+CAPTION: *À gauche, logarithme du PIB par tête de la France (multiplié par $100$), tendance de Hodrick-Prescott avec $\lambda = 100$ et tendance lissée du modèle de tendance lisse avec $\kappa = 10^9$. À droite, écart entre les deux tendances.* #+LABEL: fig:pib-france [[file:kalman-pib-france.svg]] Le contrôle externe utilise la classe =UnobservedComponents= de =statsmodels=, qui met en œuvre l'initialisation diffuse exacte. On compare les tendances lissées et les log-vraisemblances pour trois couples de variances\nbsp{}: #+begin_src python :session kalman :exports code :results none from statsmodels.tsa.statespace.structural import UnobservedComponents modele_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_l, _) = lisse_tendance(y, s2e, s2z, kappa) tendance_sm = modele_sm.smooth([s2e, s2z]).smoothed_state[0] print(np.max(np.abs(tendance_sm - a_l[:, 0])), modele_sm.loglike([s2e, s2z]) - log_vraisemblance(res, d=2)) #+end_src Les tendances diffèrent de moins de $1{,}4\times 10^{-5}$, et la différence des log-vraisemblances vaut $-1{,}83787$ à $2\times 10^{-5}$ près pour les trois couples, soit $-\log 2\pi$\nbsp{}: comme annoncé, la vraisemblance diffuse et la vraisemblance diffuse exacte ne diffèrent que par une constante, qui n'affecte ni les estimations ni les tests. ** Estimer le paramètre de lissage :PROPERTIES: :CUSTOM_ID: estimer-lambda :END: Le paramètre $\lambda$ est l'inverse du rapport signal sur bruit $q = \sigma_\zeta^2/\sigma_\varepsilon^2$. On concentre la vraisemblance par rapport à $\sigma_\varepsilon^2$ et on la maximise en $\log q$, en écartant les $d = 2$ premières innovations, puisque le niveau et la pente initiaux sont diffus. On calcule aussi la vraisemblance avec $\lambda = 100$, et l'estimation sur deux sous-périodes. #+begin_src python :session kalman :exports code :results none def logLc_tendance_lisse(log_q, y): Z, T, R, H, Q = tendance_lisse(1.0, np.exp(log_q)) res = filtre_kalman(y, Z, T, R, H, Q, np.zeros(2), kappa*np.eye(2)) return vraisemblance_concentree(res, d=2) def estime_lambda(y): opt = minimize_scalar(lambda x: -logLc_tendance_lisse(x, y)[0], bounds=(np.log(1e-6), np.log(1e3)), method='bounded', options={'xatol': 1e-8}) return opt.x log_q = estime_lambda(y) logL_max, s2_eps = logLc_tendance_lisse(log_q, y) lam_hat = np.exp(-log_q) se = 1/np.sqrt(-hessienne(lambda x: logLc_tendance_lisse(x[0], y)[0], [log_q])[0, 0]) ic = np.exp(-log_q + np.array([-1.96, 1.96])*se) logL_100, s2_eps_100 = logLc_tendance_lisse(np.log(1/100), y) LR = 2*(logL_max - logL_100) print(lam_hat, se, ic, LR, chi2.sf(LR, 1)) print(2*np.pi/(2*np.arcsin((16*lam_hat)**-0.25))) for debut in (1870, 1950): print(debut, np.exp(-estime_lambda(y[annees >= debut]))) #+end_src L'estimation donne $\hat\lambda = 1{,}24$, avec un intervalle de confiance à 95\nbsp{}% de $[0{,}65\,;\,2{,}36]$, construit à partir de l'écart-type de $\log\hat q$, qui vaut $0{,}33$. On est très loin des valeurs usuelles\nbsp{}: le test du rapport de vraisemblance rejette $\lambda = 100$ avec une statistique de $93{,}8$, pour une valeur critique de $3{,}84$ au seuil de 5\nbsp{}%. D'après la remarque sur la réponse en fréquence, le gain de la tendance vaut un demi pour une période de $6{,}4$ ans avec $\hat\lambda$, contre une vingtaine d'années avec $\lambda = 100$\nbsp{}: la tendance estimée absorbe l'essentiel des fluctuations de moyen terme. Le résultat n'est pas dû aux seules guerres\nbsp{}: l'estimation vaut $0{,}86$ à partir de 1870 et $0{,}53$ à partir de 1950. La figure [[fig:vraisemblance][ci-dessous]] montre la log-vraisemblance concentrée en fonction de $\lambda$\nbsp{}: $\lambda = 6{,}25$ est déjà nettement moins vraisemblable que $\hat\lambda$, et $\lambda = 100$ ou $\lambda = 1\,600$ le sont beaucoup moins encore. Pour comprendre ce rejet, on compare les deux modèles, avec les variances estimées, à l'aide des innovations standardisées et du cycle $y_t - a_{\mu, t|n}$\nbsp{}: #+begin_src python :session kalman :exports code :results none def ajuste(y, s2_eps, lam): res, (a_l, P_l) = lisse_tendance(y, s2_eps, s2_eps/lam, kappa) e = res['v'][2:, 0]/np.sqrt(res['F'][2:, 0, 0]) return res, a_l, P_l, e, y - a_l[:, 0] resultats = {'100': ajuste(y, s2_eps_100, 100), 'hat': ajuste(y, s2_eps, lam_hat)} for nom, (res, a_l, P_l, e, cycle) in resultats.items(): print(nom, 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-vraisemblance | $-704{,}4$ | $-657{,}5$ | | Autocorrélation des innovations (1) | $0{,}66$ | $0{,}08$ | | Autocorrélation des innovations (2) | $0{,}34$ | $-0{,}26$ | | Écart-type du cycle | $5{,}94$ | $2{,}70$ | | Autocorrélation du cycle (1) | $0{,}57$ | $-0{,}22$ | #+begin_src python :session kalman :exports none :results none lam_grille = np.logspace(-2, 4, 121) profil = [logLc_tendance_lisse(-np.log(l), y)[0] for l in lam_grille] fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(lam_grille, profil, 'b', linewidth=1) for l, c, etiquette 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=etiquette) axes[0].set_xscale('log') axes[0].set_xlabel(r'$\lambda$') axes[0].set_ylabel('log-vraisemblance concentrée') axes[0].legend() k = np.arange(1, 11) axes[1].bar(k - 0.17, acf(resultats['100'][3], 10), width=0.34, color='k', label=r'$\lambda = 100$') axes[1].bar(k + 0.17, acf(resultats['hat'][3], 10), width=0.34, color='r', label=r'$\hat\lambda$') bande = 1.96/np.sqrt(len(resultats['hat'][3])) for s in (-1, 1): axes[1].axhline(s*bande, color='k', linestyle='--', linewidth=0.8) axes[1].axhline(0, color='k', linewidth=0.5) axes[1].set_xticks(k) axes[1].set_xlabel('retard') axes[1].set_ylabel('autocorrélation des innovations') axes[1].legend() fig.tight_layout() fig.savefig("kalman-vraisemblance.svg", transparent=True) #+end_src #+CAPTION: *À gauche, log-vraisemblance concentrée du modèle de tendance lisse en fonction de $\lambda$ (échelle logarithmique). À droite, autocorrélations des innovations standardisées pour $\lambda = 100$ et pour $\hat\lambda$, avec la bande $\pm 1{,}96/\sqrt{n-2}$.* #+LABEL: fig:vraisemblance [[file:kalman-vraisemblance.svg]] Les deux valeurs de $\lambda$ racontent des histoires très différentes (figure [[fig:cycles][ci-dessous]]). Avec $\lambda = 100$, la tendance est régulière et le cycle, d'écart-type $5{,}9$ points, est très persistant\nbsp{}: il enregistre les guerres comme des récessions de plusieurs années. Mais le modèle suppose que l'écart à la tendance, $\varepsilon_t$, est un bruit blanc, et ses innovations standardisées ont une autocorrélation de $0{,}66$ à l'ordre un\nbsp{}: la spécification est rejetée par ses propres diagnostics. Le maximum de vraisemblance lève la contradiction en rendant la tendance flexible. Avec $\hat\lambda$, la tendance absorbe les fluctuations persistantes et le cycle résiduel, d'écart-type $2{,}7$ points, n'a plus de persistance. Le modèle de tendance lisse n'a pas de place pour un cycle persistant\nbsp{}: il ne connaît que la tendance et le bruit. Fixer $\lambda$, c'est donc choisir une définition du cycle, pas l'estimer, conclusion à laquelle arrivaient déjà Harvey et Jaeger (1993). #+begin_src python :session kalman :exports none :results none fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(annees, y, 'k', linewidth=0.8, label=r'$y_t$') axes[0].plot(annees, resultats['100'][1][:, 0], 'b', linewidth=1.2, label=r'tendance, $\lambda = 100$') axes[0].plot(annees, resultats['hat'][1][:, 0], 'r', linewidth=1.2, label=r'tendance, $\hat\lambda$') axes[0].set_xlabel('année') axes[0].legend() axes[1].plot(annees, resultats['100'][4], 'b', linewidth=1, label=r'$\lambda = 100$') axes[1].plot(annees, resultats['hat'][4], 'r', linewidth=1, label=r'$\hat\lambda$') axes[1].axhline(0, color='k', linewidth=0.5) axes[1].set_xlabel('année') 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: *Tendances lissées (à gauche) et cycles (à droite) pour $\lambda = 100$ et pour $\hat\lambda$.* #+LABEL: fig:cycles [[file:kalman-cycles.svg]] *** Libérer la variance du niveau :PROPERTIES: :CUSTOM_ID: tendance-lineaire-locale :END: Le modèle de tendance lisse impose $\sigma_\eta^2 = 0$. Le modèle de tendance linéaire locale, qui estime les trois variances, tranche autrement. On concentre la vraisemblance par rapport à $\sigma_\eta^2$, les rapports des deux autres variances à $\sigma_\eta^2$ étant estimés en logarithme et bornés inférieurement par $e^{-20}$. #+begin_src python :session kalman :exports code :results none def logLc_tll(x, y): Z, T, R, H, Q = tendance_lineaire_locale(np.exp(x[0]), 1.0, np.exp(x[1])) res = filtre_kalman(y, Z, T, R, H, Q, np.zeros(2), kappa*np.eye(2)) return vraisemblance_concentree(res, d=2) opt_tll = minimize(lambda x: -logLc_tll(x, y)[0], [0.0, -2.0], method='L-BFGS-B', bounds=[(-20.0, 5.0)]*2) logL_tll, s2_eta = logLc_tll(opt_tll.x, y) print(opt_tll.x, s2_eta, logL_tll, 2*(logL_tll - logL_max)) #+end_src L'optimum est atteint au bord\nbsp{}: les deux rapports sont à leur borne inférieure ou presque, et $\hat\sigma_\eta = 5{,}99$. Le modèle se réduit à une marche aléatoire avec dérive, $\Delta y_t = \beta + \eta_t$, dont la pente constante est la croissance moyenne, $1{,}54$ point par an, et dont les chocs ont l'écart-type de la croissance annuelle. La log-vraisemblance, $-635{,}0$, dépasse de $22{,}6$ celle du modèle de tendance lisse. La statistique du rapport de vraisemblance, $45{,}1$, ne suit pas une loi du $\chi^2$, puisque l'hypothèse $\sigma_\eta^2 = 0$ place un paramètre au bord de l'espace des paramètres. Dans le cadre d'observations indépendantes étudié par Self et Liang (/Journal of the American Statistical Association/, 1987), sa loi asymptotique serait un mélange à parts égales d'une masse en zéro et d'une loi $\chi^2(1)$\nbsp{}; pour une variance d'état dans un modèle de séries temporelles, elle peut être différente. Avec une statistique de cette ampleur, la conclusion ne fait cependant aucun doute. La marche aléatoire avec dérive, qui n'a qu'un paramètre de variance, est bien plus vraisemblable que le modèle de tendance lisse, qui en a deux\nbsp{}: chaque variation annuelle du PIB par tête y est un choc permanent sur le niveau, y compris les chutes et les rattrapages des guerres. #+BEGIN_remarque On retrouve le phénomène d'accumulation à la frontière signalé dans la note sur le [[https://stephane-adjemian.fr/posts/biais-de-l-estimateur-mco-d-un-processus-ar1/#etat-mesure][biais de l'estimateur des moindres carrés d'un AR(1)]]\nbsp{}: dans les modèles à composantes inobservées, les variances estimées sont exactement nulles avec une probabilité positive. Surtout, si la série est une marche aléatoire, le filtre de Hodrick-Prescott fabrique un cycle persistant qu'elle ne contient pas (Cogley et Nason, /Journal of Economic Dynamics and Control/, 1995)\nbsp{}: le cycle obtenu avec $\lambda = 100$ relève en partie de cet artefact. #+END_remarque ** Tendance filtrée et tendance lissée :PROPERTIES: :CUSTOM_ID: filtre-vs-lisseur :END: La tendance filtrée $a_{\mu, t|t}$ n'utilise que les observations passées\nbsp{}: c'est la version unilatérale, en temps réel, du filtre de Hodrick-Prescott (Stock et Watson, /Journal of Monetary Economics/, 1999). Son gain rejoint rapidement sa valeur stationnaire, que l'on calcule par l'équation de Riccati algébrique\nbsp{}: #+begin_src python :session kalman :exports code :results none for nom, l in (('100', 100), ('hat', lam_hat)): Z, T, R, H, Q = tendance_lisse(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 = resultats[nom][0]['K'][:, :, 0] print(nom, K_bar.ravel(), annees[np.argmax(np.all(np.abs(K_t - K_bar.ravel()) < 1e-3, axis=1))]) res, a_l, P_l, e, cycle = resultats['100'] revision = res['a_filt'][:, 0] - a_l[:, 0] print(annees[np.argmax(np.abs(revision))], revision.min()) print(np.sqrt(res['P_filt'][annees == 1930, 0, 0]), np.sqrt(P_l[annees == 1930, 0, 0])) print(acf(y - res['a_filt'][:, 0], 1), acf(y - a_l[:, 0], 1)) #+end_src Avec $\lambda = 100$, le gain stationnaire vaut $(0{,}362\,;\,0{,}080)'$, atteint à $10^{-3}$ près dès 1834\nbsp{}: chaque année, la tendance filtrée corrige son niveau de 36\nbsp{}% de l'innovation et sa pente de 8\nbsp{}%. Avec $\hat\lambda$, il vaut $(0{,}750\,;\,0{,}450)'$, atteint dès 1825. La figure [[fig:filtre-vs-lisseur][ci-dessous]] compare les deux tendances pour $\lambda = 100$ autour des deux guerres. En temps réel, les effondrements sont en partie attribués à la tendance, qui décroche\nbsp{}; avec le recul, le lisseur, qui voit la reconstruction, les attribue au cycle. En 1945, la tendance filtrée est inférieure de $31{,}5$ points à la tendance lissée. L'écart-type conditionnel de la tendance est de $4{,}1$ points pour le filtre et de $2{,}3$ points pour le lisseur, loin des extrémités de l'échantillon. Le cycle unilatéral est plus persistant que le cycle bilatéral, avec une autocorrélation d'ordre un de $0{,}66$ contre $0{,}57$. #+begin_src python :session kalman :exports none :results none res, a_l, P_l, e, cycle = resultats['100'] fenetre = (annees >= 1900) & (annees <= 1960) af, sf = res['a_filt'][:, 0], np.sqrt(res['P_filt'][:, 0, 0]) al, sl = a_l[:, 0], np.sqrt(P_l[:, 0, 0]) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(annees[fenetre], y[fenetre], 'k', linewidth=0.8, label=r'$y_t$') axes[0].plot(annees[fenetre], af[fenetre], 'b', linewidth=1.2, label=r'tendance filtrée $a_{\mu, t|t}$') axes[0].fill_between(annees[fenetre], (af-2*sf)[fenetre], (af+2*sf)[fenetre], color='b', alpha=0.12, linewidth=0) axes[0].plot(annees[fenetre], al[fenetre], 'r', linewidth=1.2, label=r'tendance lissée $a_{\mu, t|n}$') axes[0].fill_between(annees[fenetre], (al-2*sl)[fenetre], (al+2*sl)[fenetre], color='r', alpha=0.2, linewidth=0) axes[0].set_xlabel('année') axes[0].legend(loc='upper left') axes[1].plot(annees, y - af, 'b', linewidth=1, label='cycle unilatéral') axes[1].plot(annees, y - al, 'r', linewidth=1, label='cycle bilatéral') axes[1].axhline(0, color='k', linewidth=0.5) axes[1].set_xlabel('année') axes[1].legend() fig.tight_layout() fig.savefig("kalman-filtre-vs-lisseur.svg", transparent=True) #+end_src #+CAPTION: *Modèle de tendance lisse avec $\lambda = 100$. À gauche, tendances filtrée et lissée de 1900 à 1960, avec des bandes de plus ou moins deux écarts-types conditionnels. À droite, cycles unilatéral, $y_t - a_{\mu, t|t}$, et bilatéral, $y_t - a_{\mu, t|n}$, qui est le cycle de Hodrick-Prescott.* #+LABEL: fig:filtre-vs-lisseur [[file:kalman-filtre-vs-lisseur.svg]] #+BEGIN_remarque Le traitement des observations manquantes permet d'interpoler la tendance. Si l'on retire les années 1939 à 1945, le lisseur prolonge la tendance à travers la guerre\nbsp{}: il l'estime à $880{,}8$ points en 1942, avec un écart-type de $3{,}7$ points, contre $846{,}6$ observés et une tendance de $861{,}4$ avec toutes les données. #+END_remarque #+begin_src python :session kalman :exports code :results none y_manquant = y.copy() y_manquant[(annees >= 1939) & (annees <= 1945)] = np.nan _, (a_m, P_m) = lisse_tendance(y_manquant, s2_eps_100, s2_eps_100/100, kappa) print(a_m[annees == 1942, 0], np.sqrt(P_m[annees == 1942, 0, 0])) print(y[annees == 1942], resultats['100'][1][annees == 1942, 0]) #+end_src ** Une tendance et un cycle :PROPERTIES: :CUSTOM_ID: tendance-cycle :END: Pour obtenir un cycle persistant qui ne soit pas imposé par le choix de $\lambda$, il faut le modéliser. On remplace le bruit $\varepsilon_t$ par un processus autorégressif d'ordre deux stationnaire. C'est une variante du modèle de Clark (/Quarterly Journal of Economics/, 1987), dont la tendance comporte en outre un choc sur le niveau\nbsp{}: \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*} avec $u_t \sim \mathcal N(0, \sigma_u^2)$. L'état est $\alpha_t = (\mu_t, \beta_t, \psi_t, \psi_{t-1})'$, la mesure est parfaite ($H = 0$) et \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*} L'initialisation est mixte\nbsp{}: diffuse pour la tendance, stationnaire, par l'équation de Lyapunov, pour le cycle. On concentre la vraisemblance par rapport à $\sigma_u^2$. Pour garantir la stationnarité du cycle, on paramètre $(\phi_1, \phi_2)$ par les autocorrélations partielles $r_1, r_2 \in \left]-1, 1\right[$, avec $\phi_2 = r_2$ et $\phi_1 = r_1(1 - r_2)$ (Barndorff-Nielsen et Schou, /Journal of Multivariate Analysis/, 1973), elles-mêmes écrites $r_i = \tanh(x_i)$ avec $x_i$ réel. #+begin_src python :session kalman :exports code :results none def tendance_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_tendance_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 parametres(x): r1, r2 = np.tanh(x[1]), np.tanh(x[2]) return np.exp(x[0]), r1*(1 - r2), r2 def logLc_tendance_cycle(x, y): Z, T, R, H, Q = tendance_cycle(*parametres(x), 1.0) res = filtre_kalman(y, Z, T, R, H, Q, *init_tendance_cycle(T, R, Q)) return vraisemblance_concentree(res, d=2) #+end_src Les guerres interdisent d'estimer ce modèle sur l'échantillon complet\nbsp{}: le cycle y absorbe les chutes et les rattrapages, sa racine autorégressive dominante vaut $0{,}95$ et la variance de la pente devient presque nulle, de sorte que le modèle se rapproche d'une marche aléatoire avec dérive. On l'estime donc sur la période 1950--2018, à partir de quatre points de départ, car la vraisemblance des modèles à composantes inobservées a souvent plusieurs maxima locaux. #+begin_src python :session kalman :exports code :results none y50, annees50 = y[annees >= 1950], annees[annees >= 1950] meilleur = 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_tendance_cycle(x, y50)[0], x0, method='Nelder-Mead', options={'xatol': 1e-8, 'fatol': 1e-10, 'maxiter': 5000}) print(x0, -opt.fun) if meilleur is None or opt.fun < meilleur.fun: meilleur = opt logL_tc, s2_u = logLc_tendance_cycle(meilleur.x, y50) q_zeta, phi1, phi2 = parametres(meilleur.x) racines = np.roots([1, -phi1, -phi2]) print(q_zeta*s2_u, phi1, phi2, s2_u, np.abs(racines), 2*np.pi/np.abs(np.angle(racines[0]))) logL_st50 = logLc_tendance_lisse(estime_lambda(y50), y50)[0] print(logL_tc, logL_st50, chi2.sf(2*(logL_tc - logL_st50), 2)) #+end_src L'une des quatre optimisations s'arrête sur un maximum local où la première autocorrélation partielle du cycle tend vers un\nbsp{}; les trois autres convergent vers le même point. L'estimation donne $\hat\phi_1 = 0{,}807$ et $\hat\phi_2 = -0{,}306$\nbsp{}: les racines du polynôme autorégressif sont complexes, de module $0{,}55$, et le cycle a une période de $8{,}3$ ans. L'écart-type des chocs du cycle est de $1{,}03$ point, celui de la pente de la tendance de $0{,}40$ point. Le modèle améliore nettement la log-vraisemblance du modèle de tendance lisse estimé sur la même période, qui passe de $-126{,}4$ à $-119{,}1$\nbsp{}: la statistique du rapport de vraisemblance vaut $14{,}6$ pour les deux restrictions $\phi_1 = \phi_2 = 0$, soit une probabilité critique de $0{,}0007$. #+begin_src python :session kalman :exports code :results none Z_tc, T_tc, R_tc, H_tc, Q_tc = tendance_cycle(q_zeta*s2_u, phi1, phi2, s2_u) res_tc = filtre_kalman(y50, Z_tc, T_tc, R_tc, H_tc, Q_tc, *init_tendance_cycle(T_tc, R_tc, Q_tc)) a_tc, P_tc = lisseur_kalman(res_tc, T_tc) ecart_hp = lambda l: np.sqrt(np.mean((filtre_hp(y50, l) - a_tc[:, 0])**2)) opt = minimize_scalar(lambda x: ecart_hp(np.exp(x)), bounds=(np.log(0.1), np.log(1e4)), method='bounded') print(np.exp(opt.x), ecart_hp(np.exp(opt.x)), ecart_hp(100)) print(acf(res_tc['v'][2:, 0]/np.sqrt(res_tc['F'][2:, 0, 0]), 3)) print(a_tc[-1, 2], (y50 - filtre_hp(y50, 100))[-1]) #+end_src Les innovations standardisées ne présentent plus d'autocorrélation. La tendance du modèle est très proche d'une tendance de Hodrick-Prescott, mais pas de paramètre $100$\nbsp{}: la valeur de $\lambda$ qui la reproduit le mieux est $27$, avec un écart quadratique moyen de $0{,}10$ point, contre $0{,}39$ pour $\lambda = 100$. Le cycle lissé (figure [[fig:tendance-cycle][ci-dessous]]) retrouve les récessions de 1975, 1993 et 2009, et place l'économie légèrement au-dessus de sa tendance en 2018, de $0{,}8$ point contre $1{,}6$ pour le cycle de Hodrick-Prescott. #+begin_src python :session kalman :exports none :results none hp50 = filtre_hp(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(annees50, y50, 'k', linewidth=0.8, label=r'$y_t$') axes[0].plot(annees50, a_tc[:, 0], 'r', linewidth=1.2, label='tendance du modèle') axes[0].plot(annees50, hp50, 'b--', linewidth=1.2, label=r'Hodrick-Prescott, $\lambda = 100$') axes[0].set_xlabel('année') axes[0].legend() axes[1].plot(annees50, c, 'r', linewidth=1.2, label=r'cycle $a_{\psi, t|n}$') axes[1].fill_between(annees50, c-2*sc, c+2*sc, color='r', alpha=0.2, linewidth=0) axes[1].plot(annees50, y50 - hp50, 'b--', linewidth=1, label=r'cycle HP, $\lambda = 100$') axes[1].axhline(0, color='k', linewidth=0.5) axes[1].set_xlabel('année') axes[1].legend() fig.tight_layout() fig.savefig("kalman-tendance-cycle.svg", transparent=True) #+end_src #+CAPTION: *Modèle de tendance et de cycle estimé sur 1950--2018. À gauche, tendance lissée et tendance de Hodrick-Prescott avec $\lambda = 100$. À droite, cycle lissé avec une bande de plus ou moins deux écarts-types conditionnels, et cycle de Hodrick-Prescott.* #+LABEL: fig:tendance-cycle [[file:kalman-tendance-cycle.svg]] ** Prévisions :PROPERTIES: :CUSTOM_ID: previsions-pib :END: On compare pour finir les prévisions de 2019 à 2028 de deux modèles estimés sur 1950--2018\nbsp{}: le modèle de tendance lisse avec $\lambda = 100$, la variance $\sigma_\varepsilon^2$ étant estimée, et le modèle de tendance et de cycle. #+begin_src python :session kalman :exports code :results none h = 10 s2_100_50 = logLc_tendance_lisse(np.log(1/100), y50)[1] modele_st = tendance_lisse(s2_100_50, s2_100_50/100) res_st = filtre_kalman(y50, *modele_st, np.zeros(2), kappa*np.eye(2)) prev = {'st': previsions(res_st['a_filt'][-1], res_st['P_filt'][-1], *modele_st, h), 'tc': previsions(res_tc['a_filt'][-1], res_tc['P_filt'][-1], Z_tc, T_tc, R_tc, H_tc, Q_tc, h)} for nom, (y_prev, F_prev) in prev.items(): print(nom, y_prev[0] - y50[-1], np.sqrt(F_prev[[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 Avec $\lambda = 100$, la prévision pour 2019 est inférieure de $1{,}0$ point au niveau de 2018\nbsp{}: le cycle, supposé sans mémoire, disparaît immédiatement et la prévision repart de la tendance filtrée, qu'elle prolonge avec une pente de $0{,}60$ point par an. Le modèle avec cycle prévoit au contraire une croissance de $0{,}8$ point en 2019, puisque le cycle ne se résorbe que progressivement, puis une pente de $1{,}04$ point. Les écarts-types de prévision sont de $2{,}2$, $3{,}3$ et $5{,}6$ points à un, cinq et dix ans avec $\lambda = 100$, contre $1{,}4$, $4{,}9$ et $10{,}2$ points avec le cycle (figure [[fig:previsions][ci-dessous]]). La précision apparente du premier modèle à long terme tient à la faible variance qu'il impose à la pente de la tendance, $\sigma_\zeta = 0{,}17$ point contre $0{,}40$\nbsp{}: fixer $\lambda$, c'est aussi fixer l'incertitude sur la croissance future. #+begin_src python :session kalman :exports none :results none annees_prev = np.arange(annees50[-1] + 1, annees50[-1] + h + 1) fenetre = annees50 >= 1990 fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True) for ax, nom, titre in ((axes[0], 'st', r'tendance lisse, $\lambda = 100$'), (axes[1], 'tc', 'tendance et cycle AR(2)')): y_prev, F_prev = prev[nom] ax.plot(annees50[fenetre], y50[fenetre], 'k', linewidth=1) for z, alpha in ((1.645, 0.15), (0.674, 0.3)): ax.fill_between(annees_prev, y_prev - z*np.sqrt(F_prev), y_prev + z*np.sqrt(F_prev), color='b', alpha=alpha, linewidth=0) ax.plot(annees_prev, y_prev, 'b', linewidth=1.2) ax.set_title(titre) ax.set_xlabel('année') fig.tight_layout() fig.savefig("kalman-previsions.svg", transparent=True) #+end_src #+CAPTION: *Prévisions du logarithme du PIB par tête (multiplié par $100$) de 2019 à 2028, avec des intervalles de prévision à 50\nbsp{}% et 90\nbsp{}%, pour le modèle de tendance lisse avec $\lambda = 100$ (à gauche) et pour le modèle de tendance et de cycle (à droite), estimés sur 1950--2018.* #+LABEL: fig:previsions [[file:kalman-previsions.svg]] Au terme de l'application, le filtre de Hodrick-Prescott apparaît comme un cas particulier d'un modèle explicite, dont la vraisemblance juge les hypothèses. Sur le PIB par tête français, elle rejette la valeur usuelle du paramètre de lissage, préfère une marche aléatoire avec dérive au modèle de tendance lisse sur longue période, et, depuis 1950, retrouve un cycle de huit ans une fois celui-ci modélisé explicitement. * Exercices :PROPERTIES: :CUSTOM_ID: exercices :END: #+NAME: ex-ar2-mesure #+BEGIN_exercice On observe avec une erreur de mesure un processus autorégressif d'ordre deux stationnaire\nbsp{}: $y_t = x_t + \varepsilon_t$, avec $x_t = \phi_1x_{t-1} + \phi_2x_{t-2} + \eta_t$, où $\varepsilon_t$ et $\eta_t$ sont des bruits blancs gaussiens indépendants, de variances $\sigma_\varepsilon^2$ et $\sigma_\eta^2$. Écrivez le modèle sous forme état-mesure et indiquez comment initialiser le filtre. Simulez $n = 200$ observations avec $\phi_1 = 1{,}2$, $\phi_2 = -0{,}5$ et $\sigma_\varepsilon = \sigma_\eta = 1$, puis comparez les erreurs quadratiques moyennes de $y_t$, de l'état filtré et de l'état lissé comme estimateurs de $x_t$. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary L'état est $\alpha_t = (x_t, x_{t-1})'$, et \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*} Les racines du polynôme $1 - 1{,}2z + 0{,}5z^2$ sont de module $\sqrt 2 > 1$\nbsp{}: le processus est stationnaire, et l'on initialise le filtre par sa loi stationnaire, $a_0 = 0$ et $P_0$ solution de l'équation de Lyapunov. #+begin_src python :session kalman :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 = filtre_kalman(y_ex, Z, T, R, H, Q, np.zeros(2), P0) a_ex, _ = lisseur_kalman(res_ex, T) for estimateur in (y_ex, res_ex['a_filt'][:, 0], a_ex[:, 0]): print(np.sqrt(np.mean((estimateur - x)**2))) #+end_src Les erreurs quadratiques moyennes valent $1{,}03$ pour l'observation, $0{,}78$ pour l'état filtré et $0{,}72$ pour l'état lissé. La variance stationnaire de $x_t$ étant de $3{,}7$, le signal domine le bruit, et le gain à filtrer est moindre que dans le modèle de niveau local de la section [[#illustration][d'illustration]]. Les $100$ premières valeurs simulées sont écartées pour que $x_t$ soit tiré dans sa loi stationnaire. #+END_details #+NAME: ex-arima022 #+BEGIN_exercice Montrez que dans le modèle de tendance lisse, la différence seconde $\Delta^2y_t$ est un processus MA(2), dont vous calculerez les autocorrélations en fonction de $q = \sigma_\zeta^2/\sigma_\varepsilon^2$. Comparez-les, pour $\lambda = 100$ et pour $\hat\lambda$, aux autocorrélations empiriques de la différence seconde du logarithme du PIB par tête français. Qu'en concluez-vous\nbsp{}? #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary Puisque $\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*} Les autocovariances sont $\gamma_0 = \sigma_\zeta^2 + 6\sigma_\varepsilon^2$, $\gamma_1 = -4\sigma_\varepsilon^2$, $\gamma_2 = \sigma_\varepsilon^2$ et $\gamma_k = 0$ pour $k \geq 3$\nbsp{}: $\Delta^2 y_t$ est un MA(2), et $y_t$ un ARIMA(0,2,2), d'autocorrélations \begin{equation*} \rho_1 = -\frac{4}{6 + q},\qquad \rho_2 = \frac{1}{6 + q},\qquad \rho_k = 0 \text{ pour } k \geq 3. \end{equation*} #+begin_src python :session kalman :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 Le modèle implique $\rho_1 = -0{,}67$ et $\rho_2 = 0{,}17$ pour $\lambda = 100$, $\rho_1 = -0{,}59$ et $\rho_2 = 0{,}15$ pour $\hat\lambda$. Les autocorrélations empiriques sont $-0{,}41$, $-0{,}22$, $0{,}24$ et $0{,}01$ aux quatre premiers ordres. Le signe négatif de $\hat\rho_2$ et la valeur positive de $\hat\rho_3$, tous deux hors de la bande $\pm 1{,}96/\sqrt{197} \approx \pm 0{,}14$, sont incompatibles avec le modèle quel que soit $q$, qui impose $\rho_2 > 0$ et $\rho_3 = 0$\nbsp{}: le modèle de tendance lisse est mal spécifié pour cette série, ce que confirme la supériorité de la marche aléatoire avec dérive. #+END_details #+NAME: ex-h-nul #+BEGIN_exercice On considère la représentation état-mesure d'un processus AR($p$) donnée dans les exemples, avec $H = 0$, et une initialisation quelconque. Montrez que pour $t \geq p$, l'état filtré est $a_{t|t} = (y_t, \dots, y_{t-p+1})'$ et que $P_{t|t} = 0$. Déduisez-en que, pour $t > p$, l'innovation est l'erreur de prévision $y_t - \phi_1y_{t-1} - \dots - \phi_py_{t-p}$, de variance $\sigma^2$, et que la maximisation des termes de la log-vraisemblance correspondant aux dates $t > p$ donne l'estimateur des moindres carrés ordinaires. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary À chaque date, l'observation $y_t$ est égale, sans erreur, à la première composante de l'état. La mise à jour annule donc la variance de cette composante\nbsp{}: si $F_t > 0$, le gain vaut $K_t = P_{t|t-1}Z'/F_t$ avec $F_t = Z P_{t|t-1} Z'$, et la première ligne de $I - K_tZ$ est nulle puisque $ZK_t = 1$, de sorte que $P_{t|t}$ a une première ligne et une première colonne nulles. L'étape de prédiction décale l'état\nbsp{}: la composante $j+1$ de $\alpha_{t+1}$ est la composante $j$ de $\alpha_t$, connue exactement. Par récurrence, après $p$ observations, toutes les composantes de l'état sont des observations passées, connues exactement\nbsp{}: $a_{p|p} = (y_p, \dots, y_1)'$ et $P_{p|p} = 0$. Pour $t > p$, on a alors $a_{t|t-1} = Ta_{t-1|t-1}$, dont la première composante est $\phi_1y_{t-1} + \dots + \phi_py_{t-p}$, et $P_{t|t-1} = RQR'$, dont seul le premier élément, $\sigma^2$, est non nul. L'innovation est $v_t = y_t - \phi_1y_{t-1} - \dots - \phi_py_{t-p}$, de variance $F_t = \sigma^2$, le gain est le premier vecteur de la base canonique et, à nouveau, $a_{t|t} = (y_t, \dots, y_{t-p+1})'$ et $P_{t|t} = 0$. Les termes de la log-vraisemblance pour $t > p$ sont \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*} Pour tout $\sigma^2$, ils sont maximaux en $(\phi_1, \dots, \phi_p)$ lorsque la somme des carrés est minimale\nbsp{}: c'est l'estimateur des moindres carrés ordinaires de la régression de $y_t$ sur ses $p$ valeurs retardées. La vraisemblance exacte, obtenue avec l'initialisation stationnaire, ajoute les $p$ premiers termes, qui dépendent des coefficients à travers la loi stationnaire des premières observations. C'est la seule différence entre les deux estimateurs, négligeable asymptotiquement mais pas en petit échantillon (voir la note sur le [[https://stephane-adjemian.fr/posts/biais-de-l-estimateur-mco-d-un-processus-ar1/][biais de l'estimateur des moindres carrés d'un AR(1)]]). #+END_details #+NAME: ex-limites-hp #+BEGIN_exercice Déterminez les limites de la tendance de Hodrick-Prescott quand $\lambda$ tend vers zéro et vers l'infini, et interprétez la seconde dans le modèle de tendance lisse. Vérifiez numériquement sur le PIB par tête français, et expliquez pourquoi la convergence est si lente quand $\lambda$ augmente. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary Quand $\lambda$ tend vers zéro, $I_n + \lambda D'D$ tend vers $I_n$ et la tendance vers la série. Quand $\lambda$ tend vers l'infini, la pénalité impose à la limite $D\mu = 0$, c'est-à-dire des différences secondes nulles\nbsp{}: la tendance est une suite affine $\mu_t = a + bt$, et minimiser $\sum(y_t - a - bt)^2$ donne la droite des moindres carrés. Dans le modèle de tendance lisse, $\lambda \to \infty$ correspond à $\sigma_\zeta^2 \to 0$\nbsp{}: la pente est constante et la tendance déterministe, avec un niveau et une pente initiaux diffus, dont l'estimation est celle des moindres carrés. #+begin_src python :session kalman :exports code :results none t_ = np.arange(len(y)) droite = np.polyval(np.polyfit(t_, y, 1), t_) for l in (1e2, 1e4, 1e6, 1e8, 1e10): print(l, np.max(np.abs(filtre_hp(y, l) - droite))) 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 L'écart maximal à la droite vaut $62{,}1$ points pour $\lambda = 100$, $36{,}5$ pour $10^4$, $27{,}2$ pour $10^6$, $1{,}10$ pour $10^8$ et $0{,}012$ pour $10^{10}$. Pour le comprendre, décomposons $y$ sur une base orthonormée de vecteurs propres de $D'D$, de valeurs propres $\nu_k$\nbsp{}: la tendance multiplie la composante de $y$ sur le $k$-ème vecteur par $1/(1 + \lambda\nu_k)$. Le noyau de $D$, de dimension deux, est l'ensemble des suites affines, dont les composantes sont conservées. Les autres ne disparaissent que lorsque $\lambda\nu_k$ est grand, et la plus petite valeur propre non nulle de $D'D$, qui décroît comme $1/n^4$, ne vaut que $3{,}2\times 10^{-7}$ pour $n = 199$\nbsp{}: la composante correspondante, une oscillation très lente, est encore multipliée par $0{,}76$ pour $\lambda = 10^6$, et il faut $\lambda$ bien supérieur à $3\times 10^6$ pour qu'elle disparaisse. #+END_details #+NAME: ex-lisseur-regime-permanent #+BEGIN_exercice Dans le modèle de niveau local en régime permanent, montrez que le gain de lissage vaut $\bar J = -\theta$, où $\theta$ est le coefficient de la forme réduite, que le lisseur s'écrit $a_{t|n} = (1 + \theta)\,a_{t|t} - \theta\,a_{t+1|n}$, et que la variance lissée vaut $\bar P_{t|n} = \sigma_\varepsilon^2(1 + \theta)/(1 - \theta)$. Vérifiez ce résultat sur la série simulée de la section [[#illustration][d'illustration]]. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary En régime permanent, $P_{t|t} = \sigma_\varepsilon^2\bar K$ et $P_{t+1|t} = \sigma_\varepsilon^2\bar p$, d'où, avec $T = 1$, \begin{equation*} \bar J = \frac{\bar K}{\bar p} = \frac{1}{1 + \bar p} = -\theta, \end{equation*} d'après les propriétés de la forme réduite. Comme $a_{t+1|t} = a_{t|t}$, la récurrence du lisseur s'écrit $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}$\nbsp{}: le lisseur applique à rebours à la série filtrée un lissage exponentiel de même constante $\bar K = 1 + \theta$. La variance lissée stationnaire est le point fixe de $P = \sigma_\varepsilon^2\bar K + \theta^2(P - \sigma_\varepsilon^2\bar p)$, soit \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*} en utilisant $\bar p = -(1 + \theta)/\theta$. Le lissage divise donc la variance filtrée $\sigma_\varepsilon^2(1 + \theta)$ par $1 - \theta$. #+begin_src python :session kalman :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_liss[50, 0, 0]) #+end_src Avec $q = 0{,}25$, $\theta = -0{,}610$, et la formule donne $0{,}2425$, la valeur calculée par le lisseur au milieu de l'échantillon simulé, contre $0{,}3904$ pour la variance filtrée. #+END_details