gibbs_ar: Gibbs sampler for an autoregressive model with PACF...

In beyondWhittle: Bayesian Spectral Inference for Time Series

Gibbs sampler for an autoregressive model with PACF parametrization.

Description

Obtain samples of the posterior of a Bayesian autoregressive model of fixed order.

Usage

gibbs_ar(
  data,
  ar.order,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 500,
  numerical_thresh = 1e-07,
  adaption.N = burnin,
  adaption.batchSize = 50,
  adaption.tar = 0.44,
  full_lik = F,
  rho.alpha = rep(1, ar.order),
  rho.beta = rep(1, ar.order),
  sigma2.alpha = 0.001,
  sigma2.beta = 0.001
)

Arguments

data	numeric vector; NA values are interpreted as missing values and treated as random
ar.order	order of the autoregressive model (integer >= 0)
Ntotal	total number of iterations to run the Markov chain
burnin	number of initial iterations to be discarded
thin	thinning number (postprocessing)
print_interval	Number of iterations, after which a status is printed to console
numerical_thresh	Lower (numerical pointwise) bound for the spectral density
adaption.N	total number of iterations, in which the proposal variances (of rho) are adapted
adaption.batchSize	batch size of proposal adaption for the rho_i's (PACF)
adaption.tar	target acceptance rate for the rho_i's (PACF)
full_lik	logical; if TRUE, the full likelihood for all observations is used; if FALSE, the partial likelihood for the last n-p observations
rho.alpha, rho.beta	prior parameters for the rho_i's: 2*(rho-0.5)~Beta(rho.alpha,rho.beta), default is Uniform(-1,1)
sigma2.alpha, sigma2.beta	prior parameters for sigma2 (inverse gamma)

Details

Partial Autocorrelation Structure (PACF, uniform prior) and the residual variance sigma2 (inverse gamma prior) is used as model parametrization. The DIC is computed with two times the posterior variance of the deviance as effective number of parameters, see (7.10) in the referenced book by Gelman et al. Further details can be found in the simulation study section in the referenced paper by C. Kirch et al. For more information on the PACF parametrization, see the referenced paper by Barndorff-Nielsen and Schou.

Value

list containing the following fields:

rho	matrix containing traces of the PACF parameters (if p>0)
sigma2	trace of sigma2
DIC	a list containing the numeric value DIC of the Deviance Information Criterion (DIC) and the effective number of parameters ENP
psd.median, psd.mean	psd estimates: (pointwise) posterior median and mean
psd.p05, psd.p95	pointwise credibility interval
psd.u05, psd.u95	uniform credibility interval
lpost	trace of log posterior

References

C. Kirch et al. (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

A. Gelman et al. (2013) Bayesian Data Analysis, Third Edition

O. Barndorff-Nielsen and G. Schou On the parametrization of autoregressive models by partial autocorrelations Journal of Multivariate Analysis (3),408-419 <doi:10.1016/0047-259X(73)90030-4>

Examples

## Not run: 

##
## Example 1: Fit an AR(p) model to sunspot data:
##

# Use this variable to set the AR model order
p <- 2

data <- sqrt(as.numeric(sunspot.year))
data <- data - mean(data)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_ar(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)


##
## Example 2: Fit an AR(p) model to high-peaked AR(1) data
##

# Use this variable to set the AR model order
p <- 1

n <- 256
data <- arima.sim(n=n, model=list(ar=0.95)) 
data <- data - mean(data)
omega <- fourier_freq(n)
psd_true <- psd_arma(omega, ar=0.95, ma=numeric(0), sigma2=1)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_ar(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Compare estimate with true function (green)
plot(mcmc, log=F, pdgrm=F, credib="uniform")
lines(x=omega, y=psd_true, col=3, lwd=2)

# Compute the Integrated Absolute Error (IAE) of posterior median
cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")

## End(Not run)

beyondWhittle documentation built on June 22, 2024, 11:35 a.m.