gibbs_npc {beyondWhittle}R Documentation

Gibbs sampler for Bayesian semiparametric inference with the corrected AR likelihood


Obtain samples of the posterior of the corrected autoregressive likelihood in conjunction with a Bernstein-Dirichlet prior on the correction.


gibbs_npc(data, ar.order, Ntotal, burnin, thin = 1,
  print_interval = 100, 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)))



numeric vector; NA values are interpreted as missing values and treated as random


order of the autoregressive model (integer > 0)


total number of iterations to run the Markov chain


number of initial iterations to be discarded


thinning number (postprocessing)


Number of iterations, after which a status is printed to console


Lower (numerical pointwise) bound for the spectral density


total number of iterations, in which the proposal variances (of rho) are adapted


batch size of proposal adaption for the rho_i's (PACF)


target acceptance rate for the rho_i's (PACF)


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)


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


DP base measure constant (> 0)

g0.alpha, g0.beta

parameters of Beta base measure of DP


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


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


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)


truncation parameter of DP in stick breaking representation


Partial Autocorrelation Structure (PACF, uniform prior) and the residual variance sigma2 (inverse gamma prior) is used as model parametrization. A Bernstein-Dirichlet prior for c_eta with base measure Beta(g0.alpha, g0.beta) is used. Further details can be found in the simulation study section in the referenced paper by Kirch et al. For more information on the PACF parametrization, see the referenced paper by Barndorff-Nielsen and Schou.


list containing the following fields:


psd estimates: (pointwise) posterior median and mean


pointwise credibility interval


uniform credibility interval


posterior traces of nonparametric correction


posterior trace of model AR parameters (PACF parametrization)


posterior trace of model confidence eta


trace of log posterior


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>

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

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>


## Not run: 

## Example 1: Fit a nonparametrically corrected 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_npc(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 a nonparametrically corrected AR(p) model to high-peaked AR(1) data

# Use this variable to set the autoregressive 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_npc(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)

[Package beyondWhittle version 1.1.1 Index]