Type:	Package
Title:	Dynamic Shrinkage Process and Change Point Detection
Version:	1.2.0
Description:	Provides efficient Markov chain Monte Carlo (MCMC) algorithms for dynamic shrinkage processes, which extend global-local shrinkage priors to the time series setting by allowing shrinkage to depend on its own past. These priors yield locally adaptive estimates, useful for time series and regression functions with irregular features. The package includes full MCMC implementations for trend filtering using dynamic shrinkage on signal differences, producing locally constant or linear fits with adaptive credible bands. Also included are models with static shrinkage and normal-inverse-Gamma priors for comparison. Additional tools cover dynamic regression with time-varying coefficients and B-spline models with shrinkage on basis differences, allowing for flexible curve-fitting with unequally spaced data. Some support for heteroscedastic errors, outlier detection, and change point estimation. Methods in this package are described in Kowal et al. (2019) <doi:10.1111/rssb.12325>, Wu et al. (2024) <doi:10.1080/07350015.2024.2362269>, Schafer and Matteson (2024) <doi:10.1080/00401706.2024.2407316>, and Cho and Matteson (2024) <doi:10.48550/arXiv.2408.11315>.
License:	GPL (≥ 3)
Depends:	R (≥ 4.1.0)
Imports:	coda, fda, graphics, grDevices, Matrix, MCMCpack, methods, msm, pgdraw, Rcpp, RcppZiggurat, spam, progress, stats, stochvol, BayesLogit, truncdist, mgcv, purrr, rlang, lifecycle, glue
Suggests:	ggplot2, testthat (≥ 3.0.0)
LinkingTo:	Rcpp, RcppArmadillo, RcppEigen
Encoding:	UTF-8
RoxygenNote:	7.3.2
Config/testthat/edition:	3
URL:	https://github.com/schafert/dsp
BugReports:	https://github.com/schafert/dsp/issues
NeedsCompilation:	yes
Packaged:	2025-07-31 16:57:44 UTC; torynschafer
Author:	Daniel R. Kowal [aut, cph], Haoxuan Wu [aut], Toryn Schafer [aut, cre], Jason B. Cho [aut], David S. Matteson [aut]
Maintainer:	Toryn Schafer <toryn27@gmail.com>
Repository:	CRAN
Date/Publication:	2025-08-19 15:00:13 UTC

dsp: Dynamic Shrinkage Process and Change Point Detection

Description

Provides efficient Markov chain Monte Carlo (MCMC) algorithms for dynamic shrinkage processes, which extend global-local shrinkage priors to the time series setting by allowing shrinkage to depend on its own past. These priors yield locally adaptive estimates, useful for time series and regression functions with irregular features. The package includes full MCMC implementations for trend filtering using dynamic shrinkage on signal differences, producing locally constant or linear fits with adaptive credible bands. Also included are models with static shrinkage and normal-inverse-Gamma priors for comparison. Additional tools cover dynamic regression with time-varying coefficients and B-spline models with shrinkage on basis differences, allowing for flexible curve-fitting with unequally spaced data. Some support for heteroscedastic errors, outlier detection, and change point estimation. Methods in this package are described in Kowal et al. (2019) doi:10.1111/rssb.12325, Wu et al. (2024) doi:10.1080/07350015.2024.2362269, Schafer and Matteson (2024) doi:10.1080/00401706.2024.2407316, and Cho and Matteson (2024) doi:10.48550/arXiv.2408.11315.

Author(s)

Maintainer: Toryn Schafer toryn27@gmail.com (ORCID)

Authors:

• Daniel R. Kowal daniel.kowal@rice.edu [copyright holder]

• Haoxuan Wu hw399@cornell.edu

• Jason B. Cho bc454@cornell.edu

• David S. Matteson matteson@cornell.edu

See Also

Useful links:

• https://github.com/schafert/dsp

• Report bugs at https://github.com/schafert/dsp/issues

Adaptive Bayesian Changepoint with Outliers

Description

Run the MCMC sampler for ABCO with a penalty on first (D = 1), or second (D = 2) differences of the conditional expectation. The penalty utilizes the dynamic horseshoe prior on the evolution errors. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

abco(
  y,
  D = 1,
  useAnom = TRUE,
  obsSV = "const",
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "omega", "yhat", "evol_sigma_t2", "r", "zeta", "obs_sigma_t2",
    "zeta_sigma_t2", "dhs_phi", "dhs_mean", "h", "h_smooth"),
  verbose = TRUE,
  D_asv = 1,
  evol_error_asv = "HS",
  nugget_asv = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

MCMC Sampler for Bayesian Trend Filtering

Description

Run the MCMC for Bayesian trend filtering with a penalty on zeroth (D = 0), first (D = 1), or second (D = 2) differences of the conditional expectation. The penalty is determined by the prior on the evolution errors, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

btf(
  y,
  evol_error = "DHS",
  D = 2,
  obsSV = "const",
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "evol_sigma_t2", "obs_sigma_t2", "dhs_phi",
    "dhs_mean", "h", "h_smooth"),
  computeDIC = TRUE,
  verbose = TRUE,
  D_asv = 1,
  evol_error_asv = "HS",
  nugget_asv = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

MCMC Sampler for Bayesian Trend Filtering: D = 0

Description

Run the MCMC for Bayesian trend filtering with a penalty on the conditional expectation. The penalty is determined by the prior on the evolution errors, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

btf0(
  y,
  evol_error = "DHS",
  obsSV = "const",
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "evol_sigma_t2", "obs_sigma_t2", "dhs_phi",
    "dhs_mean", "h", "h_smooth"),
  computeDIC = TRUE,
  verbose = TRUE,
  D_asv = 1,
  evol_error_asv = "HS",
  nugget_asv = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

MCMC Sampler for B-spline Bayesian Trend Filtering

Description

Run the MCMC for B-spline fitting with a Bayesian trend filtering model on the coefficients, i.e., a penalty on zeroth (D=0), first (D=1), or second (D=2) differences of the B-spline basis coefficients. The penalty is determined by the prior on the evolution errors, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

