Package 'dsp'

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

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", "ypred", "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.

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", "ypred", "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.

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", "ypred", "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.

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", "ypred", "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.

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", "ypred", "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.

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", "ypred", "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.

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

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.

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.

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.

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:

MCMC Output Parameters

Shared parameters across wrappers include:

• mu: Conditional mean.

• yhat: Posterior predictive $\hat{y}_t$.

• evol_sigma_t2: Variance of the state variable.

• obs_sigma_t2: Observation variance $\sigma^2_{\epsilon}$ (family = "gaussian").

• dhs_phi, dhs_mean: DHS hyperparameter draws (evol_error = "DHS").

• h: Time-varying log-volatility component, log(obs_sigma_t2) (obsSV = "ASV").

Model-specific parameters include:

• changepoint (Gaussian):

  • omega: $D$-th differenced posterior mean draws, $\mu_t$.

  • zeta: anomaly component (when useAnom = TRUE).

  • zeta_sigma_t2: anomaly variance (when useAnom = TRUE).

  • r: threshold parameter in the thresholded DHS log-volatility dynamics.

• regression (Gaussian):

  • beta: time-varying regression coefficient draws.

• bspline (Gaussian):

  • beta: B-spline basis coefficient draws.

• smoothing (Negative Binomial):

  • r: overdispersion parameter draws.

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)

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)

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])

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.

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.

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.

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)

Description

Compute the design matrix X for AR(p) model

Usage

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

Arguments

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().

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).

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.

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.

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

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

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)

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)

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.

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).

Description

Compute the inverse log-odds

Usage

invlogit(x)

Arguments

Value

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

Description

Compute the log-odds

Usage

logit(x)

Arguments

Value

A scalar or vector of log-odds

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}

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,
  times = NULL,
  y_obs = NULL,
  include_joint_bands = FALSE,
  alpha = 0.05,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  mar = NULL,
  par_args = list(),
  legend = TRUE,
  legend_cex = 1,
  legend_pt_cex = 2,
  nr = NULL,
  nc = NULL,
  cp_thres = 0.5,
  ...
)

Arguments

Details

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

• 1D (scalar parameter): A density-style summary is produced using a histogram with an overlaid kernel density estimate. The posterior mean and equal-tailed $(1-\alpha)100\%$ credible interval are marked, and the true value is added if supplied through true_values.

• 2D (time-varying parameter): A time-series plot is produced showing the posterior mean together with equal-tailed pointwise credible intervals. If include_joint_bands = TRUE and the parameter is one of zeta, omega, ypred, or mu, simultaneous credible bands are also displayed. If true_values is provided, the ground truth is overlaid.

• 3D (multi-parameter time-varying quantity): A multi-panel collection of time-series plots is produced, one panel for each slice of the third dimension. Each panel shows the posterior mean, pointwise credible intervals, optional joint bands when supported, and optional ground-truth values. The panel layout is controlled by nr and nc; if omitted, a near-square layout is chosen automatically.

Axis labels, titles, plotting limits, margins, legend display, and additional graphical settings can be customized through xlab, ylab, main, xlim, ylim, mar, legend, legend_cex, legend_pt_cex, and par_args.

The x-axis values are taken from times. If times is not supplied, evenly spaced points on $[0,1]$ are used. For differenced variance parameters such as "evol_sigma_t2" and "zeta_sigma_t2", the initial time points associated with prior initialization are automatically removed before plotting.

For fitted changepoint models, changepoint annotations may be added when supported by the plotted parameter and the corresponding latent components are present in the MCMC output.

Value

No return value, called for side effects

References

Krivobokova, T., Kneib, T., and Claeskens, G. (2010). Simultaneous confidence bands for penalized spline estimators. Journal of the American Statistical Association, 105(490), 852–863. doi:10.1198/jasa.2010.tm09165

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)

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)

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.

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

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).

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.

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.

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.

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.

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.

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.

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

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.

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

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

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.

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.

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.

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.

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.

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).

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).

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)

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

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"]

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

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.

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).

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).

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)

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.

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.

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).

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

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.

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

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.

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.