beyondWhittle (version 1.0)

gibbs_NPC_nuisance: Bayesian semiparametric inference in nuisance model

Description

Obtain samples of the posterior of the Bayesian model X_t=g(theta)+e_t, where the time series part e_t is modeled semiparametrically (by the corrected AR likelihood in conjunction with a Bernstein-Dirichlet prior on the spectral density correction) and the link function g (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_NPC_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), eta = T,
  M = 1, g0.alpha = 1, g0.beta = 1, k.theta = 0.01, tau.alpha = 0.001,
  tau.beta = 0.001, trunc_l = 0.1, trunc_r = 0.9, coars = F,
  kmax = 100 * coars + 500 * (!coars), L = max(20, length(data)^(1/3)))

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)

eta

logical variable indicating whether the model confidence eta should be included in the inference (eta=T) or fixed to 1 (eta=F)

M

DP base measure constant (> 0)

g0.alpha, g0.beta

parameters of Beta base measure of DP

k.theta

prior parameter for polynomial degree k (propto exp(-k.theta*k*log(k)))

tau.alpha, tau.beta

prior parameters for tau (inverse gamma)

trunc_l, trunc_r

left and right truncation of Bernstein polynomial basis functions, 0<=trunc_l<trunc_r<=1

coars

flag indicating whether coarsened or default bernstein polynomials are used (see Appendix E.1 in Ghosal and van der Vaart 2017)

kmax

upper bound for polynomial degree of Bernstein-Dirichlet mixture (can be set to Inf, algorithm is faster with kmax<Inf due to pre-computation of basis functions, but values 500<kmax<Inf are very memory intensive)

L

truncation parameter of DP in stick breaking representation

Value

list containing the following fields:

theta

posterior traces of the parameter of interest

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

k,tau,V,W

posterior traces of nonparametric correction

rho

posterior trace of model AR parameters (PACF parametrization)

eta

posterior trace of model confidence eta

Details

See gibbs_NPC for further details on the corrected parametric time series part. See nuisanceModel_mean or nuisanceModel_linearTrend for two examplary nuisance models.

References

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

S. Ghosal and A. van der Vaart (2017) Fundamentals of Nonparametric Bayesian Inference <DOI:10.1017/9781139029834>

Examples

# NOT RUN {
##
## Example 1: Fit Linear trend model with semiparametric 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 <- 50000; burnin <- 20000; thin <- 4; Neffective <- (Ntotal-burnin)/thin
mcmc <- gibbs_NPC_nuisance(data, nuisanceModel, Ntotal=Ntotal,
                          burnin=burnin, thin=thin, ar.order=ar.order, eta=F)

# 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/ NPC (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 NPC (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 semiparametric 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 <- 50000; burnin <- 20000; thin <- 4; Neffective <- (Ntotal-burnin)/thin
mcmc <- gibbs_NPC_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) NPC (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))
# }