14

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.

Gabriel
  • 3,072
  • 1
  • 22
  • 49
Yang
  • 263
  • 2
  • 8

2 Answers2

4

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)
tchakravarty
  • 8,442
  • 2
  • 36
  • 50
0

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

xuhdev
  • 153
  • 6