Monte Carlo Standard Errors for MCMC Estimates

Why MCMC Estimates Need a Dependence-Aware Standard Error

An MCMC estimate of a target expectation is an average of values along a Markov chain. Unlike independent draws, successive chain states are usually dependent. That dependence changes the sampling variance of the average: a chain that moves slowly through the state space can provide much less information than the same number of independent draws. A Monte Carlo standard error (MCSE) estimates the resulting uncertainty in the reported average.

Let \((X_n)_{n\ge0}\) be a Markov chain on a finite, countable, or general measurable state space \(S\), on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). Its filtration is \(\mathcal{F}_n=\sigma(X_0,\ldots,X_n)\). Suppose the chain has stationary distribution \(\pi\), and for now assume it is started in stationarity, so \(X_n\sim\pi\) for every \(n\). For a function \(g:S\to\mathbb{R}\) with \(\mathbb{E}_{\pi}[g(X_0)^2]<\infty\), the target is \(\theta=\mathbb{E}_{\pi}[g(X_0)]\), estimated from \(n\) draws by \(\bar g_n=n^{-1}\sum_{t=1}^n g(X_t)\).

Definition: Write \(Y_t=g(X_t)-\theta\), and define the lag-\(k\) autocovariance by \(\gamma_k=\operatorname{Cov}_{\pi}(g(X_0),g(X_k))=\mathbb{E}_{\pi}[Y_0Y_k]\). If \(\sum_{k=1}^{\infty}|\gamma_k|<\infty\), the long-run variance is \(\sigma^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k\). When a central limit theorem holds, \(\sqrt{n}(\bar g_n-\theta)\) converges in distribution to a mean-zero normal variable with variance \(\sigma^2\), and the asymptotic MCSE is \(\sigma/\sqrt{n}\).

The variance \(\gamma_0\) is the one-draw variance. The additional covariance terms account for serial dependence. An independent-draw calculation drops them and uses an estimate of \(\sqrt{\gamma_0/n}\); that is appropriate for independent draws, but generally not for MCMC output.

How Dependence Changes the Variance

Result: If the chain is stationary and \(\sum_{k=1}^{\infty}|\gamma_k|<\infty\), then \(n\operatorname{Var}(\bar g_n)\to\sigma^2\), where \(\sigma^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k\).

Proof. By stationarity, \(\operatorname{Cov}(Y_s,Y_t)=\gamma_{|t-s|}\). Expanding the variance of the sum gives \(n\) terms with lag zero and \(2(n-k)\) terms with each lag \(k\), for \(1\le k<n\). Therefore

n\operatorname{Var}(\bar g_n) =\gamma_0+2\sum_{k=1}^{n-1}\left(1-\frac{k}{n}\right)\gamma_k.

For each fixed \(k\), the multiplier \(1-k/n\) tends to \(1\) as \(n\to\infty\), and its absolute value is at most \(1\) for \(1\le k<n\). Extend the summand by zero for \(k\ge n\). Its absolute value is then bounded by \(|\gamma_k|\), whose sum is finite by assumption. The dominated convergence theorem for series gives convergence of the sum to \(\sum_{k=1}^{\infty}\gamma_k\), proving the result. \(\square\)

This variance identity does not itself assert a central limit theorem. For a finite-state irreducible, aperiodic chain and bounded \(g\), standard finite-state Markov chain results give the relevant central limit theorem and summable covariances. For general state spaces, the central limit theorem requires additional mixing or drift conditions; it should not be inferred from a finite run alone. The estimation methods below estimate long-run variance, not whether all the assumptions behind an asymptotic approximation hold.

Worked Example: Exact Long-Run Variance in a Two-State Chain

Consider a chain on \(S=\{0,1\}\) with transition matrix \(P\), and estimate the stationary expectation of \(g(i)=i\). The matrix is row-stochastic: the first row sums to \(\tfrac34+\tfrac14=1\), and the second sums to \(\tfrac14+\tfrac34=1\).

P=\begin{bmatrix} \tfrac34&\tfrac14\\ \tfrac14&\tfrac34 \end{bmatrix}.

Solving \(\pi P=\pi\) and \(\pi_0+\pi_1=1\) gives \(\pi=(\tfrac12,\tfrac12)\): each coordinate of \(\pi P\) is \(\tfrac12\). Thus \(\theta=\tfrac12\), and \(Y_t=X_t-\tfrac12\) takes values \(-\tfrac12\) and \(\tfrac12\). In particular, \(\gamma_0=\mathbb{E}_{\pi}[Y_0^2]=\tfrac14\).

