gibbs_vnp {beyondWhittle}

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)