Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
14,705 characters · 13 sections · 10 citation commands
srvar-toolkit: A Python Implementation of Shadow-Rate Vector Autoregressions with Stochastic Volatility
Vector autoregressions have anchored macroeconomic forecasting since sims1980macroeconomics demonstrated their power for capturing dynamic relationships among economic variables. Yet when central banks push policy rates to zero---the effective lower bound (ELB)---the VAR framework faces a challenge. The observed rate no longer reflects the true stance of monetary policy. Shadow-rate models resolve this tension by positing a latent rate that can turn negative, with the observed rate censored at the bound wu2016measuring. This paper introduces srvar-toolkit, a Python implementation that makes shadow-rate VAR estimation accessible to researchers who prefer open-source tools.
grammatikopoulos2025forecasting develops the statistical framework we implement here, combining shadow-rate data augmentation with Bayesian shrinkage priors and stochastic volatility. The original MATLAB code serves the econometrics community well, but many applied researchers now work primarily in Python. We bridge that gap. Our toolkit offers a transparent, tested implementation with a clean separation between model specification, prior configuration, and MCMC sampling. Users can fit a basic Bayesian VAR in a few lines of code, then progressively add complexity---Minnesota shrinkage, ELB constraints, time-varying volatility---without changing the core workflow.
The architecture reflects four design commitments. We use immutable configuration objects to prevent subtle bugs from mutable state. We inject random number generators explicitly to ensure reproducibility. We limit core dependencies to NumPy, SciPy, and Pandas, avoiding the complexity of probabilistic programming frameworks; plotting and the configuration-driven CLI are provided via optional extras. We maintain a test suite that pins down numerical behaviour across releases. The result is a toolkit suited to both exploratory analysis and production forecasting pipelines.
(ref) summarises the statistical methodology. (ref) describes the software architecture. (ref) presents an empirical illustration. (ref) discusses limitations and future directions.
This section sketches the statistical models we implement. For derivations and theoretical foundations, see grammatikopoulos2025forecasting.
Consider $N$ macroeconomic variables observed over $T$ periods. The VAR($p$) model expresses each variable as a linear function of its own lags and lags of all other variables in the system. We write the model as
where $y_t$ collects the $N$ variables, $c$ is an intercept vector, and $B_1, \ldots, B_p$ are coefficient matrices. The error covariance $\Sigma$ captures contemporaneous relationships. Stacking observations into matrices $Y$ (responses) and $X$ (lagged values plus intercept), we obtain the compact form $Y = XB + E$ with $B$ a $K \times N$ coefficient matrix and $K = 1 + Np$.
Bayesian estimation requires priors on $(B, \Sigma)$. The conjugate choice is Normal-Inverse-Wishart, where
The matrix-normal distribution $\mathcal{MN}$ places a prior mean $M_0$ on the coefficients with row precision governed by $V_0$ and column covariance inherited from $\Sigma$. The inverse-Wishart prior on $\Sigma$ has degrees of freedom $\nu_0$ and scale matrix $S_0$. Conjugacy means the posterior takes the same NIW form with updated parameters $(M_n, V_n, S_n, \nu_n)$, enabling fast closed-form sampling.
Unrestricted VARs with many variables and lags suffer from proliferating parameters. Minnesota-style shrinkage doan1984forecasting addresses this by tightening the prior around parsimonious specifications. We construct $V_0$ so that own lags receive more prior weight than cross-variable lags, and recent lags receive more weight than distant ones. The key hyperparameter $\lambda_1$ controls overall tightness, while $\lambda_3$ governs lag decay. Larger values of $\lambda_1$ pull coefficients towards the prior mean; smaller values let the data speak. We estimate residual variances from univariate AR regressions to scale the prior appropriately across variables with different units.
When a policy rate hits the ELB, the observed value $y_{t,j}^{\text{obs}}$ no longer equals the latent “shadow rate” $y_{t,j}^*$ that would prevail without the constraint. We model the censoring as
where $b$ is the bound. At periods when the observed rate lies at or near the bound, we treat the shadow rate as a latent variable and sample it from its conditional distribution. Given the VAR parameters and the other variables, the shadow rate follows a truncated normal distribution on $(-\infty, b]$. We draw from this distribution at each MCMC iteration, alternating with draws of the VAR coefficients. This data augmentation approach, detailed in grammatikopoulos2025forecasting, respects the censoring while propagating uncertainty about the latent values into forecasts.
Macroeconomic volatility changes over time. The Great Moderation saw declining variance; the 2008 crisis brought a spike. A constant $\Sigma$ misses this variation. We adopt diagonal stochastic volatility, where each variable has its own time-varying log-variance following a random walk:
with $\eta_{t,i} \sim \mathcal{N}(0,1)$ and $\xi_{t,i} \sim \mathcal{N}(0, \sigma_{\eta,i}^2)$. The log-volatility $h_{t,i}$ evolves smoothly, capturing persistent shifts in uncertainty.
Estimation poses a challenge because the observation equation is nonlinear in $h_{t,i}$. kim1998stochastic show that squaring and taking logs transforms the problem into a linear state-space model contaminated by $\log \chi^2_1$ errors. They approximate this non-Gaussian error with a seven-component mixture of normals, enabling standard Kalman filtering and smoothing. We implement their precision-based sampler, which exploits the banded structure of the state-space precision matrix for efficient draws.
With many potential predictors, some coefficients may be effectively zero. Stochastic search variable selection george1993variable formalises this intuition through spike-and-slab priors. In our implementation, each predictor row of the VAR coefficient matrix (including lagged regressors and, optionally, the intercept) has an inclusion indicator that selects either a tight “spike” variance or a diffuse “slab” variance shared across equations. We sample these row-wise inclusion indicators within the Gibbs sampler, letting the data determine which predictors matter. The intercept can be forced to remain included to avoid shrinking the unconditional mean.
We organise srvar-toolkit into modules that separate concerns cleanly. The srvar.api module exposes the two functions most users need: fit() runs the MCMC sampler and returns posterior draws, while forecast() generates predictive simulations from a fitted model. Configuration lives in srvar.spec, where dataclasses define model structure (ModelSpec), prior parameters (\texttt{PriorSpec}), and sampler settings (\texttt{SamplerConfig}). The statistical workhorses reside in specialised modules: \texttt{srvar.bvar} handles NIW posteriors, \texttt{srvar.elb} implements shadow-rate sampling, \texttt{srvar.sv} provides the KSC mixture sampler for stochastic volatility, and \texttt{srvar.ssvs} manages spike-and-slab variable selection. Data ingestion and transformation utilities live in \texttt{srvar.data}, while \texttt{srvar.plotting} offers visualisations for shadow rates, volatility paths, and forecast fan charts.
This modular structure lets users replace or extend components without touching unrelated code. A researcher who wants a different volatility specification can subclass or replace srvar.sv while keeping the rest of the pipeline intact.
Four principles guide our implementation. First, we make configuration objects immutable. Once you create a ModelSpec, you cannot accidentally mutate it mid-estimation. This prevents a class of subtle bugs where shared mutable state leads to irreproducible results. Second, we inject random number generators explicitly. Every sampling function accepts an rng argument of type numpy.random.Generator. Callers control the seed; the library never touches global state. Third, we limit core dependencies to NumPy, SciPy, and Pandas, avoiding the complexity of probabilistic programming frameworks; plotting and the configuration-driven CLI are provided via optional extras. Fourth, we test aggressively. The test suite exercises core samplers with fixed seeds, checking that outputs match expected values across releases and platforms.
A minimal workflow loads data, configures the model, and runs estimation. The listing below fits a four-lag VAR with Minnesota shrinkage to a three-variable system:
Adding shadow-rate constraints and stochastic volatility requires only minor changes to the model specification:
The ElbSpec names the constrained variable and sets the bound; VolatilitySpec activates time-varying variances. The same fit() and forecast() calls handle the expanded model.
We illustrate the toolkit with a small, fully reproducible example consistent with the repository workflows. The demo uses a two-variable monthly system in which an observed short rate is censored at an ELB and a second macro variable responds to lagged rates. This setup mirrors the shadow-rate mechanism while remaining lightweight enough to run quickly on a laptop.
(ref) plots the posterior median shadow rate alongside the observed censored rate. When the ELB binds, the inferred shadow rate can fall below the bound, capturing accommodation that is not visible in the observed series.
(ref) shows forecast fan charts generated from the combined ELB-SV model. The widening bands reflect both parameter uncertainty and time-varying volatility.
We have introduced srvar-toolkit, a Python package for Bayesian VAR estimation with shadow-rate constraints and stochastic volatility. The software implements the methodology of grammatikopoulos2025forecasting in a transparent, tested framework that separates model specification from sampling machinery.
The current implementation has limitations worth noting. Stochastic volatility is diagonal, so we do not model time-varying correlations between variables. Steady-state VAR parameterisation, useful for imposing long-run equilibrium conditions, remains unimplemented. We plan to add alternative shrinkage priors---Bayesian LASSO, Dirichlet-Laplace---in future releases. Despite these gaps, the toolkit covers the core workflow needed for applied shadow-rate VAR analysis.
The software is actively maintained. Documentation, including installation instructions and worked examples, is available at \url{https://srvar-toolkit.charlesshaw.net/}. We welcome contributions and bug reports through the GitHub repository.