btf_bspline(
  y,
  times = NULL,
  num_knots = NULL,
  evol_error = "DHS",
  D = 2,
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "beta", "evol_sigma_t2", "obs_sigma_t2", "dhs_phi",
    "dhs_mean"),
  computeDIC = TRUE,
  verbose = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

The primary advantages of btf_bspline over btf are

1. Unequally-spaced points are handled automatically and

2. Computations are linear in the number of basis coefficients, which may be substantially fewer than the number of time points.

MCMC Sampler for B-spline Bayesian Trend Filtering: D = 0

Description

Run the MCMC for B-spline fitting with a penalty the B-spline basis coefficients. The penalty is determined by the prior on the evolution errors, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

btf_bspline0(
  y,
  times = NULL,
  num_knots = NULL,
  evol_error = "DHS",
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "beta", "evol_sigma_t2", "obs_sigma_t2", "dhs_phi",
    "dhs_mean"),
  computeDIC = TRUE,
  verbose = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

The primary advantages of btf_bspline over btf are

1. Unequally-spaced points are handled automatically and

2. Computations are linear in the number of basis coefficients, which may be substantially fewer than the number of time points.

MCMC Sampler for Bayesian Trend Filtering with Negative Binomial Observation Model

Description

Run the MCMC for Bayesian trend filtering with a penalty on first (D = 1) or second (D = 2) differences of the log conditional expectation of a Negative Binomial distribution. The penalty is determined by the prior on the evolution errors, currently only the following options are available:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

btf_nb(
  y,
  evol_error = "DHS",
  D = 2,
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "evol_sigma_t2", "r", "dhs_phi", "dhs_mean"),
  r_init = NULL,
  r_sample = TRUE,
  step = 1,
  evol0_sample = FALSE,
  evol0_sd = 10,
  sigma_e = 1/sqrt(Nt),
  chol0 = NULL,
  computeDIC = TRUE,
  offset = 0,
  verbose = TRUE,
  seed = NULL
)

Arguments

Value

A named list with the following. One named element for each of the parameters specified in mcmc_params containing the nsave posterior draws as a vector for scalar parameters and matrix for multivariate parameters. An element named loglike containing a vector of the nsave evaluations of the log likelihood of the data y with a Negative Binomial distribution defined by the current iteration conditional expectation + offset and overdispersion. If computeDIC is TRUE, a named element for the computed value of DIC and effective parameters p_d for the data y with a Negative Binomial distribution defined by the posterior mean of the conditional expectation + offset and posterior mean overdispersion.

MCMC Sampler for Bayesian Trend Filtering: Regression

Description

Run the MCMC for Bayesian trend filtering regression with a penalty on first (D=1) or second (D=2) differences of each dynamic regression coefficient. The penalty is determined by the prior on the evolution errors, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

btf_reg(
  y,
  X = NULL,
  evol_error = "DHS",
  D = 1,
  obsSV = "const",
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "beta", "evol_sigma_t2", "obs_sigma_t2", "dhs_phi",
    "dhs_mean", "h", "h_smooth"),
  use_backfitting = FALSE,
  computeDIC = TRUE,
  verbose = TRUE,
  D_asv = 1,
  evol_error_asv = "HS",
  nugget_asv = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

Run the MCMC for sparse Bayesian trend filtering

Description

Sparse Bayesian trend filtering has two penalties: (1) a penalty on the first (D = 1) or second (D = 2) differences of the conditional expectation and (2) a penalty on the conditional expectation, i.e., shrinkage to zero.

Usage

btf_sparse(
  y,
  evol_error = "DHS",
  zero_error = "DHS",
  D = 2,
  obsSV = "const",
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "evol_sigma_t2", "obs_sigma_t2", "zero_sigma_t2",
    "dhs_phi", "dhs_mean", "dhs_phi_zero", "dhs_mean_zero", "h", "h_smooth"),
  computeDIC = TRUE,
  verbose = TRUE,
  D_asv = 1,
  evol_error_asv = "HS",
  nugget_asv = TRUE
)

Arguments

Details

Each penalty is determined by a prior, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the prior is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

Compute the quadratic term in Bayesian trend filtering

Description

Compute the quadratic term arising in the full conditional distribution of a Bayesian trend filtering model with D = 1 or D = 2. This function exploits the known D-banded structure of Q to compute the matrix directly, using objects in the Matrix package.

Usage

build_Q(obs_sigma_t2, evol_sigma_t2, D = 1)

Arguments

Value

Banded T x T Matrix (object) Q

Compute X'X

Description

Build the Tp x Tp matrix XtX using the Matrix() package

Usage

build_XtX(X)

Arguments

Value

Block diagonal Tp x Tp Matrix (object) where each p x p block is tcrossprod(matrix(X[t,]))

Note

X'X is a one-time computing cost. Special cases may have more efficient computing options, but the Matrix representation is important for efficient computations within the sampler.

Function for calculating DIC and Pb (Bayesian measures of model complexity and fit by Spiegelhalter et al. 2002)

Description

Function for calculating DIC and Pb (Bayesian measures of model complexity and fit by Spiegelhalter et al. 2002)

Usage

computeDIC_ASV(y, beta, post_sigma2, post_loglike)

Arguments

Value

a list containing DIC and p_d. Two options for estimating both DIC and p_d, which are both included.

Compute Simultaneous Credible Bands

Description

Compute (1-alpha)\

Usage

credBands(sampFuns, alpha = 0.05)

Arguments

Value

m x 2 matrix of credible bands; the first column is the lower band, the second is the upper band

Note

The input needs not be curves: the simultaneous credible "bands" may be computed for vectors. The resulting credible intervals will provide joint coverage at the (1-alpha)% level across all components of the vector.

MCMC Sampler for Models with Dynamic Shrinkage Processes

Description

Wrapper function for fitting models with Dynamic Shrinkage Processes (DSP), including:

