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

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)

