gibbs_pspline: Metropolis-within-Gibbs sampler for spectral inference of a...

In pmat747/psplinePsd: Spectral density estimation using P-spline priors

Metropolis-within-Gibbs sampler for spectral inference of a stationary time series using a P-spline prior

Description

This function uses the Whittle likelihood and obtains samples from the pseudo-posterior to infer the spectral density of a stationary time series. A P-spline prior is allocated on the spectral density function.

Usage

gibbs_pspline(
  data,
  Ntotal,
  burnin,
  thin = 1,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  phi.alpha = 1,
  phi.beta = 1,
  delta.alpha = 1e-04,
  delta.beta = 1e-04,
  k = NULL,
  eqSpacedKnots = FALSE,
  degree = 3,
  diffMatrixOrder = 2,
  printIter = 100,
  psd = NULL,
  add = FALSE
)

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 (post-processing)
tau.alpha, tau.beta	prior parameters for tau (Inverse-Gamma)
phi.alpha, phi.beta	prior parameters for phi (Gamma)
delta.alpha, delta.beta	prior parameters for delta (Gamma)
k	number of B-spline densities in the mixture
eqSpacedKnots	logical value indicating whether the knots are equally spaced or defined according to the periodogram
degree	positive integer specifying the degree of the B-spline densities (default is 3)
diffMatrixOrder	positive integer specifying the order of the difference penalty matrix in the P-splines (default is 2)
printIter	positive integer specifying the periodicity of the iteration number to be printed on screen (default 100)
psd	output from gibbs_pspline function
add	logical value indicating whether to add pilot posterior samples in the "psd" object to the current analysis

Details

The function gibbs_pspline is an implementation of the (serial version of the) MCMC algorithm presented in Maturana-Russel et al. (2019). This algorithm uses a P-spline prior to estimate the spectral density of a stationary time series and is similar to the B-spline prior algorithm of Edwards et al. (2018), which used a B-spline prior allowing the number of B-spline densities and knot locations to be variable. We define the prior on the spectral density as

f(w) = \tau \sum_{j=1}^{k}w_{j}B_{j}(w)

where B_{j} is the B-spline density. The following prior is allocated indirectly on the weights w_j:

v|\phi, \delta \sim N_{k-1}(0, (\phi D^\top D)^{-1})

\phi|\delta \sim Gamma(\alpha_{\phi}, \delta \beta_{\phi})

\delta \sim Gamma(\alpha_{\delta}, \beta_{\delta})

where

v_{j} = \log \left( \frac{w_{j}}{1-\sum_{j=1}^{k-1} w_{j}} \right)

Value

A list with S3 class ‘psd’ containing the following components:

psd.median,psd.mean	psd estimates: (pointwise) posterior median and mean
psd.p05,psd.p95	90% pointwise credibility interval
psd.u05,psd.u95	90% uniform credibility interval
fpsd.sample	posterior spectral density estimates
anSpecif	a list with some of the specifications of the analysis
n	integer length of input time series
tau,phi,delta,V	posterior traces of model parameters
ll.trace	trace of log likelihood
pdgrm	periodogram
db.list	B-spline densities
DIC	deviance information criterion
count	acceptance probabilities for the weigths

References

Edwards, M. C., Meyer, R., and Christensen, N. (2018), Bayesian nonparametric spectral density estimation using B-spline priors, Statistics and Computing, <https://doi.org/10.1007/s11222-017-9796-9>.

Maturana-Russel, P., and Meyer, R. (2019), Spectral density estimation using P-spline priors. arXiv:1905.01832.

See Also

plot.psd

Examples

## Not run: 

set.seed(1)

# Generate AR(1) data with rho = 0.9
n = 128
data = arima.sim(n, model = list(ar = 0.9))
data = data - mean(data)

# Run MCMC (may take some time)
pilotmcmc = gibbs_pspline(data, 2500, 500); # pilot run used in mcmc1 analysis
mcmc1 = gibbs_pspline(data, 3000, 2000, psd = pilotmcmc);
mcmc2 = gibbs_pspline(data, 3000, 0, psd = mcmc1, add = TRUE); # reciclying mcmc1 samples

require(beyondWhittle)  # For psd_arma() function
freq = 2 * pi / n * (1:(n / 2 + 1) - 1)[-c(1, n / 2 + 1)]  # Remove first and last frequency
psd.true = psd_arma(freq, ar = 0.9, ma = numeric(0), sigma2 = 1)  # True PSD

plot(mcmc1)  # Plot log PSD (see documentation of plot.psd)
lines(freq, log(psd.true), col = 2, lty = 3, lwd = 2)  # Overlay true PSD

plot(mcmc2)  # Plot log PSD (see documentation of plot.psd)
lines(freq, log(psd.true), col = 2, lty = 3, lwd = 2)  # Overlay true PSD


## End(Not run)

pmat747/psplinePsd documentation built on July 7, 2023, 9:06 p.m.