gibbs_bspline {bsplinePsd}

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

Description

This function updates the B-spline prior using the Whittle likelihood and obtains samples from the pseudo-posterior to infer the spectral density of a stationary time series.

Usage

gibbs_bspline(data, Ntotal, burnin, thin = 1, k.theta = 0.01, MG = 1,
  G0.alpha = 1, G0.beta = 1, LG = 20, MH = 1, H0.alpha = 1,
  H0.beta = 1, LH = 20, tau.alpha = 0.001, tau.beta = 0.001,
  kmax = 100, k1 = 20, degree = 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 (post-processing)

k.theta

prior parameter for number of B-spline densities k (proportional to exp(-k.theta*k^2)) in mixture

MG

Dirichlet process base measure constant for weights of B-spline densities in mixture (> 0)

G0.alpha, G0.beta

parameters of Beta base measure of Dirichlet process for weights of B-spline densities in mixture (default is Uniform[0, 1])

LG

truncation parameter of Dirichlet process in stick breaking representation for weights of B-spline densities

MH

Dirichlet process base measure constant for knot placements of B-spline densities (> 0)

H0.alpha, H0.beta

parameters of Beta base measure of Dirichlet process for knot placements of B-spline densities (default is Uniform[0, 1])

LH

truncation parameter of Dirichlet process in stick breaking representation for knot placements of B-spline densities

tau.alpha, tau.beta

prior parameters for tau (Inverse-Gamma)

kmax

upper bound for number of B-spline densities in mixture

k1

starting value for k. If k1 = NA then a random starting value between degree + 2 and kmax is selected

degree

positive integer specifying the degree of the B-spline densities (default is 3)

Details

The function gibbs_bspline is an implementation of the (serial version of the) MCMC algorithm presented in Edwards et al. (2018). This algorithm uses a nonparametric B-spline prior to estimate the spectral density of a stationary time series and can be considered a generalisation of the algorithm of Choudhuri et al. (2004), which used the Bernstein polynomial prior. A Dirichlet process prior is used to find the weights for the B-spline densities used in the finite mixture and a seperate and independent Dirichlet process prior used to place knots. The algorithm therefore allows for a data-driven choice of the number of knots/mixtures and their locations.

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

k, tau, V, Z, U, X

posterior traces of model parameters

knots.trace

trace of knot placements

ll.trace

trace of log likelihood

pdgrm

periodogram

n

integer length of input time series

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>.

Choudhuri, N., Ghosal, S., and Roy, A. (2004), Bayesian estimation of the spectral density of a time series, Journal of the American Statistical Association, 99(468):1050–1059.

See Also

plot.psd

Examples

## Not run: 

set.seed(123456)

# 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)
mcmc = gibbs_bspline(data, 10000, 5000)

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(mcmc)  # 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)