• Adaptive Bayesian Changepoint analysis and local Outlier (ABCO),

• Bayesian Trend Filter for Gaussian Data

• Time-varying Regression

• Bayesian Trend Filter with B-spline for irregularly spaced or functional time-series.

• Bayesian Smoothing for Count Data

Method for printing basic information about the MCMC sampling settings for the fitted model

Usage

dsp_fit(
  y,
  model_spec,
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  computeDIC = TRUE,
  verbose = TRUE,
  ...
)

## S3 method for class 'dsp'
print(x, ...)

Arguments

Details

A brief summary of the settings used to fit the model including number of iterations, burn in, and thinning rates.

Value

dsp_fit returns an object of class "dsp".

An object of class "dsp" is defined as a list containing at least the following components:

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues when family is "gaussian".

Examples

set.seed(200)
signal = c(rep(0, 50), rep(10, 50))
noise = rep(1, 100)
noise_var = rep(1, 100)
for (k in 2:100){
  noise_var[k] = exp(0.9*log(noise_var[k-1]) + rnorm(1, 0, 0.5))
  noise[k] = rnorm(1, 0, sqrt(noise_var[k])) }

y = signal + noise
model_spec = dsp_spec(family = "gaussian", model = "changepoint",
                      D = 1, useAnom = TRUE, obsSV = "SV")
mcmc_output = dsp_fit(y, model_spec = model_spec, nsave = 500, nburn = 500)

print(mcmc_output)

Model Specification

Description

Method for creating dsp specification object prior to fitting.

Method for printing basic information about the model specification

Usage

dsp_spec(family, model, ...)

## S3 method for class 'dsp_spec'
print(x, ...)

Arguments

Value

A list containing the model specification.

Examples

model_spec <- dsp_spec(family = "gaussian",
                      model = "changepoint")


print(model_spec)

Compute the ergodic (running) mean.

Description

Compute the ergodic (running) mean.

Usage

ergMean(x)

Arguments

Value

A vector y with each element defined by y[i] = mean(x[1:i])

MCMC Sampler for Adaptive Stchoastic Volatility (ASV) model

Description

The penalty is determined by the prior on the evolution errors, which include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal-inverse-gamma prior ('NIG').

In each case, the evolution error is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

fit_ASV(
  y,
  beta = 0,
  evol_error = "DHS",
  D = 1,
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("h", "logy2hat", "sigma2", "evol_sigma_t2", "dhs_phi", "dhs_mean"),
  nugget = FALSE,
  computeDIC = TRUE,
  verbose = TRUE
)

Arguments

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

Helper function for Sampling parameters for ASV model

Description

Helper function for Sampling parameters for ASV model

Usage

fit_paramsASV(data, sParams, evol_error, D)

Arguments

Value

a list containing 4 sets of parameters

• s_p_error_term: matrix containing mean and the variance from 10-componenet gaussian mixture (Omori et al. 2007)

• s_mu: a vector containing the posterior sample of log variance h,

• s_evolParams0: a list containing posterior samples of parameters associated with the variance of first D observation of the log variance term, h.

• s_evolParams: a list containing posterior samples parameters associated with the variance of D to the last observations of the log variance temr , h.

Helper function for Sampling parameters for ASV model with a nugget Effect

Description

Helper function for Sampling parameters for ASV model with a nugget Effect

Usage

fit_paramsASV_n(data, sParams, evol_error, D)

Arguments

Value

a list containing 4 sets of parameters

• s_p_error_term: matrix containing mean and the variance from 10-componenet gaussian mixture (Omori et al. 2007)

• s_mu: a vector containing the posterior sample of log variance h,

• s_evolParams0: a list containing posterior samples of parameters associated with the variance of first D observation of the log variance term, h.

• s_evolParams: a list containing posterior samples parameters associated with the variance of D to the last observations of the log variance temr , h.

Posterior predictive sampler on the transformed y (log(y^2))

Description

Posterior predictive sampler on the transformed y (log(y^2))

Usage

generate_ly2hat(h, p_error_term)

Arguments

Value

a vector containing posterior predictive on log(y^2)

Compute the design matrix X for AR(p) model

Description

Compute the design matrix X for AR(p) model

Usage

getARpXmat(y, p = 1, include_intercept = FALSE)

Arguments

Summarize of effective sample size

Description

Compute the summary statistics for the effective sample size (ESS) across posterior samples for possibly many variables

Usage

getEffSize(postX)

Arguments

Value

Table of summary statistics using the function summary().

Compute Non-Zeros (Signals)

Description

Estimate the location of non-zeros (signals) implied by horseshoe-type thresholding.

Usage

getNonZeros(post_evol_sigma_t2, post_obs_sigma_t2 = NULL)

Arguments

Details

Thresholding is based on kappa[t] > 1/2, where kappa = 1/(1 + evol_sigma_t2/obs_sigma_t2), evol_sigma_t2 is the evolution error variance, and obs_sigma_t2 is the observation error variance. In particular, the decision rule is based on the posterior mean of kappa.

Value

A vector (or matrix) of indices identifying the signals according to the horsehoe-type thresholding rule.

Note

The thresholding rule depends on whether the prior variance for the state variable mu (i.e., evol_sigma_t2) is scaled by the observation standard deviation, obs_sigma_t2. Explicitly, if mu[t] ~ N(0, evol_sigma_t2[t]) then the correct thresholding rule is based on kappa = 1/(1 + evol_sigma_t2/obs_sigma_t2). However, if mu[t] ~ N(0, evol_sigma_t2[t]*obs_sigma_t2[t]) then the correct thresholding rule is based on kappa = 1/(1 + evol_sigma_t2). The latter case may be implemented by omitting the input for post_obs_sigma_t2 (or setting it to NULL).

Compute initial Cholesky decomposition for TVP Regression

Description

Computes the Cholesky decomposition for the quadratic term in the (Gaussian) posterior of the TVP regression coefficients. The sparsity pattern will not change during the MCMC, so we can save computation time by computing this up front.

