Diagnóstico MCMC Geweke

Estou executando um amostrador Metropolis (C ++) e quero usar os exemplos anteriores para estimar a taxa de convergência.

Um diagnóstico fácil de implementar que encontrei é o diagnóstico Geweke , que calcula a diferença entre as duas médias da amostra divididas pelo erro padrão estimado. O erro padrão é estimado a partir da densidade espectral em 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)}},$$

onde , B são duas janelas dentro da cadeia de Markov. Eu fiz algumas pesquisas sobre o que são ^ S A θ ( 0 ) e ^ S B θ ( 0 ), mas entrei em uma confusão de literatura sobre densidade espectral de energia e densidade espectral de potência, mas não sou especialista nesses tópicos; Eu só preciso de uma resposta rápida: essas quantidades são iguais à variação da amostra? Se não, qual é a fórmula para calculá-los?$A$$B$$\hat{S_{\theta}^A}(0)$$\hat{S_{\theta}^B}(0)$

Outra dúvida sobre esse diagnóstico de Geweke é como escolher ? A literatura acima disse que é algum θ ( X ) funcional e deveria implicar a existência de uma densidade espectral ^ S A θ ( 0 ) , mas por conveniência, acho que a maneira mais simples é usar a função de identidade (use as próprias amostras). Isso está correto?$\theta$$\theta(X)$$\hat{S_{\theta}^A}(0)$

O pacote R coda tem uma descrição, mas também não especifica como calcular os valores $S$

Você pode procurar no código a geweke.diagfunção no codapacote para ver como a variação é calculada, através da chamada para a spectrum.ar0função.

---

Aqui está uma breve motivação do cálculo da densidade espectral de um processo AR ( ) em zero.$p$

A densidade espectral de um processo AR ( ) na frequência $p$$\lambda$ é dada pela expressão: ondeαjsão os parâmetros autorregressivos.

$$f(\lambda) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j\exp(-2\pi\iota j\lambda))^2}$$ $\alpha_j$

Essa expressão simplifica consideravelmente ao calcular a densidade espectral de um AR ($p$$0$

$$f(0) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j)^2}$$

A computação ficaria assim (substituindo os estimadores usuais por parâmetros):

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)

$S_{xx}(\omega)$$S_{xx}(0)$