$ \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\ol}[1]{\bar{#1}} \newcommand{\Cov}{\on{Cov}} \newcommand{\Var}{\on{Var}} $ Problem Statement: Suppose that the time series data $\{x_i:i=1,\dots,N\}$ are generated from an $\on{AR}(1)$ process with parameter $\alpha.$ Show that the variance of their sample mean can be approximated as follows for large $N,$ $$\Var\!\left(\ol{X}\right)\approx\frac{\sigma^2}{N}\! \left(\frac{1+\alpha}{1-\alpha}\right)\!.$$
This is Problem 4.1 in The Analysis of Time Series: An Introduction with R, 7th Ed., by Chris Chatfield and Haipeng Xing.
My Work So Far: So the $\on{AR}(1)$ process with parameter $\alpha$ we write down as $$X_t=\alpha X_{t-1}+Z_t.$$ The sample mean is, of course, given by $$\ol{X}=\frac1N\sum_{i=1}^Nx_i.$$ Now the variance of this sample mean we cannot compute simply by adding the variances of the individual terms, because they are dependent. First, let's start at $X_1$ and ratchet it up: $$X_k=\alpha^{k-1}X_1+\sum_{m=0}^{k-2}\alpha^m Z_{k-m},\quad k=2,\dots,N.$$
Now then, we wish to compute: \begin{align*} \Var\!\left(\ol{X}\right) &=\Cov\!\left(\frac1N\sum_{i=1}^NX_i,\frac1N\sum_{j=1}^NX_j\right)\\ &=\frac{1}{N^2}\!\left[\Cov(X_1,X_1) +2\sum_{j=2}^N\Cov\!\left(X_1,X_j\right) +\sum_{i=2}^N\sum_{j=2}^N\Cov(X_i,X_j)\right]\\ &=\frac{1}{N^2}\,\Bigg[\frac{\sigma_Z^2}{1-\alpha^2}+ 2\sum_{j=2}^N\Cov\!\left(X_1, \alpha^{j-1}X_1+\sum_{m=0}^{j-2}\alpha^m Z_{j-m}\right)\\ &\qquad\qquad+ \sum_{i=2}^N\sum_{j=2}^N \Cov\!\left( \alpha^{i-1}X_1+\sum_{m=0}^{i-2}\alpha^m Z_{i-m}, \alpha^{j-1}X_1+\sum_{m=0}^{j-2}\alpha^m Z_{j-m}\right)\Bigg]. \end{align*} The $Z_t$'s are independent. So this computation becomes \begin{align*} \Var\!\left(\ol{X}\right) &=\frac{1}{N^2}\, \Bigg[\frac{\sigma_Z^2}{1-\alpha^2}+ \frac{2\sigma_Z^2}{1-\alpha^2}\sum_{j=2}^N\alpha^{j-1} +2\sum_{j=2}^N\sum_{m=0}^{j-2}\alpha^m \Cov(X_1,Z_{j-m})\\ &+ \sum_{i=2}^N\sum_{j=2}^N \Bigg\{\frac{\alpha^{i+j-2}\sigma_Z^2}{1-\alpha^2}+ \alpha^{i-1}\sum_{m=0}^{j-2}\alpha^m\Cov(X_1,Z_{j-m}) +\alpha^{j-1}\sum_{m=0}^{i-2}\alpha^m\Cov(X_1,Z_{i-m})\\ &\qquad\qquad\qquad+\sum_{m=0}^{i-2}\sum_{n=0}^{j-2}\alpha^{m+n} \Cov(Z_{i-m},Z_{j-n})\Bigg\} \Bigg]. \end{align*} Now we work on the expression: \begin{align*} \sum_{m=0}^{i-2}\sum_{n=0}^{j-2}\alpha^{m+n}\Cov(Z_{i-m},Z_{j-n}) &=\sum_{p=2}^i\sum_{q=2}^j\alpha^{i+j-(p+q)}\Cov(Z_p,Z_q), \end{align*} where $p=i-m$ and $q=j-n.$ This expression simplifies to \begin{align*} \sum_{p=2}^{\min(i,j)}\alpha^{i+j-2p}\sigma_Z^2 &=\alpha^{i+j}\cdot\sigma_Z^2\sum_{p=2}^{\min(i,j)}\alpha^{-2p}\\ &=\alpha^{i+j}\cdot\sigma_Z^2\,\frac{1-\alpha^{2-2\min(i,j)}}{\alpha^2(\alpha^2-1)}. \end{align*} So much for the last term. For the two middle terms with covariances of the form $\Cov(X_1,Z_{j-m}),$ note that $2\le j-m\le j$ as well as $2\le i-m\le i.$ We can certainly assume $X_1$ is independent of these higher $Z_t$'s, so those covariances are zero. Note also that $$\sum_{j=2}^N\alpha^{j-1}=\frac{\alpha^N-\alpha}{\alpha-1}.$$ Hence, we have that \begin{align*} \Var\!\left(\ol{X}\right) &=\frac{\sigma_Z^2}{N^2}\left[\frac{1}{1-\alpha^2} +\frac{2}{1-\alpha^2}\cdot\frac{\alpha^N-\alpha}{\alpha-1} +\sum_{i=2}^N\sum_{j=2}^N \left\{ \frac{\alpha^{i+j-2}}{1-\alpha^2} +\alpha^{i+j}\cdot \frac{1-\alpha^{2-2\min(i,j)}}{\alpha^2(\alpha^2-1)} \right\} \right]. \end{align*} We keep moving forward: \begin{align*} \sum_{i=2}^N\sum_{j=2}^N\frac{\alpha^{i+j-2}\,\sigma_Z^2}{1-\alpha^2} &=\frac{\sigma_Z^2}{\alpha^2(1-\alpha^2)}\sum_{i=2}^N\sum_{j=2}^N\alpha^{i+j}\\ &=\frac{\sigma_Z^2}{\alpha^2(1-\alpha^2)}\cdot \frac{\alpha^2(\alpha^N-\alpha)^2}{(\alpha-1)^2}\\ &=\frac{\sigma_Z^2(\alpha^N-\alpha)^2}{(1+\alpha)(1-\alpha)^3}. \end{align*} Now assuming $N>2,$ we have \begin{align*} \sum_{i=2}^N\sum_{j=2}^N\alpha^{i+j}\frac{1-\alpha^{2-2\min(i,j)}}{\alpha^2(\alpha^2-1)} &=\frac{1}{\alpha^2(\alpha+1)(\alpha-1)}\sum_{i=2}^N\sum_{j=2}^N \alpha^{i+j}\left(1-\alpha^{2-2\min(i,j)}\right)\\ &=\frac{\alpha^{2N}-2(\alpha+1)\alpha^N+\alpha(\alpha N+2)-N+1}{(\alpha-1)^3(\alpha+1)}. \end{align*} Now then, we can plug these results back into the expression for the variance, and take the limit as $N\to\infty.$ We do this in Mathematica, with the result that for large $N,$ $$\Var\!\left(\ol{X}\right) \approx\frac{\sigma_Z^2}{N(1-\alpha)^2}.$$ This is not the result hoped for, but it is very similar.
My Question: Where have I gone wrong?