beyondWhittle (version 1.0)

gibbs_AR: 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, 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)

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

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.

References

C. Kirch et al. (2017) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis <arXiv:1701.04846>

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

Examples

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

# Use this variable to set the autoregressive model order
ar.order <- 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, Ntotal=20000, burnin=8000, thin=4, ar.order=ar.order)

# Plot spectral estimates on log-scale (excluding the zero frequency).
N <- length(mcmc$psd.median)
pdgrm <- (abs(fft(data))^2 / (2*pi*length(data)))[1:N]
plot.ts(log(pdgrm[-1]), col="gray", 
  main=paste0("Sunspot AR results on logarithmic scale for p=", ar.order))
lines(log(mcmc$psd.median[-1]))
lines(log(mcmc$psd.p05[-1]),lty=2)
lines(log(mcmc$psd.p95[-1]),lty=2)
lines(log(mcmc$psd.u05[-1]),lty=3)
lines(log(mcmc$psd.u95[-1]),lty=3)
legend(x="topright", legend=c("periodogram", "pointwise median", 
  "pointwise CI", "uniform CI"), lty=c(1,1,2,3), col=c("gray", 1, 1, 1))


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

# Use this variable to set the autoregressive model order
ar.order <- 1

n <- 256
data <- arima.sim(n=n, model=list(ar=0.95)) 
data <- data - mean(data)
psd_true <- psd_arma(pi*omegaFreq(n), 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, Ntotal=20000, burnin=8000, thin=4, ar.order=1)

# Plot spectral estimates
N <- length(mcmc$psd.median)
pdgrm <- (abs(fft(data))^2 / (2*pi*length(data)))[1:N]
plot.ts(pdgrm[-1], col="gray",
  main=paste0("AR(1) data AR results for p=", ar.order))
lines(mcmc$psd.median[-1])
lines(mcmc$psd.p05[-1],lty=2)
lines(mcmc$psd.p95[-1],lty=2)
lines(mcmc$psd.u05[-1],lty=3)
lines(mcmc$psd.u95[-1],lty=3)
lines(psd_true[-1],col=2)
legend(x="topright", legend=c("periodogram", "true psd", 
  "pointwise median", "pointwise CI", "uniform CI"), lty=c(1,1,1,2,3), 
  col=c("gray", "red", 1, 1, 1))

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