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 : 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 : 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.
Le modèle et l'estimateur
On considère un processus autorégressif d'ordre un :
\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 : 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 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 :
\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 :
\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.
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.
Corrigé
É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 :
\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 : \(\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é.
Convergence et loi asymptotique
L'estimateur des MCO est convergent et asymptotiquement normal :
\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*}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\) ; 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\)1 :
\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\).
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 : 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\).
Un biais en échantillon fini
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\) : l'innovation \(\varepsilon_t\) entre dans \(y_t\), \(y_{t+1}\), …, 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 :
\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 : une innovation positive augmente \(y_t\), \(y_{t+1}\), …, 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.
Au premier ordre en \(1/T\), le biais de l'estimateur des MCO est2 :
\begin{equation*} \mathbb E[\hat\rho] - \rho = -\frac{2\rho}{T} + O\!\left(\frac{1}{T^2}\right) \end{equation*}Nous donnons la dérivation par un développement au premier ordre, dont nous admettons la validité3. 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 :
\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\) :
\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 :
\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 :
\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 :
\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*}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'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 : 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\).
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 exercice à la fin de la note.
Simulations de Monte Carlo
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 : chaque ligne du tableau y est un échantillon,
et la boucle ne porte que sur le temps.
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)
On commence par \(\rho = 0{,}9\) et des tailles d'échantillon de \(10\) à \(460\), avec \(B = 20\,000\) réplications.
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))
| \(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 : 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 ci-dessous reprend cette comparaison pour deux valeurs de \(\rho\).
Effet de la persistance
À 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\) : 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.
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]
La figure 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 : 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
Le biais n'est qu'un aspect de la distribution de \(\hat\rho\) en échantillon fini. La figure 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)\).
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)
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 : 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
Puisque \(\mathbb E[\hat\rho] \simeq \rho(1-2/T)\), un estimateur corrigé s'obtient immédiatement :
\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 : 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 :
\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.
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))
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 ; 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)
Le mécanisme à l'origine du biais ne doit rien à la dimension du modèle ni à l'ordre du processus : 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)
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) : 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 précédente s'appliquent sans changement, la formule analytique demandant seulement d'évaluer le biais au point estimé.
Processus VAR(p)
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 :
\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\) ; 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 : 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 exercice propose de vérifier la formule sur un VAR bivarié.
Processus ARMA(p,q)
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 : 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 : 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
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 : 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).
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 : 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 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.
Exercices
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\).
Corrigé
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}\)) :
\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*}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))
| \(\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\) : 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.
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.
Corrigé
Il suffit d'initialiser la première colonne du tableau avec des tirages d'écart-type \(\sqrt{\gamma_0} = \sigma/\sqrt{1-\rho^2}\).
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))
| \(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é.
Vérifiez par simulation la formule de Nicholls et Pope pour le VAR(1) bivarié avec constante défini par :
\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\).
Corrigé
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.
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())
| \(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) : 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 : 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.
Notes de bas de page:
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.
Ce résultat est dû à Marriott et Pope (Biometrika, 1954) et à White (Biometrika, 1961) ; 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.
La justification rigoureuse du développement, et le calcul des termes d'ordre supérieur, sont dans l'article de White (1961).