# Long-run variance for an AR(1), with simulation

For a stationary (and absolutely summable), mean-zero (just to make this easier) time series $$y_t$$, with $$\gamma_j=\mathrm{Cov}(y_t,y_{t-j})$$, the long-run variance

$$\mathcal{J}=\sum_{j=-\infty}^{\infty}\gamma_j,$$

and it can be shown that

$$\sqrt{T} \bar{y}_T \rightarrow \mathcal{N}(0,\mathcal{J}).$$

For a stationary AR(1) $$y_t = \phi y_{t-1} + \varepsilon_t$$, with $$\varepsilon_t \sim \mathcal{WN}(0,\sigma^2)$$,

$$\mathrm{Var}(y_t)=\gamma_0=\frac{\sigma^2}{1-\phi^2}$$

and

$$\gamma_j = \phi^j\gamma_0$$

and

$$\mathcal{J}=\frac{\sigma^2}{(1-\phi)^2}$$

Everything I said so far should be true and can be verified or found discussed in:

I'd like to check this via a simple simulation below.

I consider an AR(1) with $$\phi=0.9$$ and $$\sigma^2=1$$ so that $$\mathrm{Var}(y_t)=\gamma_0=\frac{1}{1-0.9^2}=5.263158$$ and $$\mathcal{J}=\frac{1}{(1-0.9)^2}=100$$.

You can see the sandwich::lrv() function seems to get the long-run variance right (although, this function reporting $$\frac{1}{T}\mathcal{J}$$ surprised me)

Q: But, why do I have var(mean_hat_stab[(n-window):(n)]) coming out as 0.001800713!?

This calculation is just the variance of the sample mean for very large n. This should be near 100, right? Did I just make a goof in my calculation? Did I miss something fundamental about LRV? Is this just a number of simulations thing?

set.seed(10)
n = 1e7
window = 1000
phi = 0.9 # AR(1) parameter
sig2 = 1
srv = sig2 / (1 - phi^2) # [1] 5.263158
lrv = sig2 / (1 - phi)^2 # [1] 100
eps <- rnorm(n,0,sig2)
y <- c(eps[1],rep(NA,n-1))
for (i in 2:n) {
y[i] = phi*y[i-1]+eps[i]
}
srv_hat <- var(y) # [1] 5.245835 ~ 5.263158, but not that close actually!!
mean_hat <- cumsum(y)/(1:n)
mean_hat_stab <- sqrt(1:n)*mean_hat # This stabilized mean should --> N(0,lrv)
var(mean_hat_stab[(n-window):(n)]) # [1] 0.001800713, Why not 100!?
qqnorm(mean_hat_stab[(n-window):n]) # Why not normal, on a perfect line!?
lrv_hat_a <- sandwich::lrvar(y, type='Andrews') # [1] 9.932345e-06
# Why not 100!?, Why not close to var(mean_hat_stab) above!?
lrv_hat_nw <- sandwich::lrvar(y, type='Newey-West') # [1] 9.921937e-06
# Why not 100!?, Why not close to var(mean_hat_stab) above!?
lrv_hat_a * n # [1] 99.32345
lrv_hat_nw * n # [1] 99.21937
# But these are close to 100!
# So, I guess that these LRV estimators actually return LRV/n, rather than LRV

• The values in mean_hat computed via cumsum are clearly not independent while the estimate of its variance using var assumes independence. That's why the estimate is completely off. Commented Mar 24 at 16:34
• @JarleTufto To make sure I'm understanding you right, in my simulation procedure where I take var(mean_hat_stab), I'm effectively assuming the sample variance of the sample mean $\frac{1}{w} \sum_{i=T-w}^T \left(\frac{1}{\sqrt{i}}(y_1+\ldots+y_i)\right)^2$ is a good estimator of the LRV, $\mathcal{J}$. But it's not because the values in the series $\{\frac{1}{\sqrt{i}}(y_1+\ldots+y_i)\}_i$ aren't actually independent. Is that right? Commented Mar 24 at 16:59
• Yes, that's right. You need to simulate multiple independent realisations of the stochastic process $y_1, y_2, \dots, y_T$ and the corresponding values of $\sqrt{T}\bar y_T$. Commented Mar 24 at 17:16