Usage

initCholReg_spam(obs_sigma_t2, evol_sigma_t2, XtX, D = 1)

Arguments

Compute initial Cholesky decomposition for Bayesian Trend Filtering

Description

Computes the Cholesky decomposition for the quadratic term in the (Gaussian) posterior of the Bayesian Trend Filtering coefficients. The sparsity pattern will not change during the MCMC, so we can save computation time by computing this up front.

Usage

initChol_spam(nT, D = 1)

Arguments

Initialize the evolution error variance parameters

Description

Compute initial values for evolution error variance parameters under the dynamic horseshoe prior

Usage

initDHS(omega)

Arguments

Value

List of relevant components: the T x p evolution error SD sigma_wt, the T x p log-volatility ht, the p x 1 log-vol unconditional mean(s) dhs_mean, the p x 1 log-vol AR(1) coefficient(s) dhs_phi, the T x p log-vol innovation SD sigma_eta_t from the PG priors, the p x 1 initial log-vol SD sigma_eta_0, and the mean of log-vol means dhs_mean0 (relevant when p > 1)

Initialize the parameters for the initial state variance

Description

The initial state SDs are assumed to follow half-Cauchy priors, C+(0,A), where the SDs may be common or distinct among the states.

Usage

initEvol0(mu0, commonSD = TRUE)

Arguments

Details

This function initializes the parameters for a PX-Gibbs sampler.

Value

List of relevant components: the p x 1 evolution error SD sigma_w0, the p x 1 parameter-expanded RV's px_sigma_w0, and the corresponding global scale parameters sigma_00 and px_sigma_00 (ignore if commonSD)

Initialize the evolution error variance parameters

Description

Compute initial values for evolution error variance parameters under the various options: dynamic horseshoe prior ('DHS'), horseshoe prior ('HS'), Bayesian lasso ('BL'), normal stochastic volatility ('SV'), or normal-inverse-gamma prior ('NIG').

Usage

initEvolParams(omega, evol_error = "DHS")

Arguments

Value

List of relevant components: sigma_wt, the T x p matrix of evolution standard deviations, and additional parameters associated with the DHS and HS priors.

Initialize the stochastic volatility parameters

Description

Compute initial values for normal stochastic volatility parameters. The model assumes an AR(1) for the log-volatility.

Usage

initSV(omega)

Arguments

Value

List of relevant components: sigma_wt, the T x p matrix of standard deviations, and additional parameters (unconditional mean, AR(1) coefficient, and standard deviation).

Helper function for initializing parameters for ASV model

Description

Helper function for initializing parameters for ASV model

Usage

init_paramsASV(data, evol_error, D)

Arguments

Value

a list containing 4 sets of parameters

• s_p_error_term: matrix containing mean and the variance from 10-componenet gaussian mixture (Omori et al. 2007)

• s_mu: a vector containing the posterior sample of log variance h,

• s_evolParams0: a list containing posterior samples of parameters associated with the variance of first D observation of the log variance term, h.

• s_evolParams: a list containing posterior samples parameters associated with the variance of D to the last observations of the log variance temr , h.

Helper function for initializing parameters for ASV model with a nugget effect

Description

Helper function for initializing parameters for ASV model with a nugget effect

Usage

init_paramsASV_n(data, evol_error, D)

Arguments

Value

a list containing 4 sets of parameters

• s_p_error_term: matrix containing mean and the variance from 10-componenet gaussian mixture (Omori et al. 2007)

• s_mu: a vector containing the posterior sample of log variance h,

• s_evolParams0: a list containing posterior samples of parameters associated with the variance of first D observation of the log variance term, h.

• s_evolParams: a list containing posterior samples parameters associated with the variance of D to the last observations of the log variance temr , h.

Compute the inverse log-odds

Description

Compute the inverse log-odds

Usage

invlogit(x)

Arguments

Value

A scalar or vector of values in (0,1)

Compute the log-odds

Description

Compute the log-odds

Usage

logit(x)

Arguments

Value

A scalar or vector of log-odds

Sample components from a discrete mixture of normals

Description

Sample Z from 1,2,...,k, with P(Z=i) proportional to q_iN(mu_i,sig2_i).

Usage

ncind(y, mu, sig, q)

Arguments

Value

Sample from {1,...,k}

Plot the Bayesian trend filtering fitted values

Description

Plot the BTF posterior mean of the conditional expectation with posterior credible intervals (pointwise and joint), the observed data, and true curves (if known)

Usage

## S3 method for class 'dsp'
plot(x, type, true_values = NULL, t01 = NULL, include_joint_bands = FALSE, ...)

Arguments

Details

The plotting behavior depends on the dimension of the posterior samples stored in x$mcmc_output[[type]]:

• 1D (scalar parameter): A density plot is generated using a histogram with overlaid kernel density estimate. The posterior mean and 95% credible interval are annotated, along with the true value if provided.

• 2D (vector-valued parameter over time): A time-series plot is created, showing the posterior mean and 95% pointwise credible intervals. If include_joint_bands = TRUE and the parameter is among "omega", "mu", "yhat", or "zeta", simultaneous credible bands are also drawn. Optionally, ground truth values (if supplied via true_values) are overlaid as orange dots.

• 3D (parameter array): A sequence of time-series plots is drawn, one for each slice of the third dimension (e.g., different components of a multivariate function). Posterior mean, pointwise intervals, joint bands (when applicable), and ground truth are visualized in the same style as the 2D case. The function pauses after each plot, allowing the user to interactively inspect each one. The x-axis values are given by t01. If not provided, they default to evenly spaced points in [0, 1]. For parameters with temporal differencing (e.g., "evol_sigma_t2"), initial time points used for prior initialization are automatically excluded. If the model includes change point detection (model = "changepoint"), and both omega and r are present in the MCMC output, vertical lines are drawn at the estimated change point locations for plots of "mu", "yhat", or "omega".