Given \(X_t\), the conditional mean of \(Y_{t+1}\) is \(\tfrac12Y_t\). Indeed, if \(X_t=0\), then \(\mathbb{E}[Y_{t+1}\mid\mathcal{F}_t]=\tfrac14(\tfrac12)+\tfrac34(-\tfrac12)=-\tfrac14=\tfrac12Y_t\); if \(X_t=1\), it is \(\tfrac34(\tfrac12)+\tfrac14(-\tfrac12)=\tfrac14=\tfrac12Y_t\). Iterating this conditional-mean relation gives \(\mathbb{E}[Y_k\mid\mathcal{F}_0]=(\tfrac12)^kY_0\). Hence \(\gamma_k=(\tfrac12)^k\mathbb{E}[Y_0^2]=\tfrac14(\tfrac12)^k\).

Consequently, \(\sum_{k=1}^{\infty}\gamma_k=\tfrac14\sum_{k=1}^{\infty}(\tfrac12)^k=\tfrac14\), so \(\sigma^2=\tfrac14+2(\tfrac14)=\tfrac34\). The independent-draw variance per observation would be only \(\gamma_0=\tfrac14\). In this example the long-run variance is three times the independent-draw variance, so the asymptotic MCSE is \(\sqrt{3}\) times as large. Check: Substituting the geometric sum into the definition gives \(\tfrac14+2\cdot\tfrac14=\tfrac34\), as claimed.

Batch Means

Batch means estimate dependence by averaging adjacent observations into blocks. Divide \(n=ab\) observations into \(a\) nonoverlapping batches of size \(b\), and let \(\bar g_j^{(b)}\) be the average in batch \(j\). The batch-means estimate of long-run variance is \(b\) times the sample variance of the batch averages. The multiplication by \(b\) compensates for averaging \(b\) successive, correlated observations: when batches are long enough, \(b\operatorname{Var}(\bar g_j^{(b)})\) approximates the long-run variance.

Definition: For equal batch size \(b\), \(a\ge2\) batches, and overall average \(\bar g_n\), define \(\widehat{\sigma}^2_{\mathrm{BM}}=\frac{b}{a-1}\sum_{j=1}^{a}(\bar g_j^{(b)}-\bar g_n)^2\). The batch-means MCSE is \(\sqrt{\widehat{\sigma}^2_{\mathrm{BM}}/n}\).

Batch size matters. Very short batches may fail to capture substantial dependence; very long batches leave too few batch averages to estimate their variance reliably. A common asymptotic regime takes both \(b\) and \(a=n/b\) to increase. Batch means are not guaranteed to be consistent for arbitrary choices of \(b\), chain behavior, or initial distribution. Consistency results for batch means require regularity and growth conditions; their proofs are beyond this tutorial.

Worked Example: Batch Means Compared with the Independent-Draw Formula

Suppose the recorded values of \(g(X_t)\) are \(0,0,1,1,1,0,0,1,0,0,1,1\). There are \(n=12\) observations, with seven ones, so \(\bar g_n=7/12\). Use four batches of size \(b=3\). The batch averages are \(\tfrac13,1,\tfrac13,\tfrac23\), whose deviations from \(7/12\) are \(-\tfrac14,\tfrac5{12},-\tfrac14,\tfrac1{12}\).

The sum of squared deviations is \(\tfrac1{16}+\tfrac{25}{144}+\tfrac1{16}+\tfrac1{144}=\tfrac{44}{144}=\tfrac{11}{36}\). Therefore \(\widehat{\sigma}^2_{\mathrm{BM}}=\tfrac{3}{4-1}\cdot\tfrac{11}{36}=\tfrac{11}{36}\), and the estimated MCSE is \(\sqrt{(\tfrac{11}{36})/12}=\sqrt{11}/(12\sqrt{3})\).

For comparison, the sample variance of the twelve values is \(\tfrac{12}{11}\cdot\frac{7}{12}\cdot\frac{5}{12}=\tfrac{35}{132}\). The independent-draw standard error estimate is therefore \(\sqrt{(\tfrac{35}{132})/12}=\sqrt{35/1584}\). The two estimates differ because they use different information: the independent-draw formula treats the observations as uncorrelated, whereas batch means uses their ordering. Check: The sample variance calculation also follows by summing squared deviations from \(7/12\): the sum is \(7(\tfrac5{12})^2+5(-\tfrac7{12})^2=\tfrac{35}{12}\), and dividing by \(n-1=11\) gives \(\tfrac{35}{132}\).

Spectral-Variance Estimation

A spectral-variance estimate uses sample autocovariances at several lags. Given observations \(g_t=g(X_t)\), define \(\bar g_n=n^{-1}\sum_{t=1}^n g_t\) and, for \(0\le k<n\), use the biased sample autocovariance \(\widehat{\gamma}_k=n^{-1}\sum_{t=1}^{n-k}(g_t-\bar g_n)(g_{t+k}-\bar g_n)\). The divisor is \(n\) at every lag. A lag window gives less weight to distant sample autocovariances, which are noisier.

