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 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.
 A: 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)

A: Check the wikipedia page. You'll see $S_{xx}(\omega)$, which is the spectral density. In your case, you should use $S_{xx}(0)$.
