beyondWhittle (version 0.18.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, M = 1, g0.alpha = 1, g0.beta = 1,
  k.theta = 0.01, tau.alpha = 0.001, tau.beta = 0.001, kmax = 500,
  L = max(20, length(data)^(1/3)))

Arguments

data

numeric vector

Ntotal

total number of iterations to run the Markov chain

burnin

number of initial iterations to be discarded

thin

thinning number (postprocessing)

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

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 nonparametric correction

Details

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

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>

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

Examples

# 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=50000, burnin=20000, thin=4)

# 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 NP results on logarithmic scale"))
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 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)
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_NP(data=data, Ntotal=50000, burnin=20000, thin=4)

# 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 NP results"))
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="")
# }