beyondWhittle (version 1.0)

gibbs_AR_nuisance: Bayesian parametric (AR) inference in nuisance model, with PACF parametrization

Description

Obtain samples of the posterior of the Bayesian model X_t=g_t(theta)+e_t, where the time series part e_t is AR(p) and the link function g_t (depending on the parameter of interest theta) is provided by the user. Examples for this framework include the simple mean model X_t=mu+e_t (theta=mu, g_t(theta)=theta) or the linear trend model X_t=bt+mu+e_t (theta=c(mu,b), g(theta)=bt+mu) for t=1,...,n.

Usage

gibbs_AR_nuisance(data, nuisanceModel, 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

nuisanceModel

model of class nuisanceModel, see nuisanceModel_mean and nuisanceModel_linearTrend for examples

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:

theta

posterior traces of the parameter of interest

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

See gibbs_AR for further details on the autoregressive time series part. See nuisanceModel_mean or nuisanceModel_linearTrend for two examplary nuisance models.

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 Linear trend model with parametric nuisance time series to temperature data
##

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

# Initialize the linear trend model
nuisanceModel <- nuisanceModel_linearTrend()

data <- as.numeric(nhtemp)
n <- length(data)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
Ntotal <- 20000; burnin <- 8000; thin <- 4; Neffective <- (Ntotal-burnin)/thin
mcmc <- gibbs_AR_nuisance(data, nuisanceModel, Ntotal=Ntotal,
                          burnin=burnin, thin=thin, ar.order=ar.order)

# Reconstruct linear trend lines from posterior traces of theta=c(mu,b)
trend_lines <- array(NA, dim=c(n,Neffective))
for (j in 1:Neffective) {
  mu_j <- mcmc$theta[1,j]
  b_j <- mcmc$theta[2,j]
  trend_lines[,j] <- mu_j + b_j * (1:n)
}
trend.median <- apply(trend_lines, 1, median)
trend.p05 <- apply(trend_lines, 1, quantile, 0.05) # 90% CI
trend.p95 <- apply(trend_lines, 1, quantile, 0.95) # "
trend.p025 <- apply(trend_lines, 1, quantile, 0.025) # 95% CI
trend.p975 <- apply(trend_lines, 1, quantile, 0.975) #

# Plot confidence bands for linear trend curve
par(mfcol=c(2,1),mar=c(4,4,2,2))
data_timePeriod <- start(nhtemp)[1]:end(nhtemp)[1]
plot(x=data_timePeriod, y=data,
     main=paste0("New Haven temperature w/ AR (p=", ar.order,
                 ") linear trend estimate"),
     col="gray", type="l", xlab="Year", ylab="Avg temp. (deg. F)")
lines(x=data_timePeriod, y=trend.median)
lines(x=data_timePeriod, y=trend.p05, lty=2)
lines(x=data_timePeriod, y=trend.p95, lty=2)
lines(x=data_timePeriod, y=trend.p025, lty=3)
lines(x=data_timePeriod, y=trend.p975, lty=3)
legend(x="topleft", legend=c("data", "posterior median",
                         "posterior 90% CI", "posterior 95% CI"),
       lty=c(1,1,2,3), col=c("gray", 1, 1, 1),
       cex=.75, x.intersp=.5, y.intersp=.5)

# Plot spectral estimate
N <- length(mcmc$psd.median)
pdgrm <- (abs(fft(data))^2 / (2*pi*length(data)))[1:N]
plot.ts((pdgrm[-1]), col="gray",
        main=paste0("New Haven temperature AR (p=",
                    ar.order, ") spectral estimate"))
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)
legend(x="top", legend=c("periodogram", "pointwise median",
                              "pointwise CI", "uniform CI"),
       lty=c(1,1,2,3), col=c("gray", 1, 1, 1),
       cex=.75, x.intersp=.5, y.intersp=.5)
par(mfcol=c(1,1))


##
## Example 2: Fit mean model (with autoregressive nuisance time series) to AR(1) data
##

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

n <- 256
mu_true <- 3
data <- arima.sim(n=n, model=list(ar=0.75))
data <- data + mu_true
psd_true <- psd_arma(pi*omegaFreq(n), ar=0.95, ma=numeric(0), sigma2=1)

# Initialize the mean model
nuisanceModel <- nuisanceModel_mean()

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
Ntotal <- 20000; burnin <- 8000; thin <- 4; Neffective <- (Ntotal-burnin)/thin
mcmc <- gibbs_AR_nuisance(data, nuisanceModel, Ntotal=Ntotal,
                          burnin=burnin, thin=thin, ar.order=ar.order)

# Plot posterior trace of mu
par(mfcol=c(2,1),mar=c(4,2,2,2))
plot.ts(mcmc$theta[1,], main="Posterior trace of mu", xlab="Iteration")
abline(h=mu_true, col=2, lty=2)
abline(h=mean(data), col=3, lty=2)
legend(x="topright", legend=c("true mean", "empirical mean of data"), lty=c(2,2), col=c(2,3))

# Plot spectral estimate
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) AR (p=", ar.order, ") spectral estimate"))
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)
legend(x="top", legend=c("periodogram", "pointwise median",
                         "pointwise CI", "uniform CI"),
       lty=c(1,1,2,3), col=c("gray", 1, 1, 1))
par(mfcol=c(1,1))
# }