bvarSvFit#

Purpose#

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

Format#

result = bvarSvFit(y)#
result = bvarSvFit(y, ctl)
Parameters:

  • y (TxM matrix or dataframe) – endogenous variables. If a dataframe, column names are used as variable labels in output. If a matrix, variables are labeled “Y1”, “Y2”, etc.

  • ctl (struct) –

    Optional input, an instance of a bvarSvControl structure. An instance is initialized by calling bvarSvControlCreate() and the following members can be set. Set ctl.xreg for exogenous regressors.

    ctl.p

    Scalar, lag order. Default = 1.

    ctl.include_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.

    ctl.lambda1

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

    ctl.lambda2

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

    ctl.lambda3

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

    ctl.lambda4

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

    ctl.lambda5

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

    ctl.lambda6

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

    ctl.lambda7

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

    ctl.lambda_exo

    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.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.include_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.sigma_mean

mxm matrix, time-averaged posterior mean of \(\Sigma\).

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

Scalar, mean number of included B coefficients across draws (SSVS only).

result.n_included_u

Scalar, mean number of included U elements across draws (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.b_draws

Array of n_draws Kxm matrices, raw posterior draws of B (sv_keep = "full" only).

result.h_last

mxn_draws matrix, last log-volatility state per draw (sv_keep = "last" or "online").

result.b_online_mean

Kxm matrix, running posterior mean of B (sv_keep = "online" only).

result.b_online_var

Kxm matrix, running posterior variance of B (sv_keep = "online" only).

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

Scalar, number of chains used.

result.seed

Scalar, RNG seed used.

Examples#

Default SV-BVAR#

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

// 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
================================================================================
...

Multi-Chain SV-BVAR#

Run 4 parallel chains for convergence diagnostics:

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

ctl = bvarSvControlCreate();
ctl.p = 4;
ctl.n_draws = 10000;
ctl.n_burn = 5000;
ctl.n_chains = 4;
ctl.parallel = 1;

result = bvarSvFit(data, ctl);

SV-BVAR with SSVS Variable Selection#

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

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

ctl = bvarSvControlCreate();
ctl.p = 4;
ctl.ssvs = 1;

result = bvarSvFit(data, ctl);

// 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
data = loadd(getGAUSSHome("pkgs/timeseries/examples/largeMacro.dat"));

ctl = bvarSvControlCreate();
ctl.p = 2;
ctl.sv_keep = "online";
ctl.reservoir_size = 1000;
ctl.n_draws = 50000;
ctl.n_burn = 10000;

result = bvarSvFit(data, ctl);

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

SV-BVAR with Exogenous Regressors#

new;
library timeseries;

y = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"), "gdp + cpi + ffr");
X = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"), "oil");

ctl = bvarSvControlCreate();
ctl.p = 4;
ctl.xreg = X;
ctl.quiet = 1;

result = bvarSvFit(y, ctl);

Model#

The SV-BVAR(p) model extends the standard BVAR by allowing the error covariance to change over time:

\[\begin{split}y_t &= B_1 y_{t-1} + \cdots + B_p y_{t-p} + u + \varepsilon_t \\ \varepsilon_t &= A_t^{-1} D_t^{1/2} e_t, \quad e_t \sim N(0, I_m)\end{split}\]

where \(A_t\) is a lower unitriangular matrix (Cholesky factor of the contemporaneous relationships) and \(D_t = \text{diag}(\exp(h_{1,t}), \ldots, \exp(h_{m,t}))\) contains the time-varying variances.

Each log-variance 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)\]

with stationary initialization \(h_{i,0} \sim N(\mu_i, \sigma_i^2 / (1 - \phi_i^2))\).

The SV parameters \((\mu_i, \phi_i, \sigma_i^2)\) are estimated per equation:

  • \(\mu_i\): level of log-volatility (prior: \(N(0, 100)\))

  • \(\phi_i\): persistence (prior: \((\phi_i + 1)/2 \sim \text{Beta}(20, 1.5)\), centering mass near 1)

  • \(\sigma_i^2\): volatility of volatility (prior: \(IG(0.5, 0.5)\))

Algorithm#

