bvarSvFit#

Purpose#

Fit a Bayesian VAR with stochastic volatility and optional SSVS variable selection.

Format#

result = bvarSvFit(y)#

result = bvarSvFit(y, p=1, const=1, ssvs=0, n_draws=5000, n_burn={}, seed=42, overall_tightness=0.2, xreg={}, quiet=0, ctl={})

Parameters:

y (TxM matrix or dataframe) – endogenous variables. If a dataframe, column names are used as variable names.
p (scalar) – Optional keyword, lag order. Default = 1.
const (scalar) – Optional keyword, 1 to include a constant (default), 0 to exclude.
ssvs (scalar) – Optional keyword, set to 1 to enable stochastic search variable selection (SSVS). Default = 0.
n_draws (scalar) – Optional keyword, number of posterior draws to retain after burn-in. Default = 5000.
n_burn (scalar) – Optional keyword, number of burn-in draws to discard. Pass {} to use the default (equal to n_draws).
seed (scalar) – Optional keyword, RNG seed for reproducible posterior draws. Default = 42.
overall_tightness (scalar) – Optional keyword, overall tightness of the Minnesota prior. Controls how much data vs prior matters. Default = 0.2.
xreg (TxK matrix or dataframe) – Optional keyword, exogenous regressors included in every equation. Pass {} (default) for no exogenous variables.
quiet (scalar) – Optional keyword, set to 1 to suppress printed output. Overrides ctl.quiet. Default = 0.

ctl (struct) –

Optional keyword, an instance of a bvarSvControl structure. An instance is initialized by calling bvarSvControlCreate() and the following members can be set:

ctl.p

Scalar, lag order. Default = 1.

ctl.const

Scalar, 1 to include a constant, 0 to exclude. Default = 1.

ctl.b_prior

String, prior type for B coefficients.

`"minnesota"`	Minnesota prior with lambda hyperparameters. (Default)
`"flat"`	Diffuse prior with variance ctl.b_prior_var.
`"horseshoe"`	Horseshoe prior (Carvalho, Polson & Scott 2010) for adaptive shrinkage. Each coefficient gets a local shrinkage λ²_ij and a global τ² with half-Cauchy priors. Replaces Minnesota/SSVS for large systems (m=50+). Reuses ctl.b_prior_var as initial τ².

ctl.overall_tightness

Scalar, overall tightness (Minnesota only). Default = 0.2.

ctl.cross_shrinkage

Scalar, cross-variable shrinkage (Minnesota only). Default = 0.5.

ctl.lag_decay

Scalar, lag decay (Minnesota only). Default = 1.0.

ctl.constant_tightness

Scalar, constant tightness (Minnesota only). Default = 1e5.

ctl.exogenous_tightness

Scalar, exogenous tightness (Minnesota only). Default = 1.0.

ctl.soc_tightness

Scalar, sum-of-coefficients (Minnesota only). 0 = disabled. Default = 0.

ctl.sur_tightness

Scalar, single-unit-root (Minnesota only). 0 = disabled. Default = 0.

ctl.exogenous_scale

Scalar, exogenous regressor tightness (Minnesota only). Default = 1.0.

ctl.ar

Scalar, AR(1) prior mean for own lags (Minnesota only). 1.0 = random walk, 0.0 = white noise. Default = 1.0.

ctl.b_prior_var

Scalar, B prior variance (flat prior only). Default = 10.0.

ctl.sv_mu

Scalar, SV level prior mean. Default = 0.0.

ctl.sv_phi_mean

Scalar, SV persistence prior mean. Default = 0.97.

ctl.sv_phi_std

Scalar, SV persistence prior standard deviation. Default = 0.1.

ctl.sv_sigma2

Scalar, SV innovation variance scale. Default = 0.01.

ctl.ssvs

Scalar, enable SSVS variable selection. 0 = off (default), 1 = on.

ctl.ssvs_c0

Scalar, SSVS spike multiplier. Default = 0.1.

ctl.ssvs_c1

Scalar, SSVS slab multiplier. Default = 10.0.

ctl.ssvs_pi_b

Scalar, prior inclusion probability for B. Default = 0.5.

ctl.ssvs_pi_u

Scalar, prior inclusion probability for U off-diagonals. Default = 0.5.

ctl.ssvs_hierarchical

