[proofplan]
We first center the observable and reduce the statement to a [central limit theorem](/theorems/521) for a stationary real-valued strongly mixing sequence. Davydov's covariance inequality gives an exponentially summable bound on the covariance sequence because the observable has a finite $(2+\delta)$-moment and the mixing coefficients decay geometrically. This proves absolute convergence of the long-run covariance series and identifies it as the limit of the normalized variances of partial sums. Finally, the Ibragimov [central limit theorem](/theorems/1848) for strongly mixing sequences applies under the same summability condition and yields the Gaussian limit with that variance.
[/proofplan]
[step:Center the observable and record the moment bound]
Define the centered observable $h:E\to\mathbb R$ by
\begin{align*}
h(x):=f(x)-\pi f.
\end{align*}
Since $f\in L^{2+\delta}(E,\mathcal E,\pi)$ and $\pi$ is a probability measure, [Jensen's inequality](/theorems/9) gives $|\pi f|<\infty$. Moreover,
\begin{align*}
|h(x)|^{2+\delta}\le 2^{1+\delta}\left(|f(x)|^{2+\delta}+|\pi f|^{2+\delta}\right)
\end{align*}
for every $x\in E$, hence
\begin{align*}
\int_E |h(x)|^{2+\delta}\,d\pi(x)<\infty.
\end{align*}
For each $n\in\mathbb N_0$, define the real-valued [random variable](/page/Random%20Variable) $Y_n:\Omega\to\mathbb R$ by
\begin{align*}
Y_n(\omega):=h(X_n(\omega)) \quad \text{for every } \omega\in\Omega.
\end{align*}
Because $(X_n)_{n\ge0}$ is stationary with marginal law $\pi$, the sequence $(Y_n)_{n\ge0}$ is stationary, $\mathbb E[Y_0]=0$, and
\begin{align*}
\mathbb E[|Y_0|^{2+\delta}]=\int_E |h(x)|^{2+\delta}\,d\pi(x)<\infty.
\end{align*}
Also,
\begin{align*}
\sqrt{N}\left(\frac{1}{N}\sum_{n=0}^{N-1}f(X_n)-\pi f\right)
=
\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}Y_n.
\end{align*}
Thus it is enough to prove the central limit theorem for the centered stationary sequence $(Y_n)_{n\ge0}$.
[/step]
[step:Bound the covariance sequence by Davydov's inequality]
For each $k\in\mathbb N_0$, define
\begin{align*}
\gamma_k:=\operatorname{Cov}(Y_0,Y_k)=\mathbb E[Y_0Y_k].
\end{align*}
Let $p:=2+\delta$ and let $q:=p$. Define $s>1$ by
\begin{align*}
\frac{1}{p}+\frac{1}{q}+\frac{1}{s}=1.
\end{align*}
Since $p=q=2+\delta$, this gives
\begin{align*}
\frac{1}{s}=\frac{\delta}{2+\delta}.
\end{align*}
Davydov's covariance inequality for strongly mixing sequences states that, for random variables $U$ and $V$ measurable with respect to two sigma-algebras separated by mixing coefficient $\alpha$, with $U\in L^p$, $V\in L^q$, and $1/p+1/q+1/s=1$, one has
\begin{align*}
|\operatorname{Cov}(U,V)|\le 8\alpha^{1/s}\|U\|_{L^p(\mathbb P)}\|V\|_{L^q(\mathbb P)}.
\end{align*}
This is a standard external result not yet linked in the wiki: Davydov covariance inequality.
Apply this inequality with $U=Y_0$ and $V=Y_k$. The random variable $Y_0$ is $\sigma(X_0)$-measurable, and $Y_k$ is $\sigma(X_k,X_{k+1},\dots)$-measurable. Hence the relevant mixing coefficient is bounded by $\alpha(k)$. By stationarity,
\begin{align*}
\|Y_k\|_{L^{2+\delta}(\mathbb P)}=\|Y_0\|_{L^{2+\delta}(\mathbb P)}.
\end{align*}
Therefore, for every $k\in\mathbb N_0$,
\begin{align*}
|\gamma_k|\le 8\alpha(k)^{\delta/(2+\delta)}\|Y_0\|_{L^{2+\delta}(\mathbb P)}^2.
\end{align*}
[guided]
The goal of this step is to turn the probabilistic mixing assumption into a numerical estimate on the covariance sequence. Define, for each $k\in\mathbb N_0$,
\begin{align*}
\gamma_k:=\operatorname{Cov}(Y_0,Y_k)=\mathbb E[Y_0Y_k],
\end{align*}
where the second equality uses $\mathbb E[Y_0]=\mathbb E[Y_k]=0$.
We use Davydov's covariance inequality. This standard external result, not yet linked in the wiki, says the following: if $U$ is measurable with respect to one sigma-algebra, $V$ is measurable with respect to another sigma-algebra, the two sigma-algebras have strong-mixing coefficient at most $\alpha$, and $U\in L^p$, $V\in L^q$ with
\begin{align*}
\frac{1}{p}+\frac{1}{q}+\frac{1}{s}=1,
\end{align*}
then
\begin{align*}
|\operatorname{Cov}(U,V)|\le 8\alpha^{1/s}\|U\|_{L^p(\mathbb P)}\|V\|_{L^q(\mathbb P)}.
\end{align*}
We choose $p=q=2+\delta$. Then the remaining exponent $s$ is determined by
\begin{align*}
\frac{1}{s}=1-\frac{1}{2+\delta}-\frac{1}{2+\delta}=\frac{\delta}{2+\delta}.
\end{align*}
This is exactly the exponent that will be summable under geometric mixing.
Now verify the hypotheses of Davydov's inequality. The random variable $Y_0=h(X_0)$ is measurable with respect to $\sigma(X_0)$. The random variable $Y_k=h(X_k)$ is measurable with respect to $\sigma(X_k,X_{k+1},\dots)$ because $\sigma(X_k)\subset\sigma(X_k,X_{k+1},\dots)$. By the definition of $\alpha(k)$, the mixing coefficient between these two sigma-algebras is at most $\alpha(k)$. The integrability condition also holds because the preceding step proved $Y_0\in L^{2+\delta}(\mathbb P)$, and stationarity gives $Y_k\in L^{2+\delta}(\mathbb P)$ with the same norm.
Applying Davydov's inequality with $U=Y_0$ and $V=Y_k$ gives
\begin{align*}
|\gamma_k|
=
|\operatorname{Cov}(Y_0,Y_k)|
\le
8\alpha(k)^{\delta/(2+\delta)}\|Y_0\|_{L^{2+\delta}(\mathbb P)}\|Y_k\|_{L^{2+\delta}(\mathbb P)}.
\end{align*}
By stationarity,
\begin{align*}
\|Y_k\|_{L^{2+\delta}(\mathbb P)}=\|Y_0\|_{L^{2+\delta}(\mathbb P)}.
\end{align*}
Hence, for every $k\in\mathbb N_0$,
\begin{align*}
|\gamma_k|\le 8\alpha(k)^{\delta/(2+\delta)}\|Y_0\|_{L^{2+\delta}(\mathbb P)}^2.
\end{align*}
[/guided]
[/step]
[step:Use exponential mixing to prove absolute convergence of the covariance series]
Define
\begin{align*}
\beta:=\frac{\delta}{2+\delta}.
\end{align*}
Then $\beta\in(0,1)$, and the exponential mixing hypothesis gives
\begin{align*}
\alpha(k)^\beta\le C^\beta r^{\beta k}
\end{align*}
for every $k\in\mathbb N_0$. Since $r^\beta\in(0,1)$,
\begin{align*}
\sum_{k=1}^{\infty}\alpha(k)^\beta\le C^\beta\sum_{k=1}^{\infty}r^{\beta k}<\infty.
\end{align*}
Combining this summability with the covariance estimate gives
\begin{align*}
\sum_{k=1}^{\infty}|\gamma_k|<\infty.
\end{align*}
Since centering by constants does not change covariance, $\gamma_k=\operatorname{Cov}_{\mathbb P}(f(X_0),f(X_k))$. Therefore the covariance series defining $\sigma_f^2$ converges absolutely.
[/step]
[step:Identify the long-run variance as a limit of normalized variances]
For each $N\in\mathbb N$, define the partial sum random variable $S_N:\Omega\to\mathbb R$ by
\begin{align*}
S_N(\omega):=\sum_{n=0}^{N-1}Y_n(\omega) \quad \text{for every } \omega\in\Omega.
\end{align*}
By stationarity and centering,
\begin{align*}
\operatorname{Var}(S_N)=\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}\operatorname{Cov}(Y_i,Y_j).
\end{align*}
For $i,j\in\{0,\dots,N-1\}$, stationarity gives $\operatorname{Cov}(Y_i,Y_j)=\gamma_{|i-j|}$. Counting the pairs with $|i-j|=k$ gives
\begin{align*}
\frac{1}{N}\operatorname{Var}(S_N)
=
\gamma_0+2\sum_{k=1}^{N-1}\left(1-\frac{k}{N}\right)\gamma_k.
\end{align*}
Because $\sum_{k=1}^{\infty}|\gamma_k|<\infty$, dominated convergence for series yields
\begin{align*}
\lim_{N\to\infty}\frac{1}{N}\operatorname{Var}(S_N)
=
\gamma_0+2\sum_{k=1}^{\infty}\gamma_k
=
\sigma_f^2.
\end{align*}
Each quantity $\operatorname{Var}(S_N)/N$ is nonnegative, so its limit satisfies $\sigma_f^2\ge0$.
[/step]
[step:Apply the strongly mixing central limit theorem]
The Ibragimov central limit theorem for strongly mixing stationary sequences is a standard external result not yet linked in the wiki. In the form needed here, it states that if $(Z_n)_{n\ge0}$ is a stationary centered real-valued sequence, $Z_0\in L^{2+\delta}(\mathbb P)$ for some $\delta>0$, and
\begin{align*}
\sum_{k=1}^{\infty}\alpha_Z(k)^{\delta/(2+\delta)}<\infty,
\end{align*}
where $\alpha_Z(k)$ denotes the usual strong-mixing coefficients
\begin{align*}
\alpha_Z(k):=\sup_{m\in\mathbb N_0}\sup_{A\in\sigma(Z_0,\dots,Z_m),\,B\in\sigma(Z_{m+k},Z_{m+k+1},\dots)}|\mathbb P(A\cap B)-\mathbb P(A)\mathbb P(B)|,
\end{align*}
then
\begin{align*}
\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}Z_n\xrightarrow{d}\mathcal N(0,\sigma^2),
\end{align*}
where
\begin{align*}
\sigma^2=\operatorname{Var}(Z_0)+2\sum_{k=1}^{\infty}\operatorname{Cov}(Z_0,Z_k).
\end{align*}
This version permits the degenerate case $\sigma^2=0$, in which $\mathcal N(0,0)$ is the point mass at $0$.
We apply this theorem with $Z_n=Y_n$. The sequence $(Y_n)_{n\ge0}$ is stationary and centered, and $Y_0\in L^{2+\delta}(\mathbb P)$. It remains to verify the mixing summability condition. For each $k\in\mathbb N_0$, define $\alpha_Y(k)$ by the displayed formula with $Z=Y$. Fix $m\in\mathbb N_0$, $A\in\sigma(Y_0,\dots,Y_m)$, and $B\in\sigma(Y_{m+k},Y_{m+k+1},\dots)$. Since $Y_j=h(X_j)$ for every $j\in\mathbb N_0$, measurability of composition gives
\begin{align*}
\sigma(Y_0,\dots,Y_m)\subset\sigma(X_0,\dots,X_m)
\end{align*}
and
\begin{align*}
\sigma(Y_{m+k},Y_{m+k+1},\dots)\subset\sigma(X_{m+k},X_{m+k+1},\dots).
\end{align*}
Thus the defining supremum for $\alpha(k)$ bounds the particular pair $(A,B)$, so
\begin{align*}
|\mathbb P(A\cap B)-\mathbb P(A)\mathbb P(B)|\le \alpha(k).
\end{align*}
Taking the supremum first over $A$ and $B$, and then over $m\in\mathbb N_0$, gives
\begin{align*}
\alpha_Y(k)\le \alpha(k).
\end{align*}
Therefore,
\begin{align*}
\sum_{k=1}^{\infty}\alpha_Y(k)^{\delta/(2+\delta)}
\le
\sum_{k=1}^{\infty}\alpha(k)^{\delta/(2+\delta)}
<\infty.
\end{align*}
All hypotheses of the Ibragimov central limit theorem are verified, and hence
\begin{align*}
\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}Y_n\xrightarrow{d}\mathcal N(0,\sigma_f^2).
\end{align*}
Using the identity from the centering step,
\begin{align*}
\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}Y_n
=
\sqrt{N}\left(\frac{1}{N}\sum_{n=0}^{N-1}f(X_n)-\pi f\right),
\end{align*}
which is exactly the asserted convergence.
[/step]
