btf_sparse: Run the MCMC for sparse Bayesian trend filtering

In drkowal/dsp: Dynamic Shrinkage Processes

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,
  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"),
  computeDIC = TRUE,
  verbose = TRUE
)

Arguments

y	the T x 1 vector of time series observations
evol_error	the evolution error distribution; must be one of 'DHS' (dynamic horseshoe prior), 'HS' (horseshoe prior), 'BL' (Bayesian lasso), or 'NIG' (normal-inverse-gamma prior)
zero_error	the shrinkage-to-zero distribution; must be one of 'DHS' (dynamic horseshoe prior), 'HS' (horseshoe prior), 'BL' (Bayesian lasso), or 'NIG' (normal-inverse-gamma prior)
D	degree of differencing (D = 1, or D = 2)
nsave	number of MCMC iterations to record
nburn	number of MCMC iterations to discard (burin-in)
nskip	number of MCMC iterations to skip between saving iterations, i.e., save every (nskip + 1)th draw
mcmc_params	named list of parameters for which we store the MCMC output; must be one or more of: "mu" (conditional mean) "yhat" (posterior predictive distribution) "evol_sigma_t2" (evolution error variance) "zero_sigma_t2" (shrink-to-zero error variance) "obs_sigma_t2" (observation error variance) "dhs_phi" (DHS AR(1) coefficient for evolution error) "dhs_mean" (DHS AR(1) unconditional mean for evolution error) "dhs_phi_zero" (DHS AR(1) coefficient for shrink-to-zero error) "dhs_mean_zero" (DHS AR(1) unconditional mean for shrink-to-zero error)
computeDIC	logical; if TRUE, compute the deviance information criterion DIC and the effective number of parameters p_d
verbose	logical; should R report extra information on progress?

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.

Examples

## Not run: 
# Example 1: Bumps Data
simdata = simUnivariate(signalName = "bumps", T = 128, RSNR = 7, include_plot = TRUE)
y = simdata$y

out = btf_sparse(y)
plot_fitted(y, mu = colMeans(out$mu), postY = out$yhat, y_true = simdata$y_true)

# Example 2: Doppler Data; longer series, more noise
simdata = simUnivariate(signalName = "doppler", T = 500, RSNR = 5, include_plot = TRUE)
y = simdata$y

out = btf_sparse(y)
plot_fitted(y, mu = colMeans(out$mu), postY = out$yhat, y_true = simdata$y_true)

# And examine the AR(1) parameters for the log-volatility w/ traceplots:
plot(as.ts(out$dhs_phi)) # AR(1) coefficient
plot(as.ts(out$dhs_mean)) # Unconditional mean

# Example 3: Blocks data (locally constant)
simdata = simUnivariate(signalName = "blocks", T = 1000, RSNR = 3, include_plot = TRUE)
y = simdata$y

out = btf_sparse(y, D = 1) # try D = 1 to approximate the locally constant behavior
plot_fitted(y, mu = colMeans(out$mu), postY = out$yhat, y_true = simdata$y_true)

## End(Not run)

drkowal/dsp documentation built on July 19, 2023, 11:42 a.m.