Autocorrelation and Effective Sample Size

Why Correlation Changes the Information in a Chain

An MCMC run may contain thousands of draws, but adjacent values often resemble one another. Counting every recorded value as if it were an independent observation can therefore overstate how much information the run provides about a target quantity. Autocorrelation describes this dependence across time; integrated autocorrelation time summarizes its effect on the variance of a chain average.

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, let $(X_n)_{n\ge0}$ be a discrete-time Markov chain on a finite or countable state space $S$, and let $g:S\to\mathbb{R}$ be a scalar-valued function. Write $Y_n=g(X_n)$. For the definitions and exact variance identity below, assume the chain is started in a stationary distribution, so $(Y_n)$ is stationary, and assume $\mathbb{E}[Y_0^2]<\infty$. Stationarity means that the joint distribution of a finite collection of values is unchanged by shifting all their time indices by the same nonnegative amount. Write $\mu_Y=\mathbb{E}[Y_0]$, and suppose $\gamma_0=\operatorname{Var}(Y_0)>0$.

Definition: The lag-$k$ autocovariance is $\gamma_k=\operatorname{Cov}(Y_0,Y_k)$, for an integer $k\ge0$. The lag-$k$ autocorrelation is $\rho_k=\gamma_k/\gamma_0$. In particular, $\rho_0=1$. If the series converges, the integrated autocorrelation time is $\tau_{\mathrm{int}}=1+2\sum_{k=1}^{\infty}\rho_k$.

Autocorrelation is defined only when the variance in the denominator is positive and finite. If $Y_n$ is constant almost surely, then $\gamma_0=0$, so this definition does not apply. The formula for $\tau_{\mathrm{int}}$ requires convergence of the autocorrelation series; absolute summability, $\sum_{k\ge1}|\rho_k|<\infty$, is a sufficient condition. Positive autocorrelation makes $\tau_{\mathrm{int}}$ larger, while negative autocorrelation can make it smaller. It is not correct to assume all autocorrelations are positive.

How Integrated Autocorrelation Time Enters the Variance

For $N\ge1$, define the chain average $\bar Y_N=N^{-1}\sum_{n=1}^N Y_n$. The following identity is exact for every $N$ when the output is stationary. It is the reason integrated autocorrelation time provides a useful measure of the information in a run.

Theorem: If $(Y_n)$ is stationary with finite variance $\gamma_0$, then $$ \operatorname{Var}(\bar Y_N)=\frac{\gamma_0}{N}\left[1+2\sum_{k=1}^{N-1}\left(1-\frac{k}{N}\right)\rho_k\right]. $$ If, in addition, $\sum_{k\ge1}|\rho_k|<\infty$, then $N\operatorname{Var}(\bar Y_N)/\gamma_0\to\tau_{\mathrm{int}}$.

Proof. Set $Z_n=Y_n-\mu_Y$, so $\mathbb{E}[Z_n]=0$. Expanding the variance of the sum gives $\operatorname{Var}(\bar Y_N)=N^{-2}\sum_{r=1}^N\sum_{s=1}^N\operatorname{Cov}(Y_r,Y_s)$. Stationarity gives $\operatorname{Cov}(Y_r,Y_s)=\gamma_{|s-r|}$. There are $N$ pairs with $r=s$, and, for each $k\in\{1,\ldots,N-1\}$, there are $N-k$ pairs with $s-r=k$ and $N-k$ with $r-s=k$. Thus the double sum is $N\gamma_0+2\sum_{k=1}^{N-1}(N-k)\gamma_k$. Dividing by $N^2$, then using $\gamma_k=\gamma_0\rho_k$, proves the identity. For the limit, multiply the identity by $N/\gamma_0$. For each fixed $k$, the weight $1-k/N$ tends to $1$. To justify taking the limit of the sum, choose $M$ so that $\sum_{k>M}|\rho_k|$ is as small as desired. The first $M$ weighted terms converge individually, and the absolute value of the remaining weighted terms is at most $\sum_{k>M}|\rho_k|$, since $0\le1-k/N\le1$ for $k<N$. This proves convergence to $1+2\sum_{k\ge1}\rho_k=\tau_{\mathrm{int}}$. $\square$

When the assumptions for the limit hold, the variance of a long chain average is approximately $\gamma_0\tau_{\mathrm{int}}/N$. An independent sample of size $m$ from the same marginal distribution has average variance $\gamma_0/m$. Equating these two variance expressions suggests the effective sample size $N/\tau_{\mathrm{int}}$.

