How can we numerically compute the autocorrelation of a sample from a Markov chain generated by the Metropolis-Hastings algorithm?

Question

Let $(X_n)_{n\in\mathbb N_0}$ denote a $\mathbb R^d$-valued Markov chain generated by the Metropolis-Hastings algorithm. Suppose I've run the algorithm on a computer and obtained a sample $x_0,\ldots,x_n$.

How can I compute the (truncated) autocorrelation of this sample?

I'm a bit unsure which quantity the autocorrelation actually is. What I really want to do is computing $\operatorname{eff}_{\overline\theta}$ from the paper Efficient Metropolis Jumping Rules (equation $(1)$ on page 600).

Let $$A_n:=\frac1n\sum_{i=0}^{n-1}X_i\;\;\;\text{for }n\in\mathbb N$$ and $$\tau_n:=\frac12+\frac1{\operatorname{Var}[X_0]}\sum_{k=1}^{n-1}\left(1-\frac kn\right)\operatorname{Cov}[X_0,X_k]\;\;\;\text{for }n\in\mathbb N.$$ I know that $$\operatorname{Var}[A_n]=\frac{\operatorname{Var}[X_0]}n2\tau_n\;\;\;\text{for all }n\in\mathbb N\tag1$$ and $$n\operatorname{Var}[A_n]\xrightarrow{n\to\infty}\operatorname{Var}[X_0]+2\sum_{k=1}^\infty\operatorname{Cov}[X_0,X_k]\tag2.$$

Maybe the quantity corresponding to equation $(1)$ from the paper is a truncated version from the right-hand side of $(2)$? In any case, how can we compute the right quantity for the sample $x_0,\ldots,x_n$ (in a numerically stable way)?

EDIT: I need to implement the computation in C++, but a pseudocode for the computation would be sufficient for me.

Answer 1

There are two aspects to your question: (1) the link between the definition of $\tau_n$ and equation 2, and (2) how to compute the asymptotic variance. I won't give a formal proof, but will try to give the intuition.

First, note that if the chain is initialized from the stationary distribution (or, alternatively, if you have removed warm-up/burn-in), then the $X_0$ and $X_k$ are identically distributed and you can rewrite $\tau_n$ using correlations instead of covariances:

$$\tau_n=\frac12+\sum_{k=1}^{n-1}\left(1-\frac kn\right)\operatorname{Cor}[X_0,X_k]\;\;\;\text{for }n\in\mathbb N.$$

Let us assume that you have run your chain for very large $n$, so that your chain has mixed well. For $k$ large enough, $\operatorname{Cor}[X_0,X_k]\approx 0$ (the Markov chain "forgets its past"), so you can truncate the series

$$\tau_n\approx\frac12+\sum_{k=1}^{K}\left(1-\frac kn\right)\operatorname{Cor}[X_0,X_k]$$

for some $K<<n$. But for the remaining terms in the series, $\frac kn <<1$ so that

$$\tau_n\approx\frac12+\sum_{k=1}^{K}\operatorname{Cor}[X_0,X_k]$$

which corresponds to the truncated version of the series in equation 2.

Finally, note that once the chain has reached stationarity, $\operatorname{Cor}[X_0,X_k] = \operatorname{Cor}[X_t,X_{t+k}]$ for all $t$. You have therefore an unbiased estimator of $\operatorname{Cor}[X_0,X_k]$: the empirical correlation between the vectors $(X_0, X_1, \ldots, X_{n-k})$ and $(X_k, X_{k+1}, \ldots, X_n)$. For $k<<n$, the size of these vectors is large enough that the estimate will be numerically stable.

The numerical computation is now easy as long as you choose $K$ appropriately. I think it is standard to choose $K \approx \arg\min_k \left\{\hat{Cor}[X_0,X_k]\leq 0\right\}$.

For instance, here is a plot from an MCMC run, obtained thanks to the R function acf; it shows $\operatorname{Cor}[X_t,X_{t+k}]$ against $k$. In this plot, I would truncate the sum at $K=100$ (and would usually only display the plot up to that point. In the code below, the first line draws the plot and the second computes $\tau$. (The $-0.5$ term is because acf returns the series starting at $k=0$ instead of $k=1$, so you need to substract $\frac12$ instead of adding it.)

acf(X, lag.max=500) 
tau = sum(acf(X, lag.max=100, plot=F)$acf) - 0.5