Definition: With Bartlett lag window and truncation lag \(L<n\), the spectral-variance estimate is \(\widehat{\sigma}^2_{\mathrm{SV}}=\widehat{\gamma}_0+2\sum_{k=1}^{L}(1-\frac{k}{L+1})\widehat{\gamma}_k\). Its MCSE estimate is \(\sqrt{\widehat{\sigma}^2_{\mathrm{SV}}/n}\) when the estimated variance is nonnegative.

This is a truncated, weighted estimate of \(\gamma_0+2\sum_{k\ge1}\gamma_k\). The Bartlett weight at lag \(k\) is \(1-k/(L+1)\); thus the last included lag receives weight \(1/(L+1)\), not zero. A fixed, small \(L\) estimates only a finite part of the covariance sum and need not estimate the long-run variance well. Consistency requires an increasing lag choice together with conditions on dependence and the rate of increase. Standard spectral-variance consistency results establish this under suitable assumptions; their proofs are beyond this tutorial. Neither a spectral estimate nor a batch-means estimate is automatically reliable just because it can be computed.

Worked Example: A Bartlett Spectral Estimate from the Same Run

Use the twelve binary observations in the preceding example and choose \(L=1\). Centering by \(7/12\), each one contributes \(5/12\) and each zero contributes \(-7/12\). Thus \(\widehat{\gamma}_0=[7(25/144)+5(49/144)]/12=(420/144)/12=35/144\).

For lag one, the eleven adjacent pairs contain three \(00\) pairs, three \(11\) pairs, and five mixed pairs. Their centered products are respectively \(49/144\), \(25/144\), and \(-35/144\). Therefore \(\widehat{\gamma}_1=[3(49/144)+3(25/144)+5(-35/144)]/12=(47/144)/12=47/1728\). With \(L=1\), the Bartlett weight on lag one is \(1-\tfrac12=\tfrac12\), so

\widehat{\sigma}^2_{\mathrm{SV}} =\frac{35}{144}+2\left(\frac12\right)\frac{47}{1728} =\frac{420+47}{1728} =\frac{467}{1728}.

The resulting MCSE is \(\sqrt{(\tfrac{467}{1728})/12}=\sqrt{467}/144\). Check: With denominator \(n=12\), the lag-zero calculation is \((420/144)/12=35/144\), and the lag-one calculation is \((147+75-175)/144/12=47/1728\). Substitution gives the displayed spectral estimate. Its difference from the batch-means estimate reflects their different ways of summarizing this short run, not a contradiction.

Interpreting and Reporting an MCSE

An MCSE describes Monte Carlo uncertainty in \(\bar g_n\), the estimate of \(\theta\), under the chain and approximation being used. It does not measure whether the model is appropriate, whether the target distribution represents the scientific question, or whether the chain has reached the intended target. If \(\widehat{\sigma}^2\) estimates the long-run variance, the corresponding standard error is \(\sqrt{\widehat{\sigma}^2/n}\). Increasing \(n\) reduces this asymptotic uncertainty at the \(1/\sqrt{n}\) rate only when the chain has suitable limiting behavior and the long-run variance is finite.

A useful comparison is the independent-draw estimate \(\sqrt{s_g^2/n}\), where \(s_g^2=(n-1)^{-1}\sum_{t=1}^n(g_t-\bar g_n)^2\). Comparing it with a dependence-aware MCSE can reveal the impact of correlation, but the independent-draw value is not a substitute for estimating that dependence. For the two-state chain above, \(\sigma^2/\gamma_0=3\); this ratio is sometimes interpreted as an asymptotic variance inflation factor. It is not a general guarantee that every chain has a positive, simple effective-sample-size interpretation.

Common mistake. In the twelve-value example, using only the sample variance and reporting \(\sqrt{35/1584}\) as the MCMC standard error implicitly sets all nonzero-lag covariances to zero. That is the independent-draw calculation, not a dependence-aware one. The correction is to estimate long-run variance from ordered output—for example, the batch-means estimate \(\tfrac{11}{36}\) or, with the specified Bartlett lag, the spectral estimate \(\tfrac{467}{1728}\)—then divide by \(n=12\) and take the square root. These two corrections need not agree exactly on a short run; the method, batch size or lag choice, and estimated MCSE should be reported.

Check Your Understanding

For each question, distinguish the one-draw variance from the long-run variance.

State the definition of \(\sigma^2\) in terms of the autocovariances \(\gamma_k\).
For stationary output with summable absolute autocovariances, explain why \(n\operatorname{Var}(\bar g_n)\) converges to the long-run variance.
If \(n=20\) observations are divided into \(a=4\) equal batches, what is the batch size, and what quantity is multiplied by that size in the batch-means estimate?
For a Bartlett spectral estimate with \(L=2\), calculate the weight on lag one and on lag two.
Why is \(\sqrt{s_g^2/n}\) not generally a valid MCMC standard error, even if the chain has a stationary distribution?