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

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.

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<. 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<, 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 • Thank you for your answer. (a) Why$\operatorname{Cor}[X_0,X_k]\approx 0$for large$k$? (b) How do you numerically compute$\operatorname{Cor}[X_0,X_k]$? (c) In your plot, which target and proposal did you choose? Mar 2 '19 at 22:01 • I have updated the answer. (a) for large$k$, the chain has "forgotten its past", so the correlation goes to 0; (b) I have added the R code to the answer; (c) the plot corresponds to the estimate of a regression coefficient in a probit model, estimated using Metropolis-Hastings with a gaussian kernel. But the method is not specific to this example. Mar 3 '19 at 0:02 • Unfortunately, the R code doesn't help me much since I need to implement it with C++. Could provide some pseudocode how$\operatorname{Cor}[X_0,X_k]$is usually computed (in a stable way)? Mar 3 '19 at 10:54 • I think the argument you need is that$\forall t, Cor[X_0, X_k]=Cor[X_t, X_{t+k}]$(once the chain has reached stationarity). From this, you can derive an unbiased estimate of$Cor[X_0, X_k]$based on a sample of size$n-k\$, which will be numerically stable. Mar 3 '19 at 20:55
• @Taylor I'm still unsure how I need to compute the autocorrelation. Suppose x, ..., x[n - 1] are my Markov chain samples and let mean denote their mean. Should I then simply compute s = 0; for (i = 1; i < n; ++i) { s += (x - mean) * (x[i] - mean); } and then return 1 / (1 + 2 * s)? Is this what is done for Table 1 of the paper? Mar 4 '19 at 22:21