gibbs_vnp: Gibbs sampler for multivaiate Bayesian nonparametric...

In beyondWhittle: Bayesian Spectral Inference for Time Series

Gibbs sampler for multivaiate Bayesian nonparametric inference with Whittle likelihood

Description

Obtain samples of the posterior of the multivariate Whittle likelihood in conjunction with an Hpd AGamma process prior on the spectral density matrix.

Usage

gibbs_vnp(
  data,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 100,
  numerical_thresh = 1e-07,
  adaption.N = burnin,
  adaption.batchSize = 50,
  adaption.tar = 0.44,
  eta = ncol(data),
  omega = ncol(data),
  Sigma = 10000 * diag(ncol(data)),
  k.theta = 0.01,
  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 matrix; 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 eigenvalues of the spectral density
adaption.N	total number of iterations, in which the proposal variances (of r and U) are adapted
adaption.batchSize	batch size of proposal adaption
adaption.tar	target acceptance rate for adapted parameters
eta	AGamma process parameter, real number > ncol(data)-1
omega	AGamma process parameter, positive constant
Sigma	AGamma process parameter, Hpd matrix
k.theta	prior parameter for polynomial degree k (propto exp(-k.theta*k*log(k)))
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 Gamma process

Details

A detailed description of the method can be found in Section 5 in Meier (2018).

Value

list containing the following fields:

r, x, U	traces of the AGamma process parameters
k	posterior trace of polynomial degree
psd.median, psd.mean	psd estimates: (pointwise, componentwise) posterior median and mean
psd.p05, psd.p95	pointwise credibility interval
psd.u05, psd.u95	uniform credibility interval, see (6.5) in Meier (2018)
lpost	trace of log posterior

References

A. Meier (2018) A Matrix Gamma Process and Applications to Bayesian Analysis of Multivariate Time Series PhD thesis, OvGU Magdeburg <https://opendata.uni-halle.de//handle/1981185920/13470>

Examples

## Not run: 

##
## Example: Fit multivariate NP model to SOI/Recruitment series:
##

data <- cbind(as.numeric(astsa::soi-mean(astsa::soi)), 
              as.numeric(astsa::rec-mean(astsa::rec)) / 50)
data <- apply(data, 2, function(x) x-mean(x))

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_vnp(data=data, Ntotal=10000, burnin=4000, thin=2)

# Visualize results
plot(mcmc, log=T)


##
## Example 2: Fit multivariate NP model to VMA(1) data
##

n <- 256
ma <- rbind(c(-0.75, 0.5), c(0.5, 0.75))
Sigma <- rbind(c(1, 0.5), c(0.5, 1))
data <- sim_varma(model=list(ma=ma), n=n, d=2)
data <- apply(data, 2, function(x) x-mean(x))

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_vnp(data=data, Ntotal=10000, burnin=4000, thin=2)

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

## End(Not run)

beyondWhittle documentation built on June 22, 2024, 11:35 a.m.