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} $$
Best Answer
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 :)