Scalar, 1 for hierarchical prior on inclusion probability, 0 for fixed. Default = 0.

ctl.n_draws

Scalar, number of posterior draws. Default = 5000.

ctl.n_burn

Scalar, burn-in draws. Default = 5000.

ctl.n_thin

Scalar, thinning interval. Default = 1.

ctl.seed

Scalar, RNG seed. Default = 42.

ctl.n_chains

Scalar, number of MCMC chains. Default = 1.

ctl.parallel

Scalar, 1 for parallel chains, 0 for sequential. Default = 0.

ctl.use_asis

Scalar, 1 to enable ASIS interweaving for SV (Kastner & Fruhwirth-Schnatter 2014). Default = 1.

ctl.sv_keep

String, storage mode for stochastic volatility draws.

`"full"`	Store all draws (default). Requires most memory.
`"last"`	Store only the last draw of log-volatilities per iteration.
`"online"`	Store running moments and a reservoir. Best for large systems.

ctl.reservoir_size

Scalar, reservoir size for sv_keep = "online". Default = 500.

ctl.u_bandwidth

Scalar, band-limited Cholesky U estimation. 0 = full lower-triangular U (default, m(m-1)/2 parameters). k > 0 = only first k off-diagonals per column estimated, rest fixed at zero. Reduces parameters from m(m-1)/2 to m*k. For m=50, k=3: 150 vs 1225 parameters. Note: this changes the model (approximation), not just computation.

ctl.quiet

Scalar, set to 1 to suppress printed output. Default = 0.

Returns:

result (struct) –

An instance of a bvarSvResult structure containing:

result.m	Scalar, number of endogenous variables.
result.p	Scalar, lag order.
result.n_obs	Scalar, effective number of observations (T - p).
result.n_total	Scalar, total number of observations (T).
result.const	Scalar, 1 if a constant was included.
result.var_names	Mx1 string array, variable names.
result.b_prior_type	String, `"minnesota"` or `"flat"`.
result.has_ssvs	Scalar, 1 if SSVS was active.
result.b_mean	Kxm matrix, posterior mean of B.
result.b_sd	Kxm matrix, posterior standard deviation of B.
result.b_lower	Kxm matrix, 16th percentile (lower 68% credible band).
result.b_upper	Kxm matrix, 84th percentile (upper 68% credible band).
result.sv_mu	mx1 vector, posterior mean of SV level parameter per equation.
result.sv_phi	mx1 vector, posterior mean of SV persistence per equation.
result.sv_sigma2	mx1 vector, posterior mean of SV innovation variance per equation.
result.sv_mu_sd	mx1 vector, posterior standard deviation of SV level.
result.sv_phi_sd	mx1 vector, posterior standard deviation of SV persistence.
result.sv_sigma2_sd	mx1 vector, posterior standard deviation of SV innovation variance.
result.pip_b	Kxm matrix, posterior inclusion probabilities for B coefficients (SSVS only).
result.pip_u	(m*(m-1)/2)x1 vector, posterior inclusion probabilities for U off-diagonals (SSVS only).
result.phi_accept_rate	mx1 vector, Metropolis-Hastings acceptance rates for SV persistence. Values between 0.2 and 0.6 indicate good mixing.
result.y	Txm matrix, original data.
result.xreg	TxK matrix, exogenous regressors. Empty matrix if none.
result.n_draws	Scalar, number of retained draws.
result.seed	Scalar, RNG seed used.

Examples#

Default SV-BVAR#

new;
library timeseries;

fname = getGAUSSHome("pkgs/timeseries/examples/data/us_macro_quarterly.csv");
data = loadd(fname);

// Default SV-BVAR(1)
result = bvarSvFit(data);

The summary includes stochastic volatility parameters:

================================================================================
BVAR-SV(1) with Minnesota Prior         Variables:             3
Draws: 5000  Burn: 5000  Thin: 1       Observations:        200
Chains: 1                               Effective obs:       199
================================================================================

Stochastic Volatility Parameters
           mu        phi      sigma2    phi_accept
------------------------------------------------------
GDP      -0.231    0.971     0.0089       0.47
CPI      -1.842    0.984     0.0052       0.41
FFR       0.402    0.962     0.0134       0.52
================================================================================
    ⋮  (Posterior coefficient tables for each equation follow)

