#+OPTIONS:   H:3 num:nil toc:nil \n:nil @:t ::t |:t ^:nil -:t f:t *:t TeX:t LaTeX:t skip:t d:t tags:not-in-toc creator:t timestamp:nil author:nil title:nil html5-fancy:t
#+HTML_DOCTYPE: html5
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../css/stylesheet-blog.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../css/theorems-fr.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../fontawesome/css/all.css" />
#+HTML_HEAD: <link rel="stylesheet" href="../../highlight/styles/base16/solarized-light.min.css">
#+HTML_HEAD: <script src="../../highlight/highlight.min.js"></script>
#+HTML_HEAD: <script>hljs.highlightAll();</script>
#+LANGUAGE: fr
#+STARTUP: latexpreview
#+TITLE: Biais de l'estimateur des moindres carrés d'un processus AR(1)
#+DATE: Septembre 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-ar1-biais/bin/python

#+BEGIN_QUOTE
L'estimateur des moindres carrés ordinaires du coefficient d'un processus
autorégressif d'ordre un est convergent et asymptotiquement normal, mais il
est biaisé en échantillon fini\nbsp{}: en moyenne, il sous-estime la
persistance du processus. Cette note explique l'origine du biais, qui
tient à ce que le régresseur n'est pas strictement exogène, établit
l'approximation classique $-2\rho/T$ par un développement au premier ordre,
puis la confronte à des simulations de Monte Carlo en Python\nbsp{}: biais
en fonction de la taille de l'échantillon et de la persistance, distribution
de l'estimateur, et deux façons de corriger le biais, analytique et par
bootstrap. 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-ar1-biais
  source /tmp/blog-ar1-biais/bin/activate
  pip install numpy scipy matplotlib
#+END_SRC

#+begin_src python :session ar1 :exports none :results none
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.stats import gaussian_kde, norm
#+end_src

* Le modèle et l'estimateur
:PROPERTIES:
:CUSTOM_ID: modele
:END:

On considère un processus autorégressif d'ordre un\nbsp{}:

\begin{equation*}
y_t = \rho\,y_{t-1} + \varepsilon_t,\qquad t = 1,\dots,T
\end{equation*}

où $|\rho| < 1$ et où les innovations $\varepsilon_t$ sont indépendantes et
identiquement distribuées selon une loi normale centrée de variance
$\sigma^2$. On suppose que le processus démarre en $y_0 = 0$, de sorte que
$y_1 = \varepsilon_1$. Le processus n'est alors pas stationnaire mais
asymptotiquement stationnaire\nbsp{}: la variance de $y_t$,
$\sigma^2(1-\rho^{2t})/(1-\rho^2)$, converge vers la variance stationnaire
$\gamma_0 = \sigma^2/(1-\rho^2)$. Nous verrons dans un
[[ex-stationnaire][exercice]] que ce choix de la condition initiale n'a
pratiquement pas d'incidence sur les résultats.

On dispose d'un échantillon $y_1,\dots,y_T$ et l'on cherche à estimer
$\rho$. L'estimateur des moindres carrés ordinaires (MCO), obtenu en
régressant $y_t$ sur $y_{t-1}$ sans constante, est\nbsp{}:

\begin{equation*}
\hat\rho = \frac{\sum_{t=2}^T y_t\,y_{t-1}}{\sum_{t=2}^T y_{t-1}^2}
\end{equation*}

C'est aussi l'estimateur du maximum de vraisemblance conditionnel à la
première observation, puisque la log-vraisemblance conditionnelle est, à
une constante près, $-\frac{1}{2\sigma^2}\sum_{t=2}^T(y_t - \rho y_{t-1})^2$,
que $\hat\rho$ maximise. En substituant $y_t = \rho y_{t-1} + \varepsilon_t$
au numérateur, on obtient la décomposition qui servira dans toute la
suite\nbsp{}:

\begin{equation*}
\hat\rho - \rho = \frac{\sum_{t=2}^T y_{t-1}\,\varepsilon_t}{\sum_{t=2}^T y_{t-1}^2}
\equiv \frac{N}{D}
\end{equation*}

L'erreur d'estimation est le rapport de deux variables aléatoires, un
numérateur $N$ et un dénominateur $D$, qui dépendent toutes deux des
innovations.

#+NAME: ex-sigma
#+BEGIN_exercice
Montrez que la distribution de $\hat\rho$ ne dépend pas de $\sigma^2$. Le
biais de l'estimateur est donc une fonction de $\rho$ et de $T$ seulement.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Écrivons $\varepsilon_t = \sigma\eta_t$, où les $\eta_t$ sont indépendants
et de loi normale centrée réduite. Par récurrence à partir de $y_0 = 0$,
$y_t = \sigma\sum_{j=0}^{t-1}\rho^j\eta_{t-j} \equiv \sigma z_t$, où le
processus $z_t$ ne dépend que de $\rho$ et des $\eta_t$. En substituant
dans l'estimateur\nbsp{}:
\begin{equation*}
\hat\rho = \frac{\sigma^2\sum_{t=2}^T z_t\,z_{t-1}}{\sigma^2\sum_{t=2}^T z_{t-1}^2} = \frac{\sum_{t=2}^T z_t\,z_{t-1}}{\sum_{t=2}^T z_{t-1}^2}
\end{equation*}
Le facteur $\sigma^2$ se simplifie\nbsp{}: $\hat\rho$ est une fonction de
$\rho$, de $T$ et des $\eta_t$ seulement, et il en va de même de sa loi. Dans
les simulations, on peut donc poser $\sigma = 1$ sans perte de généralité.
#+END_details

* Convergence et loi asymptotique
:PROPERTIES:
:CUSTOM_ID: asymptotique
:END:

#+BEGIN_property
L'estimateur des MCO est convergent et asymptotiquement normal\nbsp{}:
\begin{equation*}
\hat\rho \xrightarrow[T\to\infty]{p} \rho
\qquad\text{et}\qquad
\sqrt T\,(\hat\rho - \rho) \xrightarrow[T\to\infty]{d} \mathcal N\bigl(0, 1-\rho^2\bigr)
\end{equation*}
#+END_property