Value

No return value, called for side effects

Examples

set.seed(200)
signal = c(rep(0, 50), rep(10, 50))
noise = rep(1, 100)
noise_var = rep(1, 100)
for (k in 2:100){
  noise_var[k] = exp(0.9*log(noise_var[k-1]) + rnorm(1, 0, 0.5))
  noise[k] = rnorm(1, 0, sqrt(noise_var[k])) }

y = signal + noise
model_spec = dsp_spec(family = "gaussian", model = "changepoint",
                      D = 1, useAnom = TRUE, obsSV = "SV")
mcmc_output = dsp_fit(y, model_spec = model_spec, nsave = 500, nburn = 500)
# Estimated posterior mean vs ground truth
plot(mcmc_output, type = "mu", true_values = signal)
# Estimated innovation variance vs ground truth for illustration only
plot(mcmc_output, type = "obs_sigma_t2", true_values = noise^2)

Predict changepoints from the output of ABCO

Description

Predict changepoints from the output of ABCO

Usage

## S3 method for class 'dsp'
predict(object, cp_thres = 0.5, cp_prop = FALSE, ...)

Arguments

Details

The changepoint model uses a thresholding mechanism with a latent indicator variable. This function calculates the proportion of samples where the increment exceeds the threshold.

Value

If cp_prop = FALSE, a numeric vector of indices that correspond to indices of the observed data. If cp_prop = TRUE, a list containing:

- 'cp_t':  a numeric vector of indices that correspond to indices of the observed data.
- 'cp_prop': a numeric vector of length (T - D) with the pointwise proportion of samples where the increment exceeds the threshold.

If no proportions exceed cp_thres, then the vector will be a length 0 integer vector.

Examples

set.seed(200)
signal = c(rep(0, 50), rep(10, 50))
noise = rep(1, 100)
noise_var = rep(1, 100)
for (k in 2:100){
  noise_var[k] = exp(0.9*log(noise_var[k-1]) + rnorm(1, 0, 0.5))
  noise[k] = rnorm(1, 0, sqrt(noise_var[k])) }

y = signal + noise
model_spec = dsp_spec(family = "gaussian", model = "changepoint",
                      D = 1, useAnom = TRUE)
mcmc_output = dsp_fit(y, model_spec = model_spec, nsave = 500, nburn = 500)
predict(mcmc_output)

Sample the AR(1) coefficient(s)

Description

Compute one draw of the AR(1) coefficient in a model with Gaussian innovations and time-dependent innovation variances. In particular, we use the sampler for the log-volatility AR(1) process with the parameter-expanded Polya-Gamma sampler. The sampler also applies to a multivariate case with independent components.

Usage

sampleAR1(h_yc, h_phi, h_sigma_eta_t, prior_dhs_phi = NULL)

Arguments

Value

p x 1 vector of sampled AR(1) coefficient(s)

Note

For the standard AR(1) case, p = 1. However, the function applies more generally for sampling p > 1 independent AR(1) processes (jointly).

Sampler for first or second order random walk (RW) Gaussian dynamic linear model (DLM)

Description

Compute one draw of the T x 1 state variable mu in a DLM using back-band substitution methods. This model is equivalent to the Bayesian trend filtering (BTF) model, assuming appropriate (shrinkage/sparsity) priors for the evolution errors.

Usage

sampleBTF(
  y,
  obs_sigma_t2,
  evol_sigma_t2,
  D = 1,
  loc_obs = NULL,
  chol0 = NULL,
  prior_mean = NULL
)

Arguments

Value

T x 1 vector of simulated states

Note

Missing entries (NAs) are not permitted in y. Imputation schemes are available.

Sampler for first or second order random walk (RW) Gaussian dynamic linear model (DLM)

Description

Compute one draw of the p x 1 B-spline basis coefficients beta in a DLM using back-band substitution methods. The coefficients are penalized with a prior on the D = 0, D = 1, or D = 2 differences. This model is equivalent to the Bayesian trend filtering (BTF) model applied to p x 1 vector of equally-spaced B-spline coefficients, with the basis matrix serving as a design matrix in the observation equation.

Usage

sampleBTF_bspline(
  y,
  X,
  obs_sigma2,
  evol_sigma_t2,
  XtX_bands,
  Xty = NULL,
  D = 1
)

Arguments

Value

p x 1 vector of simulated basis coefficients beta

Note

Missing entries (NAs) are not permitted in y. Imputation schemes are available.

Sampler for first or second order random walk (RW) Gaussian dynamic linear model (DLM)

Description

Compute one draw of the T x p state variable beta in a DLM using back-band substitution methods. This model is equivalent to the Bayesian trend filtering (BTF) model applied to p dynamic regression coefficients corresponding to the design matrix X, assuming appropriate (shrinkage/sparsity) priors for the evolution errors.

Usage

sampleBTF_reg(y, X, obs_sigma_t2, evol_sigma_t2, XtX, D = 1, chol0 = NULL)

Arguments

Value

T x p matrix of simulated dynamic regression coefficients beta

Note

Missing entries (NAs) are not permitted in y. Imputation schemes are available.

(Backfitting) Sampler for first or second order random walk (RW) Gaussian dynamic linear model (DLM)

Description

Compute one draw of the T x p state variable beta in a DLM using back-band substitution methods. This model is equivalent to the Bayesian trend filtering (BTF) model applied to p dynamic regression coefficients corresponding to the design matrix X, assuming appropriate (shrinkage/sparsity) priors for the evolution errors. The sampler here uses a backfitting method that draws each predictor j=1,...,p conditional on the other predictors (and coefficients), which leads to a faster O(Tp) algorithm. However, the MCMC may be less efficient.

Usage

sampleBTF_reg_backfit(y, X, beta, obs_sigma_t2, evol_sigma_t2, D = 1)

Arguments

Value

T x p matrix of simulated dynamic regression coefficients beta

Note

Missing entries (NAs) are not permitted in y. Imputation schemes are available.

