beyondWhittle (version 1.1.1)

gibbs_np: Gibbs sampler for Bayesian nonparametric inference with Whittle likelihood

Description

Obtain samples of the posterior of the Whittle likelihood in conjunction with a Bernstein-Dirichlet prior on the spectral density.

Usage

gibbs_np(data, Ntotal, burnin, thin = 1, print_interval = 100,
  numerical_thresh = 1e-07, M = 1, g0.alpha = 1, g0.beta = 1,
  k.theta = 0.01, tau.alpha = 0.001, tau.beta = 0.001, kmax = 100 *
  coars + 500 * (!coars), trunc_l = 0.1, trunc_r = 0.9, coars = F,
  L = max(20, length(data)^(1/3)))

Arguments

data

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

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

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)

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)

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)

L

truncation parameter of DP in stick breaking representation

Value

list containing the following fields:

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 PSD parameters

lpost

trace of log posterior

Details

Further details can be found in the simulation study section in the references papers.

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>

N. Choudhuri et al. (2004) Bayesian Estimation of the Spectral Density of a Time Series JASA <doi:10.1198/016214504000000557>

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

Examples

Run this code
# NOT RUN {
##
## Example 1: Fit the NP model to sunspot data:
##

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_np(data=data, Ntotal=10000, burnin=4000, thin=2)

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


##
## Example 2: Fit the NP model to high-peaked AR(1) data
##

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_np(data=data, 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="")
# }

Run the code above in your browser using DataCamp Workspace