Multi-Chain SV-BVAR#

Run 4 parallel chains for convergence diagnostics:

new;
library timeseries;

fname = getGAUSSHome("pkgs/timeseries/examples/data/us_macro_quarterly.csv");
data = loadd(fname);

struct bvarSvControl ctl;
ctl = bvarSvControlCreate();
ctl.n_chains = 4;
ctl.parallel = 1;

result = bvarSvFit(data, p=4, n_draws=10000, n_burn=5000, ctl=ctl);

SV-BVAR with SSVS Variable Selection#

Enable stochastic search variable selection to identify which coefficients are nonzero:

new;
library timeseries;

fname = getGAUSSHome("pkgs/timeseries/examples/data/us_macro_quarterly.csv");
data = loadd(fname);

result = bvarSvFit(data, p=4, ssvs=1);

// Posterior inclusion probabilities
print "B inclusion probabilities:";
print result.pip_b;

// Coefficients with PIP > 0.5
mask = result.pip_b .> 0.5;
print "Number of included coefficients:" sumc(vecr(mask));

Large System with Online Storage#

For large systems where storing all draws is infeasible, use online mode:

new;
library timeseries;

// 20-variable system
fname = getGAUSSHome("pkgs/timeseries/examples/data/largeMacro.dat");
data = loadd(fname);

struct bvarSvControl ctl;
ctl = bvarSvControlCreate();
ctl.sv_keep = "online";
ctl.reservoir_size = 1000;

result = bvarSvFit(data, p=2, n_draws=50000, n_burn=10000, ctl=ctl);

// Posterior means available via streaming moments
print result.b_online_mean;

SV-BVAR with Exogenous Regressors#

new;
library timeseries;

fname = getGAUSSHome("pkgs/timeseries/examples/data/us_macro_quarterly.csv");
y = loadd(fname, "gdp_growth + cpi_inflation + fed_funds");

result = bvarSvFit(y, p=4);

Remarks#

Stochastic volatility model: Each equation \(i\) has a time-varying log-variance \(h_{i,t}\) that follows an AR(1) process:

\[h_{i,t} = \mu_i + \phi_i (h_{i,t-1} - \mu_i) + \sigma_i \eta_{i,t}, \quad \eta_{i,t} \sim N(0, 1)\]

The SV parameters \((\mu_i, \phi_i, \sigma_i^2)\) are estimated per equation. The persistence parameter \(\phi_i\) is drawn via Metropolis-Hastings; the acceptance rate is reported in result.phi_accept_rate. Rates between 0.2 and 0.6 indicate good mixing.

ASIS interweaving: By default, the sampler uses the Ancillarity-Sufficiency Interweaving Strategy (Kastner & Fruhwirth-Schnatter 2014) which alternates between centered and non-centered parameterizations of the SV equation. This dramatically improves mixing, especially when persistence \(\phi_i\) is near 1. Disable with ctl.use_asis = 0.

SSVS variable selection: When ctl.ssvs = 1, a spike-and-slab prior (George & McCulloch 1993) is placed on each coefficient:

\[b_j | \gamma_j \sim (1-\gamma_j) \cdot N(0, \tau_{0,j}^2) + \gamma_j \cdot N(0, \tau_{1,j}^2)\]

where \(\tau_0\) (spike) shrinks coefficients toward zero and \(\tau_1\) (slab) allows nonzero values. The posterior inclusion probability (PIP) in result.pip_b indicates which coefficients the data supports. A PIP > 0.5 corresponds to the median probability model.

Scaling is semiautomatic (George, Sun & Ni 2008): the spike and slab widths are scaled by each coefficient’s OLS standard error, making the prior adaptive to the natural scale.

Storage modes:

sv_keep	Memory	Use case
`"full"`	O(n_draws * T * m)	Small systems (m < 10), need full posterior
`"last"`	O(n_draws * m)	Medium systems, need last h for forecasting
`"online"`	O(reservoir * m)	Large systems (m > 10), production use

Multi-chain inference: When ctl.n_chains > 1, each chain starts from an independent random state. Draws from all chains are pooled before computing posterior summaries. Use multiple chains to assess convergence via split-R-hat (see varDiagnostics()).

Library#

timeseries

Source#

bvar.src