View source: R/dynamicWhittle_prior_and_mcmc_params.R
gibbs_bdp_dw | R Documentation |
BDP-DW method: performing posterior sampling and calculating statistics based on the posterior samples
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
)
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 |
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 |
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 |
Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.