How to compute the standard error of the mean of an AR(1) process?

Question

I try to compute the standard error of the mean for a demeaned AR(1) process $x_{t+1} = \rho x_t + \varepsilon_{t+1} =\sum\limits_{i=0}^{\infty} \rho^i \varepsilon_{t+1-i}$

Here is what I did:

$$ \begin{align*} Var(\overline{x}) &= Var\left(\frac{1}{N} \sum\limits_{t=0}^{N-1} x_t\right) \\ &= Var\left(\frac{1}{N} \sum\limits_{t=0}^{N-1} \sum\limits_{i=0}^{\infty} \rho^i \varepsilon_{t-i}\right) \\ &= \frac{1}{N^2} Var\begin{pmatrix} \rho^0 \varepsilon_0 + & \rho^1 \varepsilon_{-1} + & \rho^2 \varepsilon_{-2} + & \cdots & \rho^{\infty} \varepsilon_{-\infty} + \\ \rho^0 \varepsilon_1 + & \rho^1 \varepsilon_{0} + & \rho^2 \varepsilon_{-1} + & \cdots & \rho^{\infty} \varepsilon_{1-\infty} + \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^0\varepsilon_{N-1} + & \rho^1 \varepsilon_{N-2} + & \rho^2 \varepsilon_{N-3} + & \cdots & \rho^{\infty} \varepsilon_{N-1-\infty} + \\ \end{pmatrix} \\ &= \frac{1}{N^2} Var\begin{pmatrix} \rho^0 \varepsilon_{N-1} + \\ (\rho^0 + \rho^1) \varepsilon_{N-2} + \\ (\rho^0 + \rho^1 + \rho^2) \varepsilon_{N-3} + \\ \cdots \\ (\rho^0 + \rho^1 + \rho^2 + \dots + \rho^{N-2}) \varepsilon_{1} + \\ (\rho^0 + \rho^1 + \rho^2 + \dots + \rho^{N-1}) \varepsilon_{0} + \\ (\rho^1 + \rho^2 + \rho^3 + \dots + \rho^{N}) \varepsilon_{-1} + \\ (\rho^2 + \rho^3 + \rho^4 + \dots + \rho^{N+1}) \varepsilon_{-2} + \\ \cdots\\ \end{pmatrix} \\ &= \frac{\sigma_{\varepsilon}^2}{N^2} \begin{pmatrix} \rho^0 + \\ (\rho^0 + \rho^1) + \\ (\rho^0 + \rho^1 + \rho^2) + \\ \cdots \\ (\rho^0 + \rho^1 + \rho^2 + \dots + \rho^{N-2}) + \\ (\rho^0 + \rho^1 + \rho^2 + \dots + \rho^{N-1}) + \\ (\rho^1 + \rho^2 + \rho^3 + \dots + \rho^{N}) + \\ (\rho^2 + \rho^3 + \rho^4 + \dots + \rho^{N+1}) + \\ \cdots\\ \end{pmatrix} \\ &= \frac{N \sigma_{\varepsilon}^2}{N^2} (\rho^0 + \rho^1 + \dots + \rho^{\infty}) \\ &= \frac{\sigma_{\varepsilon}^2}{N} \frac{1}{1 - \rho} \\ \end{align*} $$

Probably, not every step is done in the most obvious way, so let me add some thoughts. In the third row, I just write out to two sum-signs. Here, the matrix has N rows. In the fourth row, I realign the matrix so that there is one row for every epsilon, so the number of rows is infinite here. Note that the last three parts in the matrix have the same number of elements, just differencing by a factor $\rho$ in each row. In the fifth row, I apply the rule that the variance of the sum of independent shocks is the sum of the variances of those shocks and notice that each $\rho^j$ element is summed up $N$ times.

The end result looks neat, but is probably wrong. Why do I think so? Because I run a MCS in R and things don't add up:

nrMCS <- 10000
N <- 100
pers <- 0.9
means <- numeric(nrMCS)
for (i in 1:nrMCS) {
   means[i] <- mean(arima.sim(list(order=c(1,0,0), ar=pers), n = N))
}
#quantile(means, probs=c(0.025, 0.05, 0.5, 0.95, 0.975))
#That is the empirical standard error
sd(means)
0.9459876
#This should be the standard error according to my formula
1/(N*(1-pers))
0.1