Sampler for first or second order random walk (RW) Gaussian dynamic linear model (DLM) with additional shrinkage to zero

Description

Compute one draw of the T x 1 state variable mu in a DLM using back-band substitution methods. This model is equivalent to the Bayesian trend filtering (BTF) model, assuming appropriate (shrinkage/sparsity) priors for the evolution errors, with an additional shrinkage-to-zero prior.

Usage

sampleBTF_sparse(
  y,
  obs_sigma_t2,
  evol_sigma_t2,
  zero_sigma_t2,
  D = 1,
  chol0 = NULL
)

Arguments

Value

T x 1 vector of simulated states

Note

Missing entries (NAs) are not permitted in y. Imputation schemes are available.

Sample the dynamic shrinkage process parameters

Description

Compute one draw for each of the parameters in the dynamic shrinkage process for the special case in which the shrinkage parameter kappa ~ Beta(alpha, beta) with alpha = beta. The primary example is the dynamic horseshoe process with alpha = beta = 1/2.

Usage

sampleDSP(
  omega,
  evolParams,
  sigma_e = 1,
  loc = NULL,
  prior_dhs_phi = c(10, 2),
  alphaPlusBeta = 1
)

Arguments

Value

List of relevant components:

• the T x p evolution error standard deviations sigma_wt,

• the T x p log-volatility ht, the p x 1 log-vol unconditional mean(s) dhs_mean,

• the p x 1 log-vol AR(1) coefficient(s) dhs_phi,

• the T x p log-vol innovation standard deviations sigma_eta_t from the Polya-Gamma priors,

• the p x 1 initial log-vol SD sigma_eta_0,

• and the mean of log-vol means dhs_mean0 (relevant when p > 1)

Note

The priors induced by prior_dhs_phi all imply a stationary (log-) volatility process.

Sample the parameters for the initial state variance

Description

The initial state SDs are assumed to follow half-Cauchy priors, C+(0,A), where the SDs may be common or distinct among the states.

Usage

sampleEvol0(mu0, evolParams0, commonSD = FALSE, A = 1)

Arguments

Details

This function samples the parameters for a PX-Gibbs sampler.

Value

List of relevant components: the p x 1 evolution error SD sigma_w0 and the p x 1 parameter-expanded RV's px_sigma_w0

Sampler evolution error variance parameters

Description

Compute one draw of evolution error variance parameters under the various options:

• dynamic horseshoe prior ('DHS');

• horseshoe prior ('HS');

• normal-inverse-gamma prior ('NIG').

Usage

sampleEvolParams(
  omega,
  evolParams,
  sigma_e = 1,
  evol_error = "DHS",
  loc = NULL
)

Arguments

Value

List of relevant components in evolParams: sigma_wt, the T x p matrix of evolution standard deviations, and additional parameters associated with the DHS and HS priors.

Note

The list evolParams is specific to each evol_error type, but in each case contains the evolution error standard deviations sigma_wt.

To avoid scaling by the observation standard deviation sigma_e, simply use sigma_e = 1 in the functional call.

Sample a Gaussian vector using the fast sampler of BHATTACHARYA et al.

Description

Sample from N(mu, Sigma) where Sigma = solve(crossprod(Phi) + solve(D)) and mu = Sigma*crossprod(Phi, alpha):

Usage

sampleFastGaussian(Phi, Ddiag, alpha)

Arguments

Value

Draw from N(mu, Sigma), which is p x 1, and is computed in O(n^2*p)

Note

Assumes D is diagonal, but extensions are available

Sample the AR(1) unconditional means

Description

Compute one draw of the unconditional means in an AR(1) model with Gaussian innovations and time-dependent innovation variances. In particular, we use the sampler for the log-volatility AR(1) process with the parameter-expanded Polya-Gamma sampler. The sampler also applies to a multivariate case with independent components.

Usage

sampleLogVolMu(h, h_mu, h_phi, h_sigma_eta_t, h_sigma_eta_0, h_log_scale = 0)

Arguments

Value

a list containing

• the sampled mean(s) dhs_mean and

• the sampled precision(s) dhs_mean_prec_j from the Polya-Gamma parameter expansion

Sample the mean of AR(1) unconditional means

Description

Compute one draw of the mean of unconditional means in an AR(1) model with Gaussian innovations and time-dependent innovation variances (for p > 1). More generally, the sampler applies to the "mean" parameter (on the log-scale) for a Polya-Gamma parameter expanded hierarchical model.

Usage

sampleLogVolMu0(h_mu, h_mu0, dhs_mean_prec_j, h_log_scale = 0)

Arguments

Value

The sampled mean parameter dhs_mean0

Note

This sampler is particularly for p > 1 and the setting in which we want hierarchical shrinkage effects, e.g. predictor- and time-dependent shrinkage, predictor-dependent shrinkage, and global shrinkage, with a natural hierarchical ordering.

Sample the latent log-volatilities

Description

Compute one draw of the log-volatilities using a discrete mixture of Gaussians approximation to the likelihood (see Omori, Chib, Shephard, and Nakajima, 2007) where the log-vols are assumed to follow an AR(1) model with time-dependent innovation variances. More generally, the code operates for p independent AR(1) log-vol processes to produce an efficient joint sampler in O(Tp) time.

Usage

sampleLogVols(
  h_y,
  h_prev,
  h_mu,
  h_phi,
  h_sigma_eta_t,
  h_sigma_eta_0,
  loc = NULL
)

Arguments

Value

T x p matrix of simulated log-vols

Note

For Bayesian trend filtering, p = 1. More generally, the sampler allows for p > 1 but assumes (contemporaneous) independence across the log-vols for j = 1,...,p.

Sampler for the stochastic volatility parameters

Description

Compute one draw of the normal stochastic volatility parameters. The model assumes an AR(1) for the log-volatility.

Usage

sampleSVparams(omega, svParams)

Arguments

Value

