gibbs_bdp_dw {beyondWhittle}

BDP-DW method: performing posterior sampling and calculating statistics based on the posterior samples

Description

BDP-DW method: performing posterior sampling and calculating statistics based on the posterior samples

Usage

gibbs_bdp_dw(
  data,
  m,
  likelihood_thinning = 1,
  monitor = TRUE,
  print_interval = 100,
  unif_CR = FALSE,
  res_time,
  freq,
  Ntotal = 110000,
  burnin = 60000,
  thin = 10,
  adaptive.batchSize = 50,
  adaptive.targetAcceptanceRate = 0.44,
  M = 1,
  g0.alpha = 1,
  g0.beta = 1,
  k1.theta = 0.01,
  k2.theta = 0.01,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  k1max = 100,
  k2max = 100,
  L = 20,
  bernstein1_l = 0.1,
  bernstein1_r = 0.9,
  bernstein2_l = 0.1,
  bernstein2_r = 0.9
)

Arguments

data

time series that needs to be analyzed

m

window size needed to calculate moving periodogram.

likelihood_thinning

the thinning factor of the dynamic Whittle likelihood.

monitor

a Boolean value (default TRUE) indicating whether to display the real-time status

print_interval

If monitor = TRUE, then this value indicates number of iterations after which a status is printed to console; If monitor = FALSE, it does not have any effect

unif_CR

a Boolean value (default FALSE) indicating whether to calculate the uniform credible region

res_time, freq

a set of grid lines in [0,1] and [0,\pi], respectively, specifying where to evaluate the estimated tv-PSD

Ntotal

total number of iterations to run the Markov chain

burnin

number of initial iterations to be discarded

thin

thinning number (for post-processing of the posterior sample)

adaptive.batchSize

the batch size for the adaptive MCMC algorithm for sampling tau

adaptive.targetAcceptanceRate

the target acceptance rate for the adaptive MCMC algorithm for sampling tau

M

DP base measure constant (> 0)

g0.alpha, g0.beta

parameters of Beta base measure of DP

k1.theta

prior parameter for polynomial corresponding to rescaled time (propto exp(-k1.theta*k1*log(k1)))

k2.theta

prior parameter for polynomial corresponding to rescaled frequency (propto exp(-k2.theta*k2*log(k2)))

tau.alpha, tau.beta

prior parameters for tau (inverse gamma)

k1max

upper bound of the degrees of Bernstein polynomial corresponding to rescaled time (for pre-computation of basis functions)

k2max

upper bound of the degrees of Bernstein polynomial corresponding to rescaled frequency (for pre-computation of basis functions)

L

truncation parameter of DP in stick breaking representation

bernstein1_l, bernstein1_r

left and right truncation of Bernstein polynomial basis functions for rescaled time, 0<=bernstein1_l<bernstein1_r<=1

bernstein2_l, bernstein2_r

left and right truncation of Bernstein polynomial basis functions for rescaled frequency, 0<=bernstein2_l<bernstein2_r<=1

Value

list containing the following fields:

k1, k2, tau, V, W1, W2

posterior traces of PSD parameters

lpost

traces log posterior

tim

total run time

bf_k1

Savage-Dickey estimate of Bayes factor of hypothesis k1=1

tvpsd.mean, tvpsd.median

posterior mean and pointwise posterior median (matrices of dimension length(rescaled_time) by length(freq))

tvpsd.p05, tvpsd.p95

90 percent pointwise credibility interval

tvpsd.u05, tvpsd.u95

90 percent uniform credibility interval if unif_CR = TRUE. Otherwise NA

References

Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Examples

## Not run: 

##
## Example: Applying BDP-DW method to a multi-peaked tvMA(1) process
##

# set seed
set.seed(200)
# set the length of time series
len_d <- 1500
# generate data from DGP LS2c defined in Section 4.2.2 of Tang et al. (2023). 
# see also ?sim_tvarma12
sim_data <- sim_tvarma12(len_d = 1500, dgp = "LS2", innov_distribution = "c")
# specify grid-points at which the tv-PSD is evaluated
res_time <- seq(0, 1, by = 0.005); freq <- pi * seq(0, 1, by = 0.01)
# calculate the true tv-PSD of DGP LS2c at the pre-specified grid
true_tvPSD <- psd_tvarma12(rescaled_time = res_time, freq = freq, dgp = "LS2")
# plot the true tv-PSD
# type ?plot.bdp_dw_tv_psd for more info
plot(true_tvPSD)

# If you run the example be aware that this may take several minutes
print("This example may take some time to run")
result <- gibbs_bdp_dw(data = sim_data, 
m = 50, 
likelihood_thinning = 2, 
rescaled_time = res_time,
 freq = freq)

# extract bayes factor and examine posterior summary
bayes_factor(result)
summary(result)

# compare estimate with true function
# type ?plot.bdp_dw_result for more info
par(mfrow = c(1,2))

plot(result, which = 1,
zlim = range(result$tvpsd.mean, true_tvPSD$tv_psd)
)
plot(true_tvPSD,
zlim = range(result$tvpsd.mean, true_tvPSD$tv_psd),
main = "true tv-PSD")

par(mfrow = c(1,1))


## End(Not run)