The sampler is a multi-block Gibbs sampler (Primiceri 2005; Kastner & Fruhwirth-Schnatter 2014):

  1. Draw B | y, A, h: Equation-by-equation WLS with weights \(\exp(-h_{i,t})\).

  2. Draw A | y, B, h: Column-by-column regression for Cholesky off-diagonals.

  3. Draw h | y, B, A: Per equation, using the Kim-Shephard-Chib (2004) 10-component mixture approximation to \(\log \chi^2(1)\), solved via the precision sampler (McCausland, Miller & Pelletier 2011) in \(O(T)\) time.

  4. Draw \((\mu, \phi, \sigma^2)\) | h: Conjugate draws for \(\mu\) and \(\sigma^2\); Metropolis-Hastings with Laplace proposal for \(\phi\).

When ASIS is enabled (default), steps 3-4 are repeated in the non-centered parameterization, which improves effective sample size by 4-5x (Kastner & Fruhwirth-Schnatter 2014).

Complexity: \(O(T m K^2)\) per iteration (dominated by the WLS draws for B). With m=3, p=4, T=200, 10K draws: typical wall-clock time is 1-2 seconds. 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 varDiagnose()).

Troubleshooting#

Low phi acceptance rates (< 0.2): The Metropolis-Hastings proposal for \(\phi_i\) is too far from the posterior. This usually means the data strongly favors either very high or very low persistence. Solutions: run more burn-in (increase n_burn), or increase draws and thin (set n_thin = 5). ASIS (enabled by default) already helps substantially.

Phi acceptance rates near 1.0 (> 0.95): The prior dominates the likelihood — the data has little information about persistence. This is common with short samples (T < 100) or near-constant volatility series. Consider whether SV is necessary; a constant-volatility bvarFit() may suffice.

R-hat > 1.05 on some parameters: The chain has not converged. Increase n_burn and n_draws. For persistent volatility (\(\phi_i > 0.95\)), 20K+ draws may be needed. Running multiple chains (n_chains = 4) helps diagnose the problem.

SV parameters look unreasonable: If \(\mu_i\) is very large (> 5) or \(\sigma_i^2\) is very small (< 0.001), the model may be overfitting volatility to outliers. Check your data for measurement errors or structural breaks. Verification ————

bvarSvFit has been verified through 30 cross-validation tests covering the core SV sampler and the full SV-BVAR system:

R ``stochvol`` (univariate KSC sampler): The KSC 10-component mixture approximation is the same algorithm used by Kastner’s R stochvol package. Verified on multiple DGPs with known parameters.

R ``bayesianVARs`` (Gruber & Kastner 2023): The full SV-BVAR sampler validated on multivariate data.

Canonical test cases:

  • Clark (2011): 3-variable quarterly macro, tests volatility break detection

  • CCM (2019): tests posterior concentration around known DGP parameters

  • GLP (2015): tests interaction between SV and hierarchical hyperparameters

Real data: FRED-MD large macro dataset, both with and without structural breaks.

All 30 tests pass. See gausslib-var/tests/sv_crossval.rs and the Verification and Cross-Validation page for the full chain of trust. References ———-

  • George, E.I. and R.E. McCulloch (1993). “Variable selection via Gibbs sampling.” Journal of the American Statistical Association, 88(423), 881-889.

  • George, E.I., D. Sun, and S. Ni (2008). “Bayesian stochastic search for VAR model restrictions.” Journal of Econometrics, 142(1), 553-580.

  • Gruber, L. and G. Kastner (2023). “Forecasting macroeconomic data with Bayesian VARs: Sparse or dense? It depends!” Journal of Applied Econometrics, 38(4), 459-482.

  • Kastner, G. and S. Fruhwirth-Schnatter (2014). “Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models.” Computational Statistics & Data Analysis, 76, 408-423.

  • Kim, S., N. Shephard, and S. Chib (1998). “Stochastic volatility: Likelihood inference and comparison with ARCH models.” Review of Economic Studies, 65(3), 361-393.

  • McCausland, W.J., S. Miller, and D. Pelletier (2011). “Simulation smoothing for state-space models: A computational efficiency analysis.” Computational Statistics & Data Analysis, 55(1), 199-212.

  • Primiceri, G.E. (2005). “Time varying structural vector autoregressions and monetary policy.” Review of Economic Studies, 72(3), 821-852.

Library#

timeseries

Source#

bvar.src

See also

Functions bvarSvControlCreate(), bvarFit(), bvarSvForecast(), irfSvCompute(), varDiagnose()

See also

Guides Choosing a VAR Model, Verification and Cross-Validation

On this page