Definition: When $\tau_{\mathrm{int}}>0$, the variance-equivalent effective sample size for $N$ draws is $N_{\mathrm{eff}}=N/\tau_{\mathrm{int}}$. This is a comparison of variances for estimating the scalar mean, not a claim that the chain literally contains that many independent draws. Negative dependence can give $N_{\mathrm{eff}}>N$.

The limit in the theorem also shows that $\tau_{\mathrm{int}}\ge0$ under its hypotheses: it is the limit of the nonnegative quantities $N\operatorname{Var}(\bar Y_N)/\gamma_0$. It may equal zero, in which case the usual ratio $N/\tau_{\mathrm{int}}$ is not defined. Stationarity matters here: a chain having a stationary distribution does not by itself mean that a run started from an arbitrary initial state is stationary.

Worked Example: A Two-State Chain with Positive Dependence

Let $S=\{0,1\}$, and let the chain have transition matrix $P$, with rows and columns in the order $0,1$: $$ P=\begin{bmatrix}\tfrac34&\tfrac14\\[2pt]\tfrac14&\tfrac34\end{bmatrix}. $$ Each row sums to $1$: $\tfrac34+\tfrac14=1$ and $\tfrac14+\tfrac34=1$. The distribution $\pi=(\tfrac12,\tfrac12)$ is stationary, since $\pi P=(\tfrac12,\tfrac12)$. Start the chain with this distribution and take $Y_n=X_n$.

The mean is $\mu_Y=\tfrac12$, and the variance is $\gamma_0=\tfrac12(0-\tfrac12)^2+\tfrac12(1-\tfrac12)^2=\tfrac14$. Define $f(i)=i-\tfrac12$. For state $0$, the conditional mean of $f(X_1)$ is $\tfrac34(-\tfrac12)+\tfrac14(\tfrac12)=-\tfrac14=\tfrac12 f(0)$. For state $1$, it is $\tfrac14(-\tfrac12)+\tfrac34(\tfrac12)=\tfrac14=\tfrac12 f(1)$. Therefore, $Pf=\tfrac12f$, and repeated conditioning gives $\mathbb{E}[f(X_k)\mid X_0]=2^{-k}f(X_0)$. Consequently, $$ \gamma_k=\mathbb{E}[f(X_0)f(X_k)]=2^{-k}\mathbb{E}[f(X_0)^2]=\frac{2^{-k}}4,\qquad \rho_k=2^{-k}. $$

The geometric series gives $\sum_{k=1}^{\infty}2^{-k}=1$, so $\tau_{\mathrm{int}}=1+2(1)=3$. For $N=300$ stationary draws, the variance-equivalent effective sample size is $300/3=100$. As a check, the theorem’s finite-$N$ expression is $N\operatorname{Var}(\bar Y_N)/\gamma_0=1+2\sum_{k=1}^{N-1}(1-k/N)2^{-k}$, which is below $3$ and tends to $3$ as $N$ increases.

Worked Example: Independent Draws as a Baseline

Consider a chain on $S=\{0,1\}$ with transition matrix $$ P=\begin{bmatrix}\tfrac34&\tfrac14\\[2pt]\tfrac34&\tfrac14\end{bmatrix}. $$ Both rows sum to $1$, since $\tfrac34+\tfrac14=1$. Its stationary distribution is $\pi=(\tfrac34,\tfrac14)$, because multiplying this row vector by $P$ gives $(\tfrac34,\tfrac14)$. In fact, each next state is drawn from this distribution regardless of the current state. Start at stationarity and set $Y_n=X_n$.

The mean is $\mu_Y=\tfrac14$, and the variance is $\gamma_0=\tfrac34(0-\tfrac14)^2+\tfrac14(1-\tfrac14)^2=\tfrac3{16}$. For $k\ge1$, independence gives $\gamma_k=0$, so $\rho_k=0$ and $\tau_{\mathrm{int}}=1$. For $N=12$, the chain-average variance is $\gamma_0/N=(\tfrac3{16})/12=\tfrac1{64}$; the corresponding effective sample size is $12/1=12$. Here there is no reduction in information because the draws are independent. This baseline makes the two-state dependent-chain example concrete: dependence increased the variance-equivalent sample requirement by a factor of three.

Estimating Autocorrelation from a Recorded Sequence

