MCMC Geweke diagnostic

I'm running a Metropolis sampler (C++) and want to use the previous samples to estimate the convergence rate.

One easy to implement diagnostic I found is the Geweke diagnostic, which computes the difference between the two sample means divided by itsestimated standard error. The standard error is estimated from the spectral density at zero.

$$Z_n=\frac{\bar{\theta}_A-\bar{\theta}_B}{\sqrt{\frac{1}{n_A}\hat{S_{\theta}^A}(0)+\frac{1}{n_B}\hat{S_{\theta}^B}(0)}},$$

where $A$, $B$ are two windows within the Markov chain. I did some research on what are $\hat{S_{\theta}^A}(0)$ and $\hat{S_{\theta}^B}(0)$ but get into a mess of literature on energy spectral density and power spectral density. (http://en.wikipedia.org/wiki/Spectral_density) but I'm not an expert on these topics; I just need a quick answer: are these quantities the same as sample variance? If not, what's the formula for computing them?

Another doubt on this Geweke diagnostic is how to pick $\theta$? The above literature said that it is some functional $\theta(X)$ and should imply an existence of a spectral density $\hat{S_{\theta}^A}(0)$, but for convenience I guess the simplest way is to use the identity function (use samples themselves). Is this correct?

The R coda package has a description but neither does it specify how to compute the $S$ values.

-
you could look in the guts of the coda function geweke.diag to see what it does ... – Ben Bolker Aug 21 '12 at 16:54

You can look through the code for the geweke.diag function in the coda package to see how the variance is computed, via the call to the spectrum.ar0 function.

Here is a short motivation of the computation of the spectral density of an AR($p$) process at zero.

The spectral density of an AR($p$) process at frequency $\lambda$ is given by the expression: $$f(\lambda) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j\exp(-2\pi\iota j\lambda))^2}$$ where $\alpha_j$ are the autoregressive parameters.

This expression simplifies considerably when computing the spectral density of an AR($p$) process at $0$: $$f(0) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j)^2}$$

The computation then would look something like this (substituting the usual estimators for parameters):

tsAR2 = arima.sim(list(ar = c(0.01, 0.03)), n = 1000)  # simulate an AR(2) process
ar2 = ar(tsAR2, aic = TRUE)  # estimate it with AIC complexity selection

# manual estimate of spectral density at zero
sdMan = ar2$var.pred/(1-sum(ar2$ar))^2

# coda computation of spectral density at zer0
sdCoda = coda::spectrum0.ar(tsAr2)$spec # assert equality all.equal(sdCoda, sdMan)  - Check the wikipedia page. You'll see$S_{xx}(\omega)$, which is the spectral density. In your case, you should use$S_{xx}(0)\$.

-