[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*}[/step]

[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: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]