In practice, the autocovariances are unknown, so they are estimated from the recorded scalar values $y_1,\ldots,y_N$. Let $\bar y=N^{-1}\sum_{t=1}^N y_t$. One common convention, used in the example below, is $$ \widehat{\gamma}_k=\frac1N\sum_{t=1}^{N-k}(y_t-\bar y)(y_{t+k}-\bar y),\qquad \widehat{\rho}_k=\frac{\widehat{\gamma}_k}{\widehat{\gamma}_0}, $$ where $\widehat{\gamma}_0=N^{-1}\sum_{t=1}^N(y_t-\bar y)^2$. This convention uses the same denominator $N$ at every lag; other conventions exist, so the convention should be stated when reporting a calculation. The ratio is defined only if the recorded values are not all identical.

Because only finitely many autocorrelations can be computed, a direct sum to infinity is impossible. A simple truncated estimate, for a chosen integer $K$ with $1\le K<N$, is $\widehat{\tau}_K=1+2\sum_{k=1}^K\widehat{\rho}_k$, with $\widehat{N}_{\mathrm{eff}}=N/\widehat{\tau}_K$ when $\widehat{\tau}_K>0$. A short run can make the estimated autocorrelations noisy, especially at larger lags. The choice of $K$ therefore affects the estimate; no finite sequence alone turns this truncated sum into an exact value of the population integrated autocorrelation time.

Worked Example: Estimating Effective Sample Size from Eight Values

Consider the sequence $y_1,\ldots,y_8=(0,1,0,1,0,1,0,1)$. Its mean is $\bar y=4/8=\tfrac12$. The centered values alternate between $-\tfrac12$ and $\tfrac12$, so $\widehat{\gamma}_0=8(\tfrac14)/8=\tfrac14$. Use the stated autocorrelation convention and truncate at $K=2$.

At lag $1$, all $7$ adjacent products equal $-\tfrac14$. Thus $\widehat{\gamma}_1=7(-\tfrac14)/8=-\tfrac7{32}$, and $\widehat{\rho}_1=(-\tfrac7{32})/(\tfrac14)=-\tfrac78$. At lag $2$, all $6$ products equal $\tfrac14$. Thus $\widehat{\gamma}_2=6(\tfrac14)/8=\tfrac3{16}$, and $\widehat{\rho}_2=(\tfrac3{16})/(\tfrac14)=\tfrac34$. Therefore, $$ \widehat{\tau}_2=1+2\left(-\frac78+\frac34\right)=1+2\left(-\frac18\right)=\frac34,\qquad \widehat{N}_{\mathrm{eff}}=\frac8{3/4}=\frac{32}{3}. $$

The estimated effective sample size exceeds eight because the first estimated autocorrelation is negative and outweighs the positive lag-two value in this particular truncated sum. This is an estimate from a very short, highly patterned sequence, not evidence that the sequence truly provides more information than eight independent observations. The calculation also illustrates why one should report the truncation rule and interpret a short-run estimate cautiously.

Interpretation and a Common Mistake

A frequent error is to use $N$ as the effective sample size simply because $N$ values were recorded. For the stationary two-state chain above, that would report $300$ effective draws. The autocorrelation calculation instead gives $\tau_{\mathrm{int}}=3$, and hence $N_{\mathrm{eff}}=300/3=100$. The unadjusted count ignores the positive dependence between nearby draws.

The opposite error is to treat a finite-sample estimate as a guaranteed population quantity. Estimated autocorrelations fluctuate, and a chosen truncation can produce an implausible or even nonpositive $\widehat{\tau}_K$. In that case $N/\widehat{\tau}_K$ is not a meaningful effective sample size under this formula. Reconsider the sequence and estimation procedure rather than interpreting the resulting ratio literally.

Effective sample size is specific to the scalar quantity being monitored: it concerns the chain output $Y_n=g(X_n)$, not the state vector in the abstract. Different functions of the same chain can have different autocorrelations and therefore different effective sample sizes. The variance identity gives the precise basis for the comparison, while the examples show why a nominal draw count alone does not describe the information in an MCMC estimate.

Check Your Understanding

Use the definitions and variance identity to answer the following questions.

For a stationary scalar process with $\rho_k=0$ for all $k\ge1$, what are $\tau_{\mathrm{int}}$ and $N_{\mathrm{eff}}$?
For the two-state chain in the first worked example, compute $\rho_3$ and explain its sign.
Why does the exact variance identity require stationarity of the scalar output?
Using $N=100$ and $\tau_{\mathrm{int}}=4$, find the variance-equivalent effective sample size.
What should you check before interpreting $N/\widehat{\tau}_K$ from a short sequence?

Autocorrelation and Effective Sample Size

What You'll Learn