List of relevant components in svParams: sigma_wt, the T x p matrix of standard deviations, and additional parameters associated with SV model.

Sampler for the stochastic volatility parameters using same functions as DHS prior

Description

Compute one draw of the normal stochastic volatility parameters. The model assumes an AR(1) for the log-volatility.

Usage

sampleSVparams0(omega, svParams)

Arguments

Value

List of relevant components in svParams: sigma_wt, the T x p matrix of standard deviations, and additional parameters associated with SV model.

Sampling from 10-component Gaussian Mixture component described in Omori et al. 2007

Description

Samples from the conditional posterior distribution of the log(\chi^2) distribution by approximating it with the mixture Gaussian distribution described in Omori et al. 2007.

Usage

sample_j_wrap(Td, obs = NULL)

Arguments

Value

Dataframe containing the posterior samples: mean and variance for the mixture component.

Note

When the obs is not not specified, the components are samples from the prior distribution.

Wrapper function for C++ call for sample mat, check pre-conditions to prevent crash

Description

Wrapper function for C++ call for sample mat, check pre-conditions to prevent crash

Usage

sample_mat_c(row_ind, col_ind, mat_val, mat_l, num_inp, linht, rd, D)

Arguments

Compute Simultaneous Band Scores (SimBaS)

Description

Compute simultaneous band scores (SimBaS) from Meyer et al. (2015, Biometrics). SimBaS uses MC(MC) simulations of a function of interest to compute the minimum alpha such that the joint credible bands at the alpha level do not include zero. This quantity is computed for each grid point (or observation point) in the domain of the function.

Usage

simBaS(sampFuns)

Arguments

Value

m x 1 vector of simBaS

Note

The input needs not be curves: the simBaS may be computed for vectors to achieve a multiplicity adjustment.

The minimum of the returned value, PsimBaS_t, over the domain t is the Global Bayesian P-Value (GBPV) for testing whether the function is zero everywhere.

Simulate noisy observations from a dynamic regression model

Description

Simulates data from a time series regression with dynamic regression coefficients. The dynamic regression coefficients are simulated as a Gaussian random walk, where jumps occur with a pre-specified probability sparsity. The coefficients are initialized by a N(0,1) simulation.

Usage

simRegression(
  nT = 200,
  p = 20,
  p_0 = 15,
  sparsity = 0.05,
  RSNR = 5,
  ar1 = 0,
  include_plot = FALSE
)

Arguments

Value

a list containing

• the simulated function y

• the simulated predictors X

• the simulated dynamic regression coefficients beta_true

• the true function mu_true

• the true observation standard deviation sigma_true

Note

The root-signal-to-noise ratio is defined as RSNR = (sd of true function)/(sd of noise).

Simulate noisy observations from a dynamic regression model

Description

Simulates data from a time series regression with dynamic regression coefficients. The dynamic regression coefficients are selected using the options from the simUnivariate() function in the wmtsa package.

Usage

simRegression0(
  signalNames = c("bumps", "blocks"),
  nT = 200,
  RSNR = 10,
  p_0 = 5,
  include_intercept = TRUE,
  scale_all = TRUE,
  include_plot = TRUE,
  ar1 = 0
)

Arguments

Value

a list containing

• the simulated function y

• the simulated predictors X

• the simulated dynamic regression coefficients beta_true

• the true function mu_true

• the true observation standard deviation sigma_true

Note

The number of predictors is p = length(signalNames) + p_0.

The root-signal-to-noise ratio is defined as RSNR = (sd of true function)/(sd of noise).

Generate univariate signals of different type

Description

Using code from the archived wmtsa package

Usage

simUnivariate(name, n = 1024, snr = Inf)

Arguments

Value

A numeric vector the same length as n.

Examples

nms <- c("blocks", "linchirp", "mishmash1", "bumps")
z <- lapply(nms, simUnivariate)

Compute the spectrum of an AR(p) model

Description

Compute the spectrum of an AR(p) model

Usage

spec_dsp(ar_coefs, sigma_e, n.freq = 500)

Arguments

Value

A (n.freq x 2) matrix where the first column is the frequencies and the second column is the spectrum evaluated at that frequency

Summarize DSP MCMC chains

Description

Summarize DSP MCMC chains

Usage

## S3 method for class 'dsp'
summary(object, pars, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)

Arguments

Value

Returns a named list of the same length as pars where within each element of the list is a numeric matrix (vector parameters) or vector (scalar parameters). For matrices, each row is a time point (or dimension) of the parameter and each column is a named summary. The names are accessible with colnames. For vectors (scalar parameters), each element is a named summary.

Examples

set.seed(200)
signal = c(rep(0, 50), rep(10, 50))
noise = rep(1, 100)
noise_var = rep(1, 100)
for (k in 2:100){
  noise_var[k] = exp(0.9*log(noise_var[k-1]) + rnorm(1, 0, 0.5))
  noise[k] = rnorm(1, 0, sqrt(noise_var[k])) }

y = signal + noise
model_spec = dsp_spec(family = "gaussian", model = "changepoint",
                      D = 1, useAnom = TRUE, obsSV = "SV")
mcmc_output = dsp_fit(y, model_spec = model_spec, nsave = 500, nburn = 500)

summary_fit <- summary(mcmc_output)
summary_fit$mu[,"mean"]
summary_fit$evol_sigma_t2[,"mean"]

Initializer for location indices for filling in band-sparse matrix

Description

Create row and column indices for locations of symmetric band-sparse matrix. Starts with the locations of the diagonal, proceed with upper-diagonals, followed by lower-diagonals.

Usage

t_create_loc(len, D)

Arguments

Value

a list containing

• the row indices r and

• the column indices c

Initialize the evolution error variance parameters

Description

Compute initial values for evolution error variance parameters under the dynamic horseshoe prior

Usage

t_initEvolParams_no(y, D, omega)

Arguments

Value

List of relevant components: sigma_wt, the T vector of evolution standard deviations, and additional parameters associated with the DHS priors.

Initialize the anomaly component parameters

