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

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}$$

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 :)

• Very good points, thank you. However, I think your formula does not work. See my edit above. I think due to the fact that you have to square the terms in the bracket now, you get more terms: $(\rho^0 + \rho^1)^2$ is not $(\rho^0 + \rho^2)$, but $\rho^0 + 2\rho^0 \rho + \rho^2$. I don't see any easy way now to simplify things. Maybe you have clue? – Christoph_J Oct 20 '12 at 15:16
• Well i modified my earlier answer to incorporate the auto covariance and did a bit of amendments to your R code .... the output i got seemed reasonably close to the expected answer .... – Nuzhi Oct 20 '12 at 16:43
• I accepted the answer; it seems indeed to work. So thanks! Only problem: I have a headache today and don't quite understand how you adjusted for the covariance terms. Thinking about those rhos just makes my head spin...So it would be great if you could add a few thoughts how you got from my updated rho matrix to your formula. – Christoph_J Oct 20 '12 at 19:40

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 :)

• Perfect, thank you very much! Well deserved bounty ;-) – Christoph_J Oct 21 '12 at 7:27

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))
}
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)


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

• No, I think my way of computing the standard error should be correct. As pointed out here, the standard error is the standard deviation of the sampling distribution of a statistic. My statistic is the mean based on $N$ observations, so I compute the standard deviation of that statistic. Another way to see it: I already consider $\sqrt{N}$ by computing the mean with $N$, instead of 1, observation. – Christoph_J Oct 20 '12 at 15:22
• On the simulation. Sorry, I missed the looping of means. Agree N is taken into account there. – pat Oct 20 '12 at 17:00