Any hints on what I am doing wrong would be great! Or maybe a hint where I can find the correct derivation (I couldn't find anything). Is the problem maybe that I assume independence between the same errors?

$$Var(X + X) = Var(2X) = 4Var(X) \neq 2Var(X)$$

I thought about that, but don't see where I make that erroneous assumption in my derivation.

UPDATE

I forgot to square the rhos, as Nuzhi correctly pointed out. Hence it should look like:

$$ Var(\overline{x}) = \frac{\sigma_{\varepsilon}^2}{N^2} \begin{pmatrix} \rho^{2\times0} + \\ (\rho^0 + \rho^1)^2 + \\ (\rho^0 + \rho^1 + \rho^2)^2 + \\ \cdots \\ (\rho^0 + \rho^1 + \rho^2 + \dots + \rho^{N-2})^2 + \\ (\rho^0 + \rho^1 + \rho^2 + \dots + \rho^{N-1})^2 + \\ (\rho^1 + \rho^2 + \rho^3 + \dots + \rho^{N})^2 + \\ (\rho^2 + \rho^3 + \rho^4 + \dots + \rho^{N+1})^2 + \\ \cdots\\ \end{pmatrix} $$

Answer 1

Well there are three things as i see it with this question :

1) In your derivation when your taking the variance of the terms inside Rho should get squared and you should end up with the expression .. i didnt consider the auto covariance earlier ..sorry about that

$$ Var(\overline{x}) = \frac{\sigma_{\varepsilon}^2}{N} \frac{1}{1 - \rho^2} + \sum\limits_{t=0}^{N-1}\sum\limits_{t\neq j}^{N-1}\frac{\sigma_{\varepsilon}^2}{N^2} \frac{1}{1 - \rho^2}\rho^{|j-t|}$$

2) In your code you have calculated the variance of xbar ... for the standard error the code should include taking the sqrt of the answer given

3) You have assumed that the white noise has been generated from a (0,1) distribution when in fact the white noise only has to have constant variance .. i dont know what values of the constant variance R uses to generate the time series ... perhaps you could check on that ..

Hope this helps you :)

Answer 2

Well actually when you take the following

\begin{align*} Var(\overline{x}) &= Var\left(\frac{1}{N} \sum\limits_{t=0}^{N-1} x_t\right) \\ \end{align*}

It is easier to derive an implicit value rather than an explicit value in this case..your answer and mine are the same ..it's just that yours is a bit more difficult to handle because of the expansion of rho's..but some algebraic manipulation should be able to do the trick i guess ....I derived the answer as follows

since $ Var(\overline{x})$ is a linear combination..... \begin{align*} Var(\overline{x}) &= \frac{1}{N^2}Cov\left(\sum\limits_{t=0}^{N-1}\sum\limits_{j=0}^{N-1} x_t x_j\right) \\ \end{align*} \begin{align*} Var(\overline{x}) &= \frac{1}{N^2}\sum\limits_{t=0}^{N-1}Var\left( x_t \right) + \frac{1}{N^2}\sum\limits_{t=0}^{N-1}\sum\limits_{j \ne t}^{N-1}Cov\left( x_t x_j\right) \\ \end{align*}

Now for an AR(1) process $Var(x_t) = \frac{{\sigma_{\varepsilon}}^2}{1 - \rho^2} $ and $Cov(x_tx_j) = \frac{{\sigma_{\varepsilon}}^2}{1 - \rho^2}\rho^{|j-t|} $....

Substituting in the above gives the required equation... hope this answers your question :)

Answer 3

This is the R code btw ..

nrMCS <- 10000
N <- 100
pers <- 0.9
means <- numeric(nrMCS)


for (i in 1:nrMCS) {
  means[i] <- mean(arima.sim(model=list(ar=c(pers)), n = N,mean=0,sd=1))
}
#Simulation answer
ans1 <-sd(means)

#This should be the standard error according to the given formula
cov <- 0
for(i in 1:N){
for(j in 1:N){
cov <- cov +(1/((N^2)*(1-pers^2)))*pers^abs(j-i)
 }
}

ans2 <- sqrt(cov)

Answer 4

Don't know if it qualifies as a formal answer, but on the simulation side, standard error of a means estimator is defined as est(sd(means))/sqrt(N), which would give:

> .9459876/sqrt(100)
[1] 0.09459876

Not sure why you were using sd(means) and calling it standard error (if I understood the code comment right). It would make more sense to call the value SE(u) rather than Var(u) in the derivation as well, since I think that's what you intended?

The variance of an AR(1) process is the variance of the error term divided by (1-phi^2), where you had N*(1-phi) (the N wouldn't be there if it was just variance).

.

I'll have to dig deeper to try to find a derivation of that.

varianceAR1 simple AR(1) derivation on p. 36