#+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: en-GB
#+STARTUP: latexpreview
#+TITLE: Bias of the least squares estimator of an AR(1) process
#+DATE: September 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-ar1-biais/bin/python
#+BEGIN_QUOTE
The ordinary least squares estimator of the coefficient of a first-order
autoregressive process is consistent and asymptotically normal, but it is
biased in finite samples: on average, it underestimates the persistence
of the process. This note explains where the bias comes from, namely the
fact that the regressor is not strictly exogenous, derives the classical
approximation $-2\rho/T$ by a first-order expansion, and then confronts
it with Monte Carlo simulations in Python: bias as a function of the
sample size and of the persistence, distribution of the estimator, and
two ways of correcting the bias, analytical and by bootstrap. The
solutions to the exercises are given in collapsible blocks.
#+END_QUOTE
\\
\\
\\
#+BEGIN_SRC bash :results silent :exports none :async t
python3 -m venv /tmp/blog-ar1-biais
source /tmp/blog-ar1-biais/bin/activate
pip install numpy scipy matplotlib
#+END_SRC
#+begin_src python :session ar1-en :exports none :results none
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde, norm
#+end_src
* The model and the estimator
:PROPERTIES:
:CUSTOM_ID: modele
:END:
Consider a first-order autoregressive process:
\begin{equation*}
y_t = \rho\,y_{t-1} + \varepsilon_t,\qquad t = 1,\dots,T
\end{equation*}
where $|\rho| < 1$ and the innovations $\varepsilon_t$ are independent and
identically distributed according to a normal distribution with mean zero
and variance $\sigma^2$. We assume that the process starts at $y_0 = 0$,
so that $y_1 = \varepsilon_1$. The process is then not stationary but
asymptotically stationary: the variance of $y_t$,
$\sigma^2(1-\rho^{2t})/(1-\rho^2)$, converges to the stationary variance
$\gamma_0 = \sigma^2/(1-\rho^2)$. We shall see in an
[[ex-stationnaire][exercise]] that this choice of initial condition has
practically no effect on the results.
We observe a sample $y_1,\dots,y_T$ and wish to estimate $\rho$. The
ordinary least squares (OLS) estimator, obtained by regressing $y_t$ on
$y_{t-1}$ without an intercept, is:
\begin{equation*}
\hat\rho = \frac{\sum_{t=2}^T y_t\,y_{t-1}}{\sum_{t=2}^T y_{t-1}^2}
\end{equation*}
It is also the maximum likelihood estimator conditional on the first
observation, since the conditional log-likelihood is, up to a constant,
$-\frac{1}{2\sigma^2}\sum_{t=2}^T(y_t - \rho y_{t-1})^2$, which $\hat\rho$
maximises. Substituting $y_t = \rho y_{t-1} + \varepsilon_t$ in the
numerator yields the decomposition that will be used throughout:
\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*}
The estimation error is the ratio of two random variables, a numerator
$N$ and a denominator $D$, both of which depend on the innovations.
#+NAME: ex-sigma
#+BEGIN_exercice
Show that the distribution of $\hat\rho$ does not depend on $\sigma^2$.
The bias of the estimator is therefore a function of $\rho$ and $T$ only.
#+END_exercice
#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
Write $\varepsilon_t = \sigma\eta_t$, where the $\eta_t$ are independent
standard normal variables. By induction from $y_0 = 0$,
$y_t = \sigma\sum_{j=0}^{t-1}\rho^j\eta_{t-j} \equiv \sigma z_t$, where the
process $z_t$ depends only on $\rho$ and on the $\eta_t$. Substituting into
the estimator:
\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*}
The factor $\sigma^2$ cancels: $\hat\rho$ is a function of $\rho$, $T$ and
the $\eta_t$ only, and so is its distribution. In the simulations we may
therefore set $\sigma = 1$ without loss of generality.
#+END_details
* Consistency and asymptotic distribution
:PROPERTIES:
:CUSTOM_ID: asymptotique
:END:
#+BEGIN_property
The OLS estimator is consistent and asymptotically normal:
\begin{equation*}
\hat\rho \xrightarrow[T\to\infty]{p} \rho
\qquad\text{and}\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
Divide the numerator and the denominator of $\hat\rho - \rho$ by $T$.
The process $y_t$ being asymptotically stationary and ergodic, the law of
large numbers gives $\frac1T\sum_{t=2}^Ty_{t-1}^2 \to \gamma_0 = \sigma^2/(1-\rho^2)$
in probability. The term $y_{t-1}\varepsilon_t$ has mean zero, since
$\varepsilon_t$ is independent of $y_{t-1}$, which depends only on past
innovations, and $\mathbb E[\varepsilon_t] = 0$; moreover the terms
$y_{t-1}\varepsilon_t$ are uncorrelated, since for $s < t$ we have
$\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$.
The law of large numbers therefore gives
$\frac1T\sum_{t=2}^Ty_{t-1}\varepsilon_t \to 0$, whence the convergence of
$\hat\rho$ to $\rho$. For the asymptotic distribution, we take for granted
that the central limit theorem applies to the sequence
$y_{t-1}\varepsilon_t$, whose terms have mean zero, are uncorrelated and
have variance $\mathbb E[y_{t-1}^2]\,\mathbb E[\varepsilon_t^2] \to \gamma_0\sigma^2$[fn:: The
sequence $y_{t-1}\varepsilon_t$ is a martingale difference, and the central
limit theorem for martingale differences applies.]:
\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*}
By Slutsky's theorem, $\sqrt T(\hat\rho - \rho)$ converges in distribution
to a centred normal with variance $\gamma_0\sigma^2/\gamma_0^2 = \sigma^2/\gamma_0 = 1-\rho^2$.
#+END_proof
#+BEGIN_remarque
The asymptotic variance $1-\rho^2$ does not depend on $\sigma^2$, in line
with the previous exercise, and it decreases as persistence increases: the
closer $\rho$ is to one, the more dispersed the regressor $y_{t-1}$ is
relative to the innovation, and the more precise the estimator. In large
samples, the distribution of $\hat\rho$ is therefore approximated by
$\mathcal N\bigl(\rho, (1-\rho^2)/T\bigr)$, a distribution centred on the
true value. This approximation says nothing about the finite-sample bias,
which is of order $1/T$ and disappears in the normalisation by $\sqrt T$.
#+END_remarque
* A finite-sample bias
:PROPERTIES:
:CUSTOM_ID: biais
:END:
In the linear regression model with deterministic regressors, the OLS
estimator is unbiased because the expectation of the error is zero
conditional on /all/ the regressors. Here the regressor is random, and
while $\mathbb E[\varepsilon_t\mid y_{t-1}] = 0$ does hold,
$\mathbb E[\varepsilon_t\mid y_1,\dots,y_{T-1}] = 0$ does not: the innovation
$\varepsilon_t$ enters $y_t$, $y_{t+1}$, \dots, which are the regressors of
the following dates. The regressor is /predetermined/, not /strictly
exogenous/. To see how this dependence produces a bias, write the
estimation error as a weighted sum of the innovations:
\begin{equation*}
\hat\rho - \rho = \sum_{t=2}^T w_t\,\varepsilon_t
\qquad\text{with}\qquad
w_t = \frac{y_{t-1}}{\sum_{s=2}^Ty_{s-1}^2}
\end{equation*}
If the weights $w_t$ were independent of the innovations, each term would
have zero expectation. But the weight $w_t$ depends on $\varepsilon_t$
through its denominator: a positive innovation raises $y_t$, $y_{t+1}$,
\dots, hence the denominator, and lowers the weight given to that
innovation. A positive innovation thus receives, on average, a smaller
weight than a negative innovation of the same size, and the weighted sum
is negative in expectation when $\rho > 0$. The following property
quantifies this effect.
#+BEGIN_property
To first order in $1/T$, the bias of the OLS estimator is[fn:: This
result is due to Marriott and Pope (/Biometrika/, 1954) and White
(/Biometrika/, 1961); the bias of the OLS estimator in dynamic models was
first pointed out by Hurwicz (1950), in the Cowles Commission volume
edited by Koopmans, and it is sometimes called the Hurwicz bias.]:
\begin{equation*}
\mathbb E[\hat\rho] - \rho = -\frac{2\rho}{T} + O\!\left(\frac{1}{T^2}\right)
\end{equation*}
#+END_property
#+BEGIN_proof
We give the derivation by a first-order expansion, whose validity we take
for granted[fn:: The rigorous justification of the expansion, and the
computation of the higher-order terms, are in White's (1961) article.].
Write $\bar D = \mathbb E[D]$ and $D = \bar D(1 + \delta)$ with
$\delta = (D-\bar D)/\bar D$, which is of order $1/\sqrt T$. The expansion
$(1+\delta)^{-1} = 1 - \delta + O(\delta^2)$ gives:
\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*}
since $\mathbb E[N] = 0$ and $\mathbb E[N\bar D] = 0$. It remains to compute
$\mathbb E[ND] = \sum_{t=2}^T\sum_{s=2}^T\mathbb E\bigl[y_{t-1}\varepsilon_t\,y_{s-1}^2\bigr]$.
If $s\leq t$, $y_{s-1}$ depends only on innovations prior to
$\varepsilon_t$ and the term vanishes. If $s > t$, we isolate in $y_{s-1}$
the contribution of $y_{t-1}$ and of $\varepsilon_t$:
\begin{equation*}
y_{s-1} = \rho^{s-t}y_{t-1} + \rho^{s-t-1}\varepsilon_t + u_{s,t}
\qquad\text{with}\qquad
u_{s,t} = \sum_{j=t+1}^{s-1}\rho^{s-1-j}\varepsilon_j
\end{equation*}
where $u_{s,t}$ has mean zero and is independent of
$(y_{t-1},\varepsilon_t)$. Expanding the square $y_{s-1}^2$ and multiplying
by $y_{t-1}\varepsilon_t$, all terms have zero expectation (they contain
$\varepsilon_t$ or $u_{s,t}$ to an odd power, or $y_{t-1}$ to an odd power
with $\varepsilon_t$ to an even power), except the cross product:
\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*}
Summing over $s = t+1,\dots,T$ gives the geometric series
$\rho(1-\rho^{2(T-t)})/(1-\rho^2)$. To leading order, we may replace
$\mathbb E[y_{t-1}^2]$ by the stationary variance $\gamma_0$ and the
geometric sum by its limit $\rho/(1-\rho^2)$, the discrepancies affecting
only a bounded number of terms at the edges of the sample, which only
changes the terms of order $1/T^2$ of the bias:
\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{and}\qquad
\bar D \simeq (T-1)\,\gamma_0
\end{equation*}
where we used $\sigma^2/(1-\rho^2) = \gamma_0$. Finally:
\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
The bias is negative for $\rho > 0$, proportional to $\rho$ and inversely
proportional to the sample size. It vanishes to first order when
$\rho = 0$. It does not depend on $\sigma^2$, in line with the
[[ex-sigma][exercise]] of the first section. For a sample of $50$
observations and $\rho = 0.9$, the approximation gives a bias of $-0.036$,
which is not negligible: the half-life of a shock, $\log(1/2)/\log\rho$,
falls from $6.6$ periods with $\rho = 0.9$ to $4.7$ periods with
$\rho = 0.864$.
#+END_remarque
#+BEGIN_remarque
When the regression includes an intercept, that is, when the mean of the
process is estimated along with $\rho$, the bias is larger. Kendall
(/Biometrika/, 1954) showed that to first order it equals $-(1+3\rho)/T$.
Even under the hypothesis $\rho = 0$, the first-order sample
autocorrelation is therefore biased downwards, by $-1/T$. This case is
treated in an [[ex-constante][exercise]] at the end of the note.
#+END_remarque
* Monte Carlo simulations
:PROPERTIES:
:CUSTOM_ID: monte-carlo
:END:
To measure the bias we proceed by simulation. We draw $B$ samples of
size $T$ from the process, compute the estimator on each of them, and the
mean of the $B$ estimates, minus the true value, approximates the bias
with an error of order $1/\sqrt B$. The following function simulates $B$
trajectories at once: each row of the array ~y~ is a sample, and the loop
only runs over time.
#+begin_src python :session ar1-en :exports code :results none
rng = np.random.default_rng(2014)
def simulate_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 estimate_ar1(y):
return (y[:, 1:]*y[:, :-1]).sum(axis=1)/(y[:, :-1]**2).sum(axis=1)
#+end_src
We start with $\rho = 0.9$ and sample sizes from $10$ to $460$, with
$B = 20\,000$ replications.
#+begin_src python :session ar1-en :exports code :results none
B = 20000
rho0 = 0.9
for T in range(10, 501, 50):
bias = estimate_ar1(simulate_ar1(T, rho0, B)).mean() - rho0
print(T, round(bias, 4), round(-2*rho0/T, 4))
#+end_src
| $T$ | Simulated bias | $-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$ |
The first-order approximation is excellent as soon as the sample has
about sixty observations. For very small samples it overstates the bias:
with $T = 10$, the simulated bias is $-0.134$ against a predicted
$-0.180$, the terms of order $1/T^2$ no longer being negligible. The
figure [[fig:biais-T][below]] repeats this comparison for two values of
$\rho$.
#+begin_src python :session ar1-en :exports none :results none
T_grid = 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')):
bias = [estimate_ar1(simulate_ar1(T, rho, B)).mean() - rho for T in T_grid]
ax.plot(T_grid, bias, c+'o', markersize=4, label=fr'$\rho = {rho}$ (Monte Carlo)')
ax.plot(T_grid, -2*rho/T_grid, 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: *Bias of the OLS estimator as a function of the sample size (logarithmic scale), for two values of the autoregressive coefficient. Dots are obtained by simulation, dashed lines are the first-order approximation.*
#+LABEL: fig:biais-T
[[file:ar1-biais-T.svg]]
** Effect of persistence
:PROPERTIES:
:CUSTOM_ID: persistance
:END:
For a given sample size, the bias grows with $\rho$, linearly according to
the approximation. To plot the bias as a function of $\rho$ without the
simulation noise blurring the curve, we use the /same/ innovations for
all values of $\rho$: the estimates obtained for two neighbouring values
of $\rho$ are then strongly correlated, and their difference is measured
much more precisely than each of them.
#+begin_src python :session ar1-en :exports code :results none
def simulate_ar1_with(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_grid = np.linspace(0, 0.98, 50)
bias_rho = {}
for T in (25, 50, 100):
epsilon = rng.standard_normal((B, T))
bias_rho[T] = [estimate_ar1(simulate_ar1_with(epsilon, rho)).mean() - rho for rho in rho_grid]
#+end_src
#+begin_src python :session ar1-en :exports none :results none
fig, ax = plt.subplots(figsize=(6, 4))
for T, c in ((25, 'b'), (50, 'g'), (100, 'r')):
ax.plot(rho_grid, bias_rho[T], c, linewidth=1, label=fr'$T = {T}$')
ax.plot(rho_grid, -2*rho_grid/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: *Bias of the OLS estimator as a function of the autoregressive coefficient, for three sample sizes. Solid lines: simulations; dashed lines: the approximation $-2\rho/T$.*
#+LABEL: fig:biais-rho
[[file:ar1-biais-rho.svg]]
The figure [[fig:biais-rho][above]] confirms the proportionality to
$\rho$ for moderate values of the persistence, and shows that the
approximation deteriorates as $\rho$ approaches one, especially in small
samples: the simulated bias bends and remains smaller in absolute value
than $2\rho/T$. In the vicinity of the unit root, the distribution of the
estimator changes nature and the expansion in powers of $1/T$ is no
longer adequate.
** Distribution of the estimator
:PROPERTIES:
:CUSTOM_ID: distribution
:END:
The bias is only one aspect of the finite-sample distribution of
$\hat\rho$. The figure [[fig:densites][below]] compares the density of
the estimator, estimated by a kernel method from the $B$ replications,
with the asymptotic approximation $\mathcal N\bigl(\rho, (1-\rho^2)/T\bigr)$.
#+begin_src python :session ar1-en :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 = estimate_ar1(simulate_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: *Density of the OLS estimator (solid lines, kernel estimate from 20,000 replications) and asymptotic distribution (dotted lines), for three values of $\rho$ and two sample sizes. The vertical line marks the true value.*
#+LABEL: fig:densites
[[file:ar1-densites.svg]]
For $\rho = 0$ the asymptotic distribution is a good approximation even
with $25$ observations. For $\rho = 0.5$ the distribution is shifted to
the left, which is the bias, and slightly asymmetric. For $\rho = 0.9$ the
asymmetry becomes pronounced: the distribution has a long left tail and
is compressed on the right, the estimator rarely exceeding one. The normal
distribution, symmetric around $\rho$, then poorly describes the
uncertainty about the estimate, and a confidence interval built on the
asymptotic distribution is too optimistic on the low-persistence side.
* Correcting the bias
:PROPERTIES:
:CUSTOM_ID: correction
:END:
Since $\mathbb E[\hat\rho] \simeq \rho(1-2/T)$, a corrected estimator
follows immediately:
\begin{equation*}
\tilde\rho = \frac{T}{T-2}\,\hat\rho
\end{equation*}
whose expectation is $\rho + O(1/T^2)$. This analytical correction has a
cost: the variance of $\tilde\rho$ is that of $\hat\rho$ multiplied by
$(T/(T-2))^2$. In small samples, where the correction is useful, this
inflation is not negligible, and the gain must be measured in mean squared
error.
A second approach, which does not assume the form of the bias to be
known, is the /parametric bootstrap/. From a sample we compute
$\hat\rho$, we simulate $R$ artificial samples of the process
$y_t = \hat\rho\,y_{t-1} + \varepsilon_t$, we re-estimate $\rho$ on each
of them, and the mean of the $R$ estimates $\hat\rho^{\star}_r$, minus
$\hat\rho$, estimates the bias. The bootstrap bias-corrected estimator is:
\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*}
The bias is thus evaluated at $\hat\rho$ rather than at $\rho$, which
introduces an error of order $1/T^2$ only, the bias function being
smooth. The following function applies the correction to $B$ samples
simultaneously, nesting the $R$ bootstrap replications in an additional
dimension.
#+begin_src python :session ar1-en :exports code :results none
def bootstrap(y, R=199):
B, T = y.shape
rhohat = estimate_ar1(y)
rhostar = np.empty((B, R))
for r in range(R):
z = simulate_ar1_with(rng.standard_normal((B, T)), rhohat)
rhostar[:, r] = estimate_ar1(z)
return 2*rhohat - rhostar.mean(axis=1)
for T in (25, 50, 100):
y = simulate_ar1(T, rho0, B)
estimators = {'OLS': estimate_ar1(y), 'corrected': T*estimate_ar1(y)/(T-2), 'bootstrap': bootstrap(y)}
for name, e in estimators.items():
print(T, name, round(e.mean() - rho0, 4), round(e.std(), 4), round(np.sqrt(((e - rho0)**2).mean()), 4))
#+end_src
The following table reports, for $\rho = 0.9$, the bias, the standard
deviation and the root mean squared error (RMSE) of the three estimators.
| $T$ | Estimator | Bias | Std. dev. | RMSE |
|-----+-----------+----------+-----------+---------|
| 25 | OLS | $-0.065$ | $0.134$ | $0.149$ |
| 25 | Corrected | $+0.007$ | $0.146$ | $0.146$ |
| 25 | Bootstrap | $-0.006$ | $0.142$ | $0.142$ |
| 50 | OLS | $-0.034$ | $0.080$ | $0.087$ |
| 50 | Corrected | $+0.002$ | $0.083$ | $0.083$ |
| 50 | Bootstrap | $-0.001$ | $0.083$ | $0.083$ |
| 100 | OLS | $-0.018$ | $0.051$ | $0.054$ |
| 100 | Corrected | $+0.000$ | $0.052$ | $0.052$ |
| 100 | Bootstrap | $-0.001$ | $0.052$ | $0.052$ |
Both corrections remove most of the bias, at the price of a modest
increase in the standard deviation, and reduce the mean squared error.
The bootstrap does slightly better than the analytical correction in the
smallest sample, where the approximation $-2\rho/T$ overstates the bias,
as we saw; the two methods are equivalent beyond $50$ observations. The
analytical correction has the advantage of being immediate, the bootstrap
that of extending without further calculation to models for which no
formula is available.
* Beyond the AR(1)
:PROPERTIES:
:CUSTOM_ID: generalisations
:END:
The mechanism behind the bias owes nothing to the dimension of the model
or to the order of the process: as soon as the regressors contain lagged
values of the dependent variable, they are predetermined without being
strictly exogenous, and the OLS or maximum likelihood estimator is biased
to order $1/T$. What changes with the model is the possibility of writing
the bias explicitly, and the form taken by the small-sample difficulties.
** AR(p) processes
:PROPERTIES:
:CUSTOM_ID: ar-p
:END:
For an autoregressive process of order $p$, the first-order expansion
proceeds in the same way, with matrices in place of scalars, and the bias
of each coefficient is a linear combination of the coefficients, divided
by $T$. The explicit formulas are due to Tjøstheim and Paulsen
(/Biometrika/, 1983) and, for the regression with an intercept, to Shaman
and Stine (/Journal of the American Statistical Association/, 1988). The
qualitative result is that of the AR(1): persistence is underestimated,
the sum of the autoregressive coefficients, which governs the half-life
of shocks, being biased downwards by an amount comparable to
$-(1+3\rho)/T$ where $\rho$ would be the sum of the coefficients. The two
corrections of the [[#correction][previous section]] apply unchanged, the
analytical formula only requiring the bias to be evaluated at the
estimated point.
** VAR(p) processes
:PROPERTIES:
:CUSTOM_ID: var-p
:END:
The vector case is the one where bias correction matters most in
practice, since VAR models are estimated on short samples, with many
parameters, and are used to compute impulse response functions that
depend non-linearly on the coefficients. For a VAR(1) in $n$ variables,
$Y_t = c + A\,Y_{t-1} + \varepsilon_t$, where $\varepsilon_t$ is a white
noise with variance matrix $\Sigma$, Nicholls and Pope (/Australian
Journal of Statistics/, 1988) generalised Kendall's formula:
\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*}
where the $\lambda_k$ are the eigenvalues of $A$ and $\Gamma_0$ the
stationary variance matrix of $Y_t$, the solution of
$\Gamma_0 = A\,\Gamma_0A^\top + \Sigma$. For $n = 1$, with
$\Sigma = \sigma^2$, $\Gamma_0 = \sigma^2/(1-\rho^2)$ and $\lambda_1 = \rho$,
the matrix $M$ equals $(1-\rho)^{-1} + 2\rho(1-\rho^2)^{-1}$ and we recover
$-(1+3\rho)/T$; without an intercept, the first term of $M$,
which comes from the estimation of the mean, disappears, and $-2\rho/T$
remains. A VAR(p) reduces to a VAR(1) through its companion
representation, and the formula applies to the companion matrix, as shown
by Pope (/Journal of Time Series Analysis/, 1990). Kilian (/Review of
Economics and Statistics/, 1998) proposed to correct the bias of the
coefficients, by this formula or by bootstrap, before computing impulse
responses and their confidence intervals: the downward bias in
persistence translates into responses that die out too quickly, and
intervals built without correction are shifted. An [[ex-var][exercise]]
proposes to check the formula on a bivariate VAR.
** ARMA(p,q) processes
:PROPERTIES:
:CUSTOM_ID: arma
:END:
When the model has a moving average part, the maximum likelihood
estimator, exact or conditional, no longer has an explicit form, and
neither does its finite-sample bias. The simulations of Ansley and
Newbold (/Journal of Econometrics/, 1980) show biases of the same order
as for autoregressive processes, with an additional difficulty specific
to moving average coefficients: the likelihood reaches its maximum on
the invertibility boundary, $\hat\theta = \pm1$, with strictly positive
probability, even when the true coefficient is well inside: for an MA(1) with
$\theta = 0.8$ and $T = 50$, Cryer and Ledolter (/Biometrika/, 1981) put
it at $0.13$. This /pile-up/ of the distribution of the estimator on the
boundary, studied by Sargan and Bhargava (/Econometrica/, 1983), makes
the distribution of $\hat\theta$ bimodal, with a point mass at one, and
no proportional correction can recentre it. The parametric bootstrap
remains usable to estimate the bias, provided the model is re-estimated
on each artificial sample, which has a numerical cost, and the invertible
solution is retained.
** State-space models
:PROPERTIES:
:CUSTOM_ID: etat-mesure
:END:
The parameters of a state-space model are estimated by maximum
likelihood, the likelihood being computed by the Kalman filter from the
prediction errors. Nothing new as far as the bias is concerned, since
ARMA models are a special case, but the pile-up phenomenon takes a form
worth knowing. In the local level model, $y_t = \mu_t + \varepsilon_t$ and
$\mu_t = \mu_{t-1} + \eta_t$, the maximum likelihood estimator of the
variance of $\eta_t$ is exactly zero with strictly positive probability,
even if this variance is positive (Shephard and Harvey, /Journal of Time
Series Analysis/, 1990). The level is then estimated as constant, and the
model reduces to a white noise around a mean. This is no coincidence: the
local level model is equivalent to an ARIMA(0,1,1), and a zero variance of
$\eta_t$ corresponds to $\theta = -1$, the invertibility boundary. For
this kind of parameter, Stock and Watson (/Journal of the American
Statistical Association/, 1998) proposed a median-unbiased estimator,
obtained by inverting a function of the test statistic for constancy of
the level, computed by simulation. The bootstrap of state-space models,
which resamples the standardised innovations of the Kalman filter, was
developed by Stoffer and Wall (/Journal of the American Statistical
Association/, 1991).
#+BEGIN_remarque
Correcting the expectation is not the only way to recentre an estimator.
Andrews (/Econometrica/, 1993) built for the AR(1) a /median-unbiased/
estimator: one computes, by simulation, the median of $\hat\rho$ for each
value of $\rho$, and inverts this function at the observed $\hat\rho$.
This construction remains valid in the vicinity of the unit root, where
the expansion in powers of $1/T$ is no longer adequate and the
distribution of the estimator is very asymmetric, as we saw on the figure
of the [[fig:densites][densities]]. It extends to models with an
intercept and a trend, and it is the one Stock and Watson take up for
the variances of state-space models.
#+END_remarque
* Exercises
:PROPERTIES:
:CUSTOM_ID: exercices
:END:
#+NAME: ex-constante
#+BEGIN_exercice
We now estimate $\rho$ by OLS in a regression with an intercept,
$y_t = c + \rho\,y_{t-1} + \varepsilon_t$. Write down the estimator of
$\rho$, simulate its bias for $\rho = 0$, $0.5$ and $0.9$ and $T = 25$,
$50$, $100$, and compare with Kendall's approximation, $-(1+3\rho)/T$.
#+END_exercice
#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
In the regression with an intercept, the estimator of $\rho$ is the
coefficient of the regression of $y_t$ on $y_{t-1}$ after centring both
variables on their sample means (computed over $t = 2,\dots,T$ for $y_t$
and over $t = 1,\dots,T-1$ for $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*}
#+begin_src python :session ar1-en :exports code :results none
def estimate_ar1_intercept(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):
bias = estimate_ar1_intercept(simulate_ar1(T, rho, B)).mean() - rho
print(rho, T, round(bias, 4), round(-(1+3*rho)/T, 4))
#+end_src
| $\rho$ | $T$ | Simulated bias | $-(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$ |
The bias is markedly larger than without an intercept, and it no longer
vanishes when $\rho = 0$: centring introduces a negative correlation
between the deviations from the mean of $y_t$ and of $y_{t-1}$, even for a
white noise. Kendall's approximation is good for moderate persistence, but
understates the bias when $\rho = 0.9$, especially in small samples,
contrary to what we observed without an intercept.
#+END_details
#+NAME: ex-stationnaire
#+BEGIN_exercice
The simulations start at $y_0 = 0$. Redo the computation of the bias of
$\hat\rho$ for $\rho = 0.9$ and $T = 25$, $50$, $100$, drawing the initial
condition from the stationary distribution,
$y_0\sim\mathcal N(0, \gamma_0)$, and compare.
#+END_exercice
#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
It suffices to initialise the first column of the array with draws of
standard deviation $\sqrt{\gamma_0} = \sigma/\sqrt{1-\rho^2}$.
#+begin_src python :session ar1-en :exports code :results none
def simulate_ar1_stationary(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 = estimate_ar1(simulate_ar1(T, rho0, B)).mean() - rho0
b1 = estimate_ar1(simulate_ar1_stationary(T, rho0, B)).mean() - rho0
print(T, round(b0, 4), round(b1, 4))
#+end_src
| $T$ | $y_0 = 0$ | Stationary $y_0$ |
|-----+-----------+------------------|
| 25 | $-0.065$ | $-0.058$ |
| 50 | $-0.034$ | $-0.032$ |
| 100 | $-0.017$ | $-0.017$ |
Starting from the stationary distribution slightly reduces the bias, by a
little less than one hundredth for $T = 25$, by two thousandths for
$T = 50$, and the difference is no longer measurable for $T = 100$, the
simulation error being about $0.001$ with $20\,000$ replications. The
difference is divided by four each time $T$ doubles, which is indeed the
behaviour of a $1/T^2$ term. With $\rho = 0.9$ the process nevertheless
takes several dozen periods to reach its stationary variance, but the
initial condition only affects the first terms of the sums, whose weight is
of order $1/T$, hence an effect of order $1/T^2$ on the bias, in line with
the proof of the property.
#+END_details
#+NAME: ex-var
#+BEGIN_exercice
Check by simulation the Nicholls and Pope formula for the bivariate
VAR(1) with intercept defined by:
\begin{equation*}
A = \begin{pmatrix} 0.6 & 0.2\\ -0.1 & 0.8\end{pmatrix}
\qquad\text{and}\qquad
\Sigma = \begin{pmatrix} 1 & 0.3\\ 0.3 & 1\end{pmatrix}
\end{equation*}
for $T = 50$, $100$ and $200$.
#+END_exercice
#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
The OLS estimator of $A$ is $\hat A = \bigl(\sum_tY_tY_{t-1}^\top\bigr)\bigl(\sum_tY_{t-1}Y_{t-1}^\top\bigr)^{-1}$,
the variables being centred on their sample means since the regression
includes an intercept. The matrix $\Gamma_0$ is obtained by solving the
discrete Lyapunov equation, which ~solve_discrete_lyapunov~ from ~scipy~
does.
#+begin_src python :session ar1-en :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 simulate_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 estimate_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
bracket = np.linalg.inv(I - At) + At @ np.linalg.inv(I - At @ At)
bracket = bracket + sum(lam*np.linalg.inv(I - lam*At) for lam in np.linalg.eigvals(A))
bias_pope = -(Sigma @ bracket @ np.linalg.inv(Gamma0)).real
for T in (50, 100, 200):
bias_mc = estimate_var1(simulate_var1(T, B)).mean(axis=0) - A
print(T, np.round(bias_mc, 3).tolist(), np.round(bias_pope/T, 3).tolist())
#+end_src
| $T$ | Simulated bias | Nicholls and 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}\) |
The agreement is good from $T = 100$ onwards, and the order of magnitude
is the same as for the AR(1): the two diagonal coefficients, which carry
the persistence of each variable, are biased downwards by about $-0.07$
for $T = 50$, a tenth of their value. The bias of the cross coefficient
$a_{1,2}$ is positive: in the vector case, the bias of a coefficient
depends on the whole matrix $A$, and its sign is not necessarily opposite
to that of the coefficient.
#+END_details