Description

Compute initial values for either a horseshoe prior or horseshoe+ prior for the anomaly component.

Usage

t_initEvolZeta_ps(zeta)

Arguments

Value

List of relevant components: sigma_wt, the T vector of standard deviations, and additional parameters for inverse gamma priors (shape and scale).

Initialize the stochastic volatility parameters

Description

Compute initial values for normal stochastic volatility parameters. The model assumes an AR(1) for the log-volatility.

Usage

t_initSV(omega)

Arguments

Value

List of relevant components: sigma_wt, the T vector of standard deviations, and additional parameters (unconditional mean, AR(1) coefficient, and standard deviation).

Sample the TAR(1) coefficients

Description

Compute one draw of the TAR(1) coefficients in a model with Gaussian innovations and time-dependent innovation variances. In particular, we use the sampler for the log-volatility TAR(1) process with the parameter-expanded Polya-Gamma sampler. The sampler also applies to a multivariate case with independent components.

Usage

t_sampleAR1(h_yc, h_phi, h_phi2, h_sigma_eta_t, h_st, prior_dhs_phi = NULL)

Arguments

Value

2 vector of sampled TAR(1) coefficient(s)

Sampler for first or second order random walk (RW) Gaussian dynamic linear model (DLM)

Description

Compute one draw of the T state variable mu in a DLM using back-band substitution methods. This model is equivalent to the Bayesian trend filtering (BTF) model, assuming appropriate (shrinkage/sparsity) priors for the evolution errors.

Usage

t_sampleBTF(y, obs_sigma_t2, evol_sigma_t2, D = 1, loc_obs)

Arguments

Value

T x 1 vector of simulated states

Note

Missing entries (NAs) are not permitted in y. Imputation schemes are available.

Sample the thresholded dynamic shrinkage process parameters

Description

Compute one draw for each of the parameters in the thresholded dynamic shrinkage process for the special case in which the shrinkage parameter kappa ~ Beta(alpha, beta) with alpha = beta = 1/2.

Usage

t_sampleEvolParams(
  omega,
  evolParams,
  D = 1,
  sigma_e = 1,
  lower_b,
  upper_b,
  loc,
  prior_dhs_phi = c(20, 1),
  alphaPlusBeta = 1
)

Arguments

Value

List of relevant components:

• the T evolution error standard deviations sigma_wt,

• the T log-volatility ht,

• the 1 log-vol unconditional mean(s) dhs_mean,

• the 1 log-vol AR(1) coefficient(s) dhs_phi,

• the 1 log-vol correction coefficient(s) dhs_phi2,

• the T log-vol innovation standard deviations sigma_eta_t from the Polya-Gamma priors,

• the 1 initial log-vol SD sigma_eta_0,

• the 1 threshold parameter r

Note

The priors induced by prior_dhs_phi all imply a stationary (log-) volatility process.

Sampler for the anomaly component parameters

Description

Compute one draw of the anomaly component parameters.

Usage

t_sampleEvolZeta_ps(omega, evolParams)

Arguments

Value

List of relevant components in evolParams: sigma_wt, the T vector of standard deviations, and additional parameters for inverse gamma priors (shape and scale).

Sample the TAR(1) unconditional means

Description

Compute one draw of the unconditional means in an TAR(1) model with Gaussian innovations and time-dependent innovation variances. In particular, we use the sampler for the log-volatility TAR(1) process with the parameter-expanded Polya-Gamma sampler. The sampler also applies to a multivariate case with independent components.

Usage

t_sampleLogVolMu(
  h,
  h_mu,
  h_phi,
  h_phi2,
  h_sigma_eta_t,
  h_sigma_eta_0,
  h_st,
  h_log_scale = 0
)

Arguments

Value

the sampled mean(s) dhs_mean

Sample the latent log-volatilities

Description

Compute one draw of the log-volatilities using a discrete mixture of Gaussians approximation to the likelihood (see Omori, Chib, Shephard, and Nakajima, 2007) where the log-vols are assumed to follow an TAR(1) model with time-dependent innovation variances. More generally, the code operates for p independent TAR(1) log-vol processes to produce an efficient joint sampler in O(Tp) time.

Usage

t_sampleLogVols(
  h_y,
  h_prev,
  h_mu,
  h_phi,
  h_phi2,
  h_sigma_eta_t,
  h_sigma_eta_0,
  h_st,
  loc
)

Arguments

Value

T x p vector of simulated log-vols

Note

For Bayesian trend filtering, p = 1. More generally, the sampler allows for p > 1 but assumes (contemporaneous) independence across the log-vols for j = 1,...,p.

Sample the threshold parameter

Description

Compute one draw of the threshold parameter in th TAR(1) model with Gaussian innovations and time-dependent innovation variances. The sampler utilizes metropolis hasting to draw from uniform prior.

Usage

t_sampleR_mh(
  h_yc,
  h_phi,
  h_phi2,
  h_sigma_eta_t,
  h_sigma_eta_0,
  h_st,
  h_r,
  lower_b,
  upper_b,
  omega,
  D
)

Arguments

Value

the sampled threshold value r

Sampler for the stochastic volatility parameters

Description

Compute one draw of the normal stochastic volatility parameters. The model assumes an AR(1) for the log-volatility.

Usage

t_sampleSVparams(omega, svParams)

Arguments

Value

List of relevant components in svParams: sigma_wt, the T vector of standard deviations, and additional parameters associated with SV model.

Univariate Slice Sampler from Neal (2008)

Description

Compute a draw from a univariate distribution using the code provided by Radford M. Neal. The documentation below is also reproduced from Neal (2008).

Usage

uni.slice(x0, g, w = 1, m = Inf, lower = -Inf, upper = +Inf, gx0 = NULL)

Arguments

Value

The point sampled, with its log density attached as an attribute.

Note

The log density function may return -Inf for points outside the support of the distribution. If a lower and/or upper bound is specified for the support, the log density function will not be called outside such limits.