#+BEGIN_proof
Divisons le numérateur et le dénominateur de $\hat\rho - \rho$ par $T$.
Le processus $y_t$ étant asymptotiquement stationnaire et ergodique, la loi
des grands nombres donne $\frac1T\sum_{t=2}^Ty_{t-1}^2 \to \gamma_0 = \sigma^2/(1-\rho^2)$
en probabilité. Le terme $y_{t-1}\varepsilon_t$ est centré, car
$\varepsilon_t$ est indépendant de $y_{t-1}$, qui ne dépend que des
innovations passées, et \(\mathbb E[\varepsilon_t] = 0\)\nbsp{}; les termes
$y_{t-1}\varepsilon_t$ sont de plus non corrélés entre eux, car pour $s < t$
on a $\mathbb E[y_{s-1}\varepsilon_sy_{t-1}\varepsilon_t] = \mathbb E[y_{s-1}\varepsilon_sy_{t-1}]\,\mathbb E[\varepsilon_t] = 0$.
La loi des grands nombres donne donc $\frac1T\sum_{t=2}^Ty_{t-1}\varepsilon_t \to 0$,
d'où la convergence de $\hat\rho$ vers $\rho$. Pour la loi asymptotique, on
admet que le théorème central limite s'applique à la suite
$y_{t-1}\varepsilon_t$, dont les termes sont centrés, non corrélés et de
variance $\mathbb E[y_{t-1}^2]\,\mathbb E[\varepsilon_t^2] \to \gamma_0\sigma^2$[fn:: La
suite $y_{t-1}\varepsilon_t$ est une différence de martingale, et le
théorème central limite pour les différences de martingale s'applique.]\nbsp{}:
\begin{equation*}
\frac{1}{\sqrt T}\sum_{t=2}^Ty_{t-1}\varepsilon_t \xrightarrow{d} \mathcal N\bigl(0, \gamma_0\sigma^2\bigr)
\end{equation*}
Par le théorème de Slutsky, $\sqrt T(\hat\rho - \rho)$ converge en loi vers
une normale centrée de variance $\gamma_0\sigma^2/\gamma_0^2 = \sigma^2/\gamma_0 = 1-\rho^2$.
#+END_proof

#+BEGIN_remarque
La variance asymptotique $1-\rho^2$ ne dépend pas de $\sigma^2$, ce qui est
cohérent avec l'exercice précédent, et elle diminue quand la persistance
augmente\nbsp{}: plus $\rho$ est proche de un, plus le régresseur $y_{t-1}$
est dispersé relativement à l'innovation, et plus l'estimateur est précis.
En grand échantillon, on approche donc la loi de $\hat\rho$ par
$\mathcal N\bigl(\rho, (1-\rho^2)/T\bigr)$, une loi centrée sur la vraie
valeur. Cette approximation ne dit rien du biais en échantillon fini, qui
est d'ordre $1/T$ et disparaît dans la normalisation par $\sqrt T$.
#+END_remarque

* Un biais en échantillon fini
:PROPERTIES:
:CUSTOM_ID: biais
:END:

Dans le modèle de régression linéaire avec régresseurs déterministes,
l'estimateur des MCO est sans biais parce que l'espérance de l'erreur est
nulle conditionnellement à /tous/ les régresseurs. Ici le régresseur est
aléatoire, et l'on a bien $\mathbb E[\varepsilon_t\mid y_{t-1}] = 0$, mais
pas \(\mathbb E[\varepsilon_t\mid y_1,\dots,y_{T-1}] = 0\)\nbsp{}: l'innovation
$\varepsilon_t$ entre dans $y_t$, $y_{t+1}$, \dots, qui sont les régresseurs
des dates suivantes. Le régresseur est /prédéterminé/, il n'est pas
/strictement exogène/. Pour voir comment cette dépendance produit un
biais, écrivons l'erreur d'estimation comme une somme pondérée des
innovations\nbsp{}:

\begin{equation*}
\hat\rho - \rho = \sum_{t=2}^T w_t\,\varepsilon_t
\qquad\text{avec}\qquad
w_t = \frac{y_{t-1}}{\sum_{s=2}^Ty_{s-1}^2}
\end{equation*}

Si les poids $w_t$ étaient indépendants des innovations, l'espérance de
chaque terme serait nulle. Mais le poids $w_t$ dépend de $\varepsilon_t$
par son dénominateur\nbsp{}: une innovation positive augmente $y_t$,
$y_{t+1}$, \dots, donc le dénominateur, et réduit le poids attribué à cette
innovation. Une innovation positive reçoit ainsi, en moyenne, un poids
plus faible qu'une innovation négative de même amplitude, et la somme
pondérée est négative en espérance lorsque $\rho > 0$. La propriété
suivante quantifie cet effet.

#+BEGIN_property
Au premier ordre en $1/T$, le biais de l'estimateur des MCO est[fn:: Ce
résultat est dû à Marriott et Pope (/Biometrika/, 1954) et à White
(/Biometrika/, 1961)\nbsp{}; le biais de l'estimateur des MCO dans les
modèles dynamiques a été mis en évidence par Hurwicz (1950), dans le volume
de la Cowles Commission édité par Koopmans, et on parle parfois de biais de
Hurwicz.]\nbsp{}:
\begin{equation*}
\mathbb E[\hat\rho] - \rho = -\frac{2\rho}{T} + O\!\left(\frac{1}{T^2}\right)
\end{equation*}
#+END_property

#+BEGIN_proof
Nous donnons la dérivation par un développement au premier ordre, dont
nous admettons la validité[fn:: La justification rigoureuse du
développement, et le calcul des termes d'ordre supérieur, sont dans
l'article de White (1961).]. Notons $\bar D = \mathbb E[D]$ et écrivons
$D = \bar D(1 + \delta)$ avec $\delta = (D-\bar D)/\bar D$, qui est d'ordre
$1/\sqrt T$. Le développement $(1+\delta)^{-1} = 1 - \delta + O(\delta^2)$
donne\nbsp{}:
\begin{equation*}
\begin{split}
\hat\rho - \rho &= \frac{N}{\bar D}(1 - \delta) + \dots\\
\mathbb E[\hat\rho] - \rho &= \frac{\mathbb E[N]}{\bar D} - \frac{\mathbb E[N D]}{\bar D^2} + \dots = -\frac{\mathbb E[ND]}{\bar D^2} + \dots
\end{split}
\end{equation*}
puisque $\mathbb E[N] = 0$ et $\mathbb E[N\bar D] = 0$. Il reste à calculer
$\mathbb E[ND] = \sum_{t=2}^T\sum_{s=2}^T\mathbb E\bigl[y_{t-1}\varepsilon_t\,y_{s-1}^2\bigr]$.
Si $s\leq t$, $y_{s-1}$ ne dépend que des innovations antérieures à
$\varepsilon_t$ et le terme est nul. Si $s > t$, on isole dans $y_{s-1}$ la
contribution de $y_{t-1}$ et de \(\varepsilon_t\)\nbsp{}:
\begin{equation*}
y_{s-1} = \rho^{s-t}y_{t-1} + \rho^{s-t-1}\varepsilon_t + u_{s,t}
\qquad\text{avec}\qquad
u_{s,t} = \sum_{j=t+1}^{s-1}\rho^{s-1-j}\varepsilon_j
\end{equation*}
où $u_{s,t}$ est centré et indépendant de $(y_{t-1},\varepsilon_t)$. En
développant le carré $y_{s-1}^2$ et en multipliant par
$y_{t-1}\varepsilon_t$, tous les termes ont une espérance nulle (ils
contiennent $\varepsilon_t$ ou $u_{s,t}$ à une puissance impaire, ou
$y_{t-1}$ à une puissance impaire avec $\varepsilon_t$ à une puissance
paire), sauf le double produit\nbsp{}:
\begin{equation*}
\mathbb E\bigl[y_{t-1}\varepsilon_t\,y_{s-1}^2\bigr] = 2\rho^{2(s-t)-1}\,\mathbb E\bigl[y_{t-1}^2\bigr]\,\mathbb E\bigl[\varepsilon_t^2\bigr] = 2\rho^{2(s-t)-1}\sigma^2\,\mathbb E\bigl[y_{t-1}^2\bigr]
\end{equation*}
En sommant sur $s = t+1,\dots,T$, on obtient la série géométrique
$\rho(1-\rho^{2(T-t)})/(1-\rho^2)$. À l'ordre dominant, on peut remplacer
$\mathbb E[y_{t-1}^2]$ par la variance stationnaire $\gamma_0$ et la somme
géométrique par sa limite $\rho/(1-\rho^2)$, les écarts ne concernant qu'un
nombre borné de termes aux bords de l'échantillon, ce qui ne modifie que les
termes d'ordre $1/T^2$ du biais\nbsp{}:
\begin{equation*}
\mathbb E[ND] \simeq (T-1)\,\frac{2\rho\sigma^2\gamma_0}{1-\rho^2} = 2\rho\,(T-1)\,\gamma_0^2
\qquad\text{et}\qquad
\bar D \simeq (T-1)\,\gamma_0
\end{equation*}
où l'on a utilisé $\sigma^2/(1-\rho^2) = \gamma_0$. Finalement\nbsp{}:
\begin{equation*}
\mathbb E[\hat\rho] - \rho \simeq -\frac{2\rho\,(T-1)\,\gamma_0^2}{(T-1)^2\gamma_0^2} = -\frac{2\rho}{T-1} = -\frac{2\rho}{T} + O\!\left(\frac{1}{T^2}\right)
\end{equation*}
#+END_proof

#+BEGIN_remarque
Le biais est négatif pour $\rho > 0$, proportionnel à $\rho$ et inversement
proportionnel à la taille de l'échantillon. Il est nul au premier ordre
lorsque $\rho = 0$. Il ne dépend pas de $\sigma^2$, conformément à
l'[[ex-sigma][exercice]] de la première section. Pour un échantillon de
$50$ observations et $\rho = 0{,}9$, l'approximation donne un biais de
$-0{,}036$, ce qui n'est pas négligeable\nbsp{}: la demi-vie d'un choc,
$\log(1/2)/\log\rho$, passe de $6{,}6$ périodes avec $\rho = 0{,}9$ à
$4{,}7$ périodes avec $\rho = 0{,}864$.
#+END_remarque

#+BEGIN_remarque
Lorsque la régression comporte une constante, c'est-à-dire lorsque la
moyenne du processus est estimée en même temps que $\rho$, le biais est
plus important. Kendall (/Biometrika/, 1954) a montré qu'au premier ordre
il vaut $-(1+3\rho)/T$. Même sous l'hypothèse $\rho = 0$, l'autocorrélation
empirique d'ordre un est donc biaisée vers le bas, de $-1/T$. Ce cas est
traité dans un [[ex-constante][exercice]] à la fin de la note.
#+END_remarque

* Simulations de Monte Carlo
:PROPERTIES:
:CUSTOM_ID: monte-carlo
:END:

Pour mesurer le biais on procède par simulation. On tire $B$ échantillons
de taille $T$ du processus, on calcule l'estimateur sur chacun d'eux, et la
moyenne des $B$ estimations, diminuée de la vraie valeur, approche le biais
avec une erreur d'ordre $1/\sqrt B$. La fonction suivante simule $B$
trajectoires à la fois\nbsp{}: chaque ligne du tableau ~y~ est un échantillon,
et la boucle ne porte que sur le temps.

#+begin_src python :session ar1 :exports code :results none
  rng = np.random.default_rng(2014)

  def simule_ar1(T, rho, B, sigma=1.0):
      epsilon = sigma*rng.standard_normal((B, T))
      y = np.zeros((B, T+1))
      for t in range(1, T+1):
          y[:, t] = rho*y[:, t-1] + epsilon[:, t-1]
      return y[:, 1:]

  def estime_ar1(y):
      return (y[:, 1:]*y[:, :-1]).sum(axis=1)/(y[:, :-1]**2).sum(axis=1)
#+end_src

On commence par $\rho = 0{,}9$ et des tailles d'échantillon de $10$ à
$460$, avec $B = 20\,000$ réplications.

#+begin_src python :session ar1 :exports code :results none
  B = 20000
  rho0 = 0.9
  for T in range(10, 501, 50):
      biais = estime_ar1(simule_ar1(T, rho0, B)).mean() - rho0
      print(T, round(biais, 4), round(-2*rho0/T, 4))
#+end_src

| $T$ | Biais simulé | $-2\rho/T$ |
|-----+--------------+------------|
|  10 |     $-0{,}134$ | $-0{,}180$   |
|  60 |     $-0{,}028$ | $-0{,}030$   |
| 110 |     $-0{,}016$ | $-0{,}016$   |
| 160 |     $-0{,}011$ | $-0{,}011$   |
| 210 |     $-0{,}008$ | $-0{,}009$   |
| 260 |     $-0{,}007$ | $-0{,}007$   |
| 310 |     $-0{,}006$ | $-0{,}006$   |
| 360 |     $-0{,}005$ | $-0{,}005$   |
| 410 |     $-0{,}004$ | $-0{,}004$   |
| 460 |     $-0{,}004$ | $-0{,}004$   |

L'approximation au premier ordre est excellente dès que l'échantillon
compte une soixantaine d'observations. Pour les très petits échantillons,
elle surestime le biais\nbsp{}: avec $T = 10$, le biais simulé est de
$-0{,}134$ contre $-0{,}180$ prédit, les termes d'ordre $1/T^2$ n'étant plus
négligeables. La figure [[fig:biais-T][ci-dessous]] reprend cette
comparaison pour deux valeurs de $\rho$.

#+begin_src python :session ar1 :exports none :results none
  T_grille = np.array([10, 15, 20, 30, 40, 50, 75, 100, 150, 200, 300, 500])
  fig, ax = plt.subplots(figsize=(6, 4))
  for rho, c in ((0.5, 'b'), (0.9, 'r')):
      biais = [estime_ar1(simule_ar1(T, rho, B)).mean() - rho for T in T_grille]
      ax.plot(T_grille, biais, c+'o', markersize=4, label=fr'$\rho = {rho}$ (Monte Carlo)')
      ax.plot(T_grille, -2*rho/T_grille, c+'--', linewidth=1, label=fr'$-2\rho/T$')
  ax.axhline(0, color='k', linewidth=0.5)
  ax.set_xscale('log')
  ax.set_xlabel(r'$T$'); ax.set_ylabel(r'$\mathbb{E}[\hat\rho] - \rho$')
  ax.legend()
  fig.tight_layout()
  fig.savefig("ar1-biais-T.svg", transparent=True)
#+end_src

#+CAPTION: *Biais de l'estimateur des MCO en fonction de la taille de l'échantillon (échelle logarithmique), pour deux valeurs du coefficient autorégressif. Les points sont obtenus par simulation, les pointillés sont l'approximation au premier ordre.*
#+LABEL: fig:biais-T
[[file:ar1-biais-T.svg]]

** Effet de la persistance
:PROPERTIES:
:CUSTOM_ID: persistance
:END:

À taille d'échantillon donnée, le biais croît avec $\rho$, linéairement
d'après l'approximation. Pour tracer le biais en fonction de $\rho$ sans
que le bruit de simulation ne brouille la courbe, on utilise les /mêmes/
innovations pour toutes les valeurs de \(\rho\)\nbsp{}: les estimations
obtenues pour deux valeurs voisines de $\rho$ sont alors fortement
corrélées, et leur différence est mesurée bien plus précisément que chacune
d'elles.

#+begin_src python :session ar1 :exports code :results none
  def simule_ar1_avec(epsilon, rho):
      B, T = epsilon.shape
      y = np.zeros((B, T+1))
      for t in range(1, T+1):
          y[:, t] = rho*y[:, t-1] + epsilon[:, t-1]
      return y[:, 1:]

  rho_grille = np.linspace(0, 0.98, 50)
  biais_rho = {}
  for T in (25, 50, 100):
      epsilon = rng.standard_normal((B, T))
      biais_rho[T] = [estime_ar1(simule_ar1_avec(epsilon, rho)).mean() - rho for rho in rho_grille]
#+end_src

#+begin_src python :session ar1 :exports none :results none
  fig, ax = plt.subplots(figsize=(6, 4))
  for T, c in ((25, 'b'), (50, 'g'), (100, 'r')):
      ax.plot(rho_grille, biais_rho[T], c, linewidth=1, label=fr'$T = {T}$')
      ax.plot(rho_grille, -2*rho_grille/T, c+'--', linewidth=1)
  ax.axhline(0, color='k', linewidth=0.5)
  ax.set_xlabel(r'$\rho$'); ax.set_ylabel(r'$\mathbb{E}[\hat\rho] - \rho$')
  ax.legend()
  fig.tight_layout()
  fig.savefig("ar1-biais-rho.svg", transparent=True)
#+end_src

#+CAPTION: *Biais de l'estimateur des MCO en fonction du coefficient autorégressif, pour trois tailles d'échantillon. Traits pleins : simulations ; pointillés : approximation $-2\rho/T$.*
#+LABEL: fig:biais-rho
[[file:ar1-biais-rho.svg]]

La figure [[fig:biais-rho][ci-dessus]] confirme la proportionnalité à
$\rho$ pour les valeurs modérées de la persistance, et montre que
l'approximation se dégrade lorsque $\rho$ approche de un, surtout dans les
petits échantillons\nbsp{}: le biais simulé s'incurve et reste inférieur en
valeur absolue à $2\rho/T$. Au voisinage de la racine unitaire, la
distribution de l'estimateur change de nature et le développement en
puissances de $1/T$ n'est plus adapté.

** Distribution de l'estimateur
:PROPERTIES:
:CUSTOM_ID: distribution
:END:

Le biais n'est qu'un aspect de la distribution de $\hat\rho$ en
échantillon fini. La figure [[fig:densites][ci-dessous]] compare la
densité de l'estimateur, estimée par noyau à partir des $B$ réplications,
avec l'approximation asymptotique $\mathcal N\bigl(\rho, (1-\rho^2)/T\bigr)$.

#+begin_src python :session ar1 :exports code :results none
  fig, axes = plt.subplots(1, 3, figsize=(12, 4))
  x = np.linspace(-0.6, 1.2, 600)
  for ax, rho in zip(axes, (0.0, 0.5, 0.9)):
      for T, c in ((25, 'b'), (100, 'r')):
          rhohat = estime_ar1(simule_ar1(T, rho, B))
          ax.plot(x, gaussian_kde(rhohat)(x), c, linewidth=1, label=fr'$T = {T}$')
          ax.plot(x, norm.pdf(x, rho, np.sqrt((1-rho**2)/T)), c+':', linewidth=1,
                  label=fr'$\mathcal{{N}}(\rho, (1-\rho^2)/{T})$')
      ax.axvline(rho, color='k', linewidth=0.5)
      ax.set_xlim(rho-0.7, min(rho+0.7, 1.15))
      ax.set_xlabel(r'$\hat\rho$'); ax.set_title(fr'$\rho = {rho}$')
      ax.legend(fontsize=8)
  fig.tight_layout()
  fig.savefig("ar1-densites.svg", transparent=True)
#+end_src

#+CAPTION: *Densité de l'estimateur des MCO (traits pleins, estimation par noyau sur 20 000 réplications) et loi asymptotique (pointillés), pour trois valeurs de $\rho$ et deux tailles d'échantillon. Le trait vertical marque la vraie valeur.*
#+LABEL: fig:densites
[[file:ar1-densites.svg]]

Pour $\rho = 0$ la loi asymptotique est une bonne approximation même avec
$25$ observations. Pour $\rho = 0{,}5$ la distribution est décalée vers la
gauche, c'est le biais, et légèrement asymétrique. Pour $\rho = 0{,}9$
l'asymétrie devient marquée\nbsp{}: la distribution a une longue queue à
gauche et est comprimée à droite, l'estimateur dépassant rarement un. La
loi normale, symétrique autour de $\rho$, décrit alors mal l'incertitude
sur l'estimation, et un intervalle de confiance construit sur la loi
asymptotique est trop optimiste du côté des faibles persistances.

* Corriger le biais
:PROPERTIES:
:CUSTOM_ID: correction
:END:

Puisque $\mathbb E[\hat\rho] \simeq \rho(1-2/T)$, un estimateur corrigé
s'obtient immédiatement\nbsp{}:

\begin{equation*}
\tilde\rho = \frac{T}{T-2}\,\hat\rho
\end{equation*}

dont l'espérance est $\rho + O(1/T^2)$. Cette correction analytique a un
coût\nbsp{}: la variance de $\tilde\rho$ est celle de $\hat\rho$ multipliée
par $(T/(T-2))^2$. Sur les petits échantillons, où la correction est utile,
cette inflation n'est pas négligeable, et le gain doit se mesurer en erreur
quadratique moyenne.

Une seconde approche, qui ne suppose pas connue la forme du biais, est le
/bootstrap paramétrique/. À partir d'un échantillon on calcule
$\hat\rho$, on simule $R$ échantillons artificiels du processus
$y_t = \hat\rho\,y_{t-1} + \varepsilon_t$, on ré-estime $\rho$ sur chacun
d'eux, et la moyenne des $R$ estimations $\hat\rho^{\star}_r$, diminuée de
$\hat\rho$, estime le biais. L'estimateur corrigé par bootstrap
est\nbsp{}:

\begin{equation*}
\hat\rho_{\text{boot}} = \hat\rho - \Bigl(\frac1R\sum_{r=1}^R\hat\rho^{\star}_r - \hat\rho\Bigr) = 2\hat\rho - \frac1R\sum_{r=1}^R\hat\rho^{\star}_r
\end{equation*}

Le biais est ainsi évalué en $\hat\rho$ plutôt qu'en $\rho$, ce qui
introduit une erreur d'ordre $1/T^2$ seulement, la fonction de biais étant
régulière. La fonction suivante applique la correction à $B$ échantillons
simultanément, en emboîtant les $R$ réplications bootstrap dans une
dimension supplémentaire.

#+begin_src python :session ar1 :exports code :results none
  def bootstrap(y, R=199):
      B, T = y.shape
      rhohat = estime_ar1(y)
      rhostar = np.empty((B, R))
      for r in range(R):
          z = simule_ar1_avec(rng.standard_normal((B, T)), rhohat)
          rhostar[:, r] = estime_ar1(z)
      return 2*rhohat - rhostar.mean(axis=1)

  for T in (25, 50, 100):
      y = simule_ar1(T, rho0, B)
      estimateurs = {'MCO': estime_ar1(y), 'corrigé': T*estime_ar1(y)/(T-2), 'bootstrap': bootstrap(y)}
      for nom, e in estimateurs.items():
          print(T, nom, round(e.mean() - rho0, 4), round(e.std(), 4), round(np.sqrt(((e - rho0)**2).mean()), 4))
#+end_src

Le tableau suivant reporte, pour $\rho = 0{,}9$, le biais, l'écart-type et
la racine de l'erreur quadratique moyenne (REQM) des trois estimateurs.

| $T$ | Estimateur | Biais    | Écart-type | REQM     |
|-----+------------+----------+------------+----------|
|  25 | MCO        | $-0{,}065$ | $0{,}134$    | $0{,}149$  |
|  25 | Corrigé    | $+0{,}007$ | $0{,}146$    | $0{,}146$  |
|  25 | Bootstrap  | $-0{,}006$ | $0{,}142$    | $0{,}142$  |
|  50 | MCO        | $-0{,}034$ | $0{,}080$    | $0{,}087$  |
|  50 | Corrigé    | $+0{,}002$ | $0{,}083$    | $0{,}083$  |
|  50 | Bootstrap  | $-0{,}001$ | $0{,}083$    | $0{,}083$  |
| 100 | MCO        | $-0{,}018$ | $0{,}051$    | $0{,}054$  |
| 100 | Corrigé    | $+0{,}000$ | $0{,}052$    | $0{,}052$  |
| 100 | Bootstrap  | $-0{,}001$ | $0{,}052$    | $0{,}052$  |

Les deux corrections éliminent l'essentiel du biais, au prix d'une
augmentation modeste de l'écart-type, et réduisent l'erreur quadratique
moyenne. Le bootstrap fait légèrement mieux que la correction analytique
dans le plus petit échantillon, où l'approximation $-2\rho/T$ surestime le
biais, comme on l'a vu\nbsp{}; les deux méthodes sont équivalentes au-delà
de $50$ observations. La correction analytique a l'avantage d'être
immédiate, le bootstrap celui de s'étendre sans calcul à des modèles pour
lesquels on ne dispose pas de formule.

* Au-delà de l'AR(1)
:PROPERTIES:
:CUSTOM_ID: generalisations
:END:

Le mécanisme à l'origine du biais ne doit rien à la dimension du modèle
ni à l'ordre du processus\nbsp{}: dès que les régresseurs contiennent des
valeurs retardées de la variable expliquée, ils sont prédéterminés sans
être strictement exogènes, et l'estimateur des MCO, ou du maximum de
vraisemblance, est biaisé à l'ordre $1/T$. Ce qui change avec le modèle,
c'est la possibilité d'écrire le biais explicitement, et la forme que
prennent les difficultés en petit échantillon.

** Processus AR(p)
:PROPERTIES:
:CUSTOM_ID: ar-p
:END:

Pour un processus autorégressif d'ordre $p$, le développement au premier
ordre se mène de la même façon, avec des matrices à la place des
scalaires, et le biais de chaque coefficient est une combinaison linéaire
des coefficients, divisée par $T$. Les formules explicites sont dues à
Tjøstheim et Paulsen (/Biometrika/, 1983) et, pour la régression avec
constante, à Shaman et Stine (/Journal of the American Statistical
Association/, 1988). Le résultat qualitatif est celui de l'AR(1)\nbsp{}: la
persistance est sous-estimée, la somme des coefficients autorégressifs,
qui gouverne la demi-vie des chocs, étant biaisée vers le bas d'une
quantité comparable à $-(1+3\rho)/T$ où $\rho$ serait la somme des
coefficients. Les deux corrections de la section
[[#correction][précédente]] s'appliquent sans changement, la formule
analytique demandant seulement d'évaluer le biais au point estimé.

** Processus VAR(p)
:PROPERTIES:
:CUSTOM_ID: var-p
:END:

Le cas vectoriel est celui où la correction du biais a le plus
d'importance pratique, car les modèles VAR sont estimés sur des
échantillons courts, avec de nombreux paramètres, et servent à calculer
des fonctions de réponse aux chocs qui dépendent de façon non linéaire des
coefficients. Pour un VAR(1) à $n$ variables, $Y_t = c + A\,Y_{t-1} + \varepsilon_t$,
où $\varepsilon_t$ est un bruit blanc de matrice de variance $\Sigma$,
Nicholls et Pope (/Australian Journal of Statistics/, 1988) ont
généralisé la formule de Kendall\nbsp{}:

\begin{equation*}
\begin{split}
\mathbb E[\hat A] - A &\simeq -\frac{1}{T}\,\Sigma\,M\,\Gamma_0^{-1}\\
M &= (I_n - A^\top)^{-1} + A^\top(I_n - A^{\top 2})^{-1} + \sum_{k=1}^n\lambda_k\,(I_n - \lambda_kA^\top)^{-1}
\end{split}
\end{equation*}

où les $\lambda_k$ sont les valeurs propres de $A$ et $\Gamma_0$ la
matrice de variance stationnaire de $Y_t$, solution de
$\Gamma_0 = A\,\Gamma_0A^\top + \Sigma$. Pour $n = 1$, avec
$\Sigma = \sigma^2$, $\Gamma_0 = \sigma^2/(1-\rho^2)$ et $\lambda_1 = \rho$,
la matrice $M$ vaut $(1-\rho)^{-1} + 2\rho(1-\rho^2)^{-1}$ et l'on retrouve
\(-(1+3\rho)/T\)\nbsp{}; sans constante, le premier terme de $M$, qui
provient de l'estimation de la moyenne, disparaît, et il reste $-2\rho/T$.
Un VAR(p) se ramène à un VAR(1) par sa représentation compagnon, et la
formule s'applique à la matrice compagnon, comme l'a montré Pope
(/Journal of Time Series Analysis/, 1990). Kilian (/Review of Economics
and Statistics/, 1998) a proposé de corriger le biais des coefficients, par
cette formule ou par bootstrap, avant de calculer les fonctions de
réponse et leurs intervalles de confiance\nbsp{}: le biais vers le bas de la
persistance se traduit par des réponses trop rapidement amorties, et les
intervalles construits sans correction sont décalés. Un
[[ex-var][exercice]] propose de vérifier la formule sur un VAR bivarié.

** Processus ARMA(p,q)
:PROPERTIES:
:CUSTOM_ID: arma
:END:

Lorsque le modèle comporte une partie moyenne mobile, l'estimateur du
maximum de vraisemblance, exact ou conditionnel, n'a plus de forme
explicite, et son biais en échantillon fini non plus. Les simulations
d'Ansley et Newbold (/Journal of Econometrics/, 1980) montrent des biais
du même ordre que pour les processus autorégressifs, avec une difficulté
supplémentaire propre aux coefficients de moyenne mobile\nbsp{}: la
vraisemblance atteint son maximum sur la frontière d'inversibilité,
$\hat\theta = \pm1$, avec une probabilité strictement positive, même
lorsque le vrai coefficient est nettement à l'intérieur\nbsp{}: pour un
MA(1) avec $\theta = 0{,}8$ et $T = 50$, Cryer et Ledolter
(/Biometrika/, 1981) l'évaluent à $0{,}13$. Ce phénomène d'/accumulation/ (/pile-up/) de la distribution
de l'estimateur sur la frontière, étudié par Sargan et Bhargava
(/Econometrica/, 1983), fait que la distribution de $\hat\theta$ est
bimodale, avec une masse ponctuelle en un, et qu'aucune correction
proportionnelle ne peut la recentrer. Le bootstrap paramétrique reste
utilisable pour estimer le biais, à condition de réestimer le modèle sur
chaque échantillon artificiel, ce qui a un coût numérique, et de
retenir la solution inversible.

** Modèles état-mesure
:PROPERTIES:
:CUSTOM_ID: etat-mesure
:END:

Les paramètres d'un modèle état-mesure sont estimés par le maximum de
vraisemblance, la vraisemblance étant calculée par le filtre de Kalman à
partir des erreurs de prévision. Rien de nouveau au regard du biais,
puisque les modèles ARMA en sont un cas particulier, mais le phénomène
d'accumulation prend une forme qui mérite d'être connue. Dans le modèle
de niveau local, $y_t = \mu_t + \varepsilon_t$ et $\mu_t = \mu_{t-1} + \eta_t$,
l'estimateur du maximum de vraisemblance de la variance de $\eta_t$ est
exactement nul avec une probabilité strictement positive, même si cette
variance est positive (Shephard et Harvey, /Journal of Time Series
Analysis/, 1990). Le niveau est alors estimé constant, et le modèle se
réduit à un bruit blanc autour d'une moyenne. Ce n'est pas un hasard\nbsp{}:
le niveau local est équivalent à un ARIMA(0,1,1), et une variance nulle
de $\eta_t$ correspond à $\theta = -1$, la frontière d'inversibilité.
Pour ce type de paramètre, Stock et Watson (/Journal of the American
Statistical Association/, 1998) ont proposé un estimateur sans biais en
médiane, obtenu en inversant une fonction de la statistique de test de
constance du niveau, calculée par simulation. Le bootstrap des modèles
état-mesure, qui rééchantillonne les innovations standardisées du filtre
de Kalman, a été développé par Stoffer et Wall (/Journal of the American
Statistical Association/, 1991).

#+BEGIN_remarque
Corriger l'espérance n'est pas la seule façon de recentrer un estimateur.
Andrews (/Econometrica/, 1993) a construit pour l'AR(1) un estimateur
/sans biais en médiane/\nbsp{}: on calcule, par simulation, la médiane de
$\hat\rho$ pour chaque valeur de $\rho$, et l'on inverse cette fonction au
point $\hat\rho$ observé. Cette construction reste valable au voisinage
de la racine unitaire, où le développement en puissances de $1/T$ n'est
plus adapté et où la distribution de l'estimateur est très asymétrique,
comme on l'a vu sur la figure des [[fig:densites][densités]]. Elle
s'étend aux modèles avec constante et tendance, et c'est elle que
reprennent Stock et Watson pour les variances des modèles état-mesure.
#+END_remarque

* Exercices
:PROPERTIES:
:CUSTOM_ID: exercices
:END:

#+NAME: ex-constante
#+BEGIN_exercice
On estime maintenant $\rho$ par les MCO dans une régression avec
constante, $y_t = c + \rho\,y_{t-1} + \varepsilon_t$. Écrivez l'estimateur
de $\rho$, simulez son biais pour $\rho = 0$, $0{,}5$ et $0{,}9$ et
$T = 25$, $50$, $100$, et comparez avec l'approximation de Kendall,
$-(1+3\rho)/T$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Dans la régression avec constante, l'estimateur de $\rho$ est le
coefficient de la régression de $y_t$ sur $y_{t-1}$ après avoir centré les
deux variables sur leurs moyennes empiriques (calculées sur $t = 2,\dots,T$
pour $y_t$ et sur $t = 1,\dots,T-1$ pour $y_{t-1}$)\nbsp{}:
\begin{equation*}
\hat\rho_c = \frac{\sum_{t=2}^T(y_t - \bar y_{(1)})(y_{t-1} - \bar y_{(0)})}{\sum_{t=2}^T(y_{t-1} - \bar y_{(0)})^2}
\end{equation*}
#+begin_src python :session ar1 :exports code :results none
  def estime_ar1_constante(y):
      y1 = y[:, 1:] - y[:, 1:].mean(axis=1, keepdims=True)
      y0 = y[:, :-1] - y[:, :-1].mean(axis=1, keepdims=True)
      return (y1*y0).sum(axis=1)/(y0**2).sum(axis=1)

  for rho in (0.0, 0.5, 0.9):
      for T in (25, 50, 100):
          biais = estime_ar1_constante(simule_ar1(T, rho, B)).mean() - rho
          print(rho, T, round(biais, 4), round(-(1+3*rho)/T, 4))
#+end_src
| $\rho$ | $T$ | Biais simulé | $-(1+3\rho)/T$ |
|--------+-----+--------------+----------------|
| $0$    |  25 | $-0{,}041$     | $-0{,}040$       |
| $0$    |  50 | $-0{,}019$     | $-0{,}020$       |
| $0$    | 100 | $-0{,}010$     | $-0{,}010$       |
| $0{,}5$  |  25 | $-0{,}103$     | $-0{,}100$       |
| $0{,}5$  |  50 | $-0{,}053$     | $-0{,}050$       |
| $0{,}5$  | 100 | $-0{,}024$     | $-0{,}025$       |
| $0{,}9$  |  25 | $-0{,}182$     | $-0{,}148$       |
| $0{,}9$  |  50 | $-0{,}087$     | $-0{,}074$       |
| $0{,}9$  | 100 | $-0{,}041$     | $-0{,}037$       |
Le biais est nettement plus grand qu'en l'absence de constante, et il ne
s'annule plus lorsque \(\rho = 0\)\nbsp{}: le centrage introduit une corrélation
négative entre les écarts à la moyenne de $y_t$ et de $y_{t-1}$, même pour
un bruit blanc. L'approximation de Kendall est bonne pour les persistances
modérées, mais sous-estime le biais lorsque $\rho = 0{,}9$, en particulier
dans les petits échantillons, à l'inverse de ce que l'on observait sans
constante.
#+END_details

#+NAME: ex-stationnaire
#+BEGIN_exercice
Les simulations démarrent en $y_0 = 0$. Reprenez le calcul du biais de
$\hat\rho$ pour $\rho = 0{,}9$ et $T = 25$, $50$, $100$ en tirant la condition
initiale dans la loi stationnaire, $y_0\sim\mathcal N(0, \gamma_0)$, et
comparez.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Il suffit d'initialiser la première colonne du tableau avec des tirages
d'écart-type $\sqrt{\gamma_0} = \sigma/\sqrt{1-\rho^2}$.
#+begin_src python :session ar1 :exports code :results none
  def simule_ar1_stationnaire(T, rho, B, sigma=1.0):
      epsilon = sigma*rng.standard_normal((B, T))
      y = np.zeros((B, T+1))
      y[:, 0] = sigma/np.sqrt(1-rho**2)*rng.standard_normal(B)
      for t in range(1, T+1):
          y[:, t] = rho*y[:, t-1] + epsilon[:, t-1]
      return y[:, 1:]

  for T in (25, 50, 100):
      b0 = estime_ar1(simule_ar1(T, rho0, B)).mean() - rho0
      b1 = estime_ar1(simule_ar1_stationnaire(T, rho0, B)).mean() - rho0
      print(T, round(b0, 4), round(b1, 4))
#+end_src
| $T$ | $y_0 = 0$ | $y_0$ stationnaire |
|-----+-----------+--------------------|
|  25 | $-0{,}065$  | $-0{,}058$           |
|  50 | $-0{,}034$  | $-0{,}032$           |
| 100 | $-0{,}017$  | $-0{,}017$           |
Partir de la loi stationnaire réduit légèrement le biais, d'un peu moins
d'un centième pour $T = 25$, de deux millièmes pour $T = 50$, et l'écart
n'est plus mesurable pour $T = 100$, l'erreur de simulation étant d'environ
$0{,}001$ avec $20\,000$ réplications. L'écart est divisé par quatre à
chaque doublement de $T$, ce qui est bien le comportement d'un terme en
$1/T^2$. Avec $\rho = 0{,}9$ le processus met pourtant plusieurs dizaines de
périodes à atteindre sa variance stationnaire, mais la condition initiale
ne modifie que les premiers termes des sommes, dont le poids est d'ordre
$1/T$, d'où un effet d'ordre $1/T^2$ sur le biais, conformément à la preuve
de la propriété.
#+END_details

#+NAME: ex-var
#+BEGIN_exercice
Vérifiez par simulation la formule de Nicholls et Pope pour le VAR(1)
bivarié avec constante défini par\nbsp{}:
\begin{equation*}
A = \begin{pmatrix} 0{,}6 & 0{,}2\\ -0{,}1 & 0{,}8\end{pmatrix}
\qquad\text{et}\qquad
\Sigma = \begin{pmatrix} 1 & 0{,}3\\ 0{,}3 & 1\end{pmatrix}
\end{equation*}
pour $T = 50$, $100$ et $200$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
L'estimateur des MCO de $A$ est $\hat A = \bigl(\sum_tY_tY_{t-1}^\top\bigr)\bigl(\sum_tY_{t-1}Y_{t-1}^\top\bigr)^{-1}$,
les variables étant centrées sur leurs moyennes empiriques puisque la
régression comporte une constante. La matrice $\Gamma_0$ s'obtient en
résolvant l'équation de Lyapunov discrète, ce que fait
~solve_discrete_lyapunov~ de ~scipy~.
#+begin_src python :session ar1 :exports code :results none
  from scipy.linalg import solve_discrete_lyapunov

  A = np.array([[0.6, 0.2], [-0.1, 0.8]])
  Sigma = np.array([[1.0, 0.3], [0.3, 1.0]])
  n = 2

  def simule_var1(T, B):
      epsilon = rng.standard_normal((B, T, n)) @ np.linalg.cholesky(Sigma).T
      y = np.zeros((B, T+1, n))
      for t in range(1, T+1):
          y[:, t] = y[:, t-1] @ A.T + epsilon[:, t-1]
      return y[:, 1:]

  def estime_var1(y):
      y1 = y[:, 1:] - y[:, 1:].mean(axis=1, keepdims=True)
      y0 = y[:, :-1] - y[:, :-1].mean(axis=1, keepdims=True)
      S10 = np.einsum('bti,btj->bij', y1, y0)
      S00 = np.einsum('bti,btj->bij', y0, y0)
      return S10 @ np.linalg.inv(S00)

  Gamma0 = solve_discrete_lyapunov(A, Sigma)
  I, At = np.eye(n), A.T
  crochet = np.linalg.inv(I - At) + At @ np.linalg.inv(I - At @ At)
  crochet = crochet + sum(lam*np.linalg.inv(I - lam*At) for lam in np.linalg.eigvals(A))
  biais_pope = -(Sigma @ crochet @ np.linalg.inv(Gamma0)).real

  for T in (50, 100, 200):
      biais_mc = estime_var1(simule_var1(T, B)).mean(axis=0) - A
      print(T, np.round(biais_mc, 3).tolist(), np.round(biais_pope/T, 3).tolist())
#+end_src
| $T$ | Biais simulé                       | Nicholls et Pope                  |
|-----+------------------------------------+-----------------------------------|
|  50 | \(\begin{pmatrix} -0{,}073 & 0{,}022\\ -0{,}019 & -0{,}067\end{pmatrix}\) | \(\begin{pmatrix} -0{,}070 & 0{,}024\\ -0{,}020 & -0{,}059\end{pmatrix}\) |
| 100 | \(\begin{pmatrix} -0{,}036 & 0{,}011\\ -0{,}009 & -0{,}033\end{pmatrix}\) | \(\begin{pmatrix} -0{,}035 & 0{,}012\\ -0{,}010 & -0{,}030\end{pmatrix}\) |
| 200 | \(\begin{pmatrix} -0{,}018 & 0{,}006\\ -0{,}005 & -0{,}015\end{pmatrix}\) | \(\begin{pmatrix} -0{,}017 & 0{,}006\\ -0{,}005 & -0{,}015\end{pmatrix}\) |
L'accord est bon dès $T = 100$, et l'ordre de grandeur est le même que
pour l'AR(1)\nbsp{}: les deux coefficients diagonaux, qui portent la
persistance de chaque variable, sont biaisés vers le bas de $-0{,}07$
environ pour $T = 50$, soit un dixième de leur valeur. Le biais du
coefficient croisé $a_{1,2}$ est positif\nbsp{}: dans le cas vectoriel, le
biais d'un coefficient dépend de toute la matrice $A$, et son signe n'est
pas nécessairement opposé à celui du coefficient.
#+END_details