"""VAR model specification."""
from collections.abc import Sequence
from typing import TYPE_CHECKING, Any, Literal, Self
import numpy as np
from pydantic import Field, model_validator
from impulso._arviz_compat import InferenceDataLike
from impulso._base import ImpulsoBaseModel
from impulso._posterior import COEFFICIENTS, EXOG_COEFFICIENTS, INTERCEPT
from impulso.data import VARData, _format_names
from impulso.observation import Gaussian, StudentT
from impulso.priors import MinnesotaPrior
from impulso.protocols import ErrorDistribution, Prior, PyMCVolatilityProcess, Sampler
from impulso.sv.spec import StochasticVolatility
from impulso.volatility import Constant
if TYPE_CHECKING:
from impulso.fitted import FittedVAR
_PRIOR_REGISTRY: dict[str, type] = {
"minnesota": MinnesotaPrior,
}
_VOLATILITY_REGISTRY: dict[str, type] = {
"constant": Constant,
"sv": StochasticVolatility,
}
_ERROR_DIST_REGISTRY: dict[str, type] = {
"gaussian": Gaussian,
"student_t": StudentT,
}
# A column whose spread is below this fraction of its own largest absolute value is
# numerically constant even though it passed VARData's exactly-constant check. Using
# its raw standard deviation would inflate the prior towards infinity, so the floor
# substitutes a scale derived from the column's level instead.
_EXOG_SD_FLOOR_FRACTION: float = 1e-3
def _exog_prior_sigma(
endog: np.ndarray,
x_exog: np.ndarray,
scale: float,
exog_names: Sequence[str] | None = None,
) -> np.ndarray:
"""Prior standard deviations for the exogenous coefficients `B_exog`.
The coefficient on an exogenous regressor is not a unit-free quantity: it
converts the regressor's units into the dependent variable's. A prior fixed
in coefficient space therefore encodes a different belief for every dataset
— crushing coefficients on small-scale regressors and leaving coefficients
on large-scale ones effectively unrestricted. This scales the prior so the
belief lives in *contribution* space instead:
sd[i, j] = scale * sigma_i / s_j
where `sigma_i` is the AR(1) residual standard deviation of endogenous
variable `i` (the same scale the Minnesota prior uses) and `s_j` is the
sample standard deviation of exogenous column `j`. One prior standard
deviation of `B_exog[i, j]` then moves variable `i` by `scale` of its own
residual standard deviation when regressor `j` moves by one of its own.
The default `scale` is deliberately loose (see `VAR.exog_prior_scale`).
Args:
endog: Endogenous data of shape `(T, n_vars)`. Passed whole — `sigma_i`
is a property of the series, not of the estimation sample.
x_exog: Exogenous regressor block of shape `(T_eff, n_exog)`, already
trimmed to the rows the likelihood sees.
scale: Multiplier in units of "residual standard deviations of the
dependent variable per standard deviation of the regressor".
exog_names: Optional column names, used only to make the
constant-column error message readable.
Returns:
Array of shape `(n_vars, n_exog)` of prior standard deviations.
Raises:
ValueError: If a column of `x_exog` is exactly constant. `VARData`
rejects columns that are constant over the whole sample, but
trimming the first `n_lags` rows can flatten a column that did
vary — a dummy that only switches inside the initial conditions,
say. What the likelihood then sees is collinear with the
intercept, so the coefficient is not identified; the floor below
would happily hand it a wide prior and hide that.
"""
# Lazy: `_conjugate` imports scipy at module level, and `spec` is on the
# package import path.
from impulso._conjugate import ar1_residual_sd
sigma = ar1_residual_sd(endog)
s = x_exog.std(axis=0, ddof=1)
# Checked before the floor is applied: the floor exists to tame columns with
# tiny-but-real variation, not to manufacture a scale for columns with none.
degenerate = np.flatnonzero(s <= 0.0)
if degenerate.size:
labels = [exog_names[j] if exog_names is not None else f"column {j}" for j in degenerate]
raise ValueError(
f"exog columns are constant over the estimation sample: {_format_names(labels)}. "
"The first n_lags rows are consumed as initial conditions, and what remains of these columns "
"does not vary, so their coefficients are collinear with the intercept and not identified. "
"Drop the columns, or reduce `lags` so the rows that do vary enter the estimation sample."
)
peak = np.abs(x_exog).max(axis=0)
s_eff = np.maximum(s, _EXOG_SD_FLOOR_FRACTION * peak)
return scale * np.outer(sigma, 1.0 / s_eff)
[docs]
class VAR(ImpulsoBaseModel):
"""Immutable VAR model specification.
`VAR` specifies the *reduced-form* model — lag order, coefficient prior,
volatility process, and observation error distribution. Nothing here says
which shock is which: structural meaning is layered on afterwards, by
applying an identification scheme to the `FittedVAR` that `fit` returns.
Attributes:
lags: Fixed lag order (int >= 1) or selection criterion string.
max_lags: Upper bound for automatic selection. Only valid with string lags.
prior: Prior shorthand string or Prior protocol instance.
volatility: Volatility shorthand string or PyMCVolatilityProcess protocol instance.
exog_prior_scale: Tightness of the prior on the exogenous coefficients
`B_exog`, read in contribution space: one prior standard deviation
moves an endogenous variable by this many of its own AR(1) residual
standard deviations when the regressor moves by one of its own. The
default of 100 is deliberately loose — deterministic and exogenous
terms are conventionally left near-uninformative (the conjugate
engine uses `Vc = 10e6` on the intercept), and the prior's job here
is to stop the scale of the regressor from silently setting the
answer, not to shrink. Lower it to shrink `B_exog` towards zero.
Applies only to `VAR.fit`; `prior` governs the lag coefficients.
error_dist: Observation error distribution — shorthand string
(`"gaussian"`, the default, or `"student_t"`) or an
`ErrorDistribution` protocol instance. The string form takes the
adapter's defaults, so `error_dist="student_t"` *infers* the
degrees of freedom; pass `StudentT(nu=5.0)` to fix them. Heavy-
tailed errors are rejected in combination with time-varying
volatility.
Governs the exogenous block only; `prior` governs the lag
coefficients. Both `VAR.fit` and `VAR.prior_predictive` build the
same graph, so it applies to either.
"""
lags: int | Literal["aic", "bic", "hq"] = Field(...)
max_lags: int | None = None
prior: Literal["minnesota"] | Prior = "minnesota"
volatility: Literal["constant", "sv"] | PyMCVolatilityProcess = "constant"
exog_prior_scale: float = Field(100.0, gt=0)
error_dist: Literal["gaussian", "student_t"] | ErrorDistribution = "gaussian"
@model_validator(mode="after")
def _validate_spec(self) -> Self:
if self.max_lags is not None and isinstance(self.lags, int):
raise ValueError("max_lags is only valid when lags is a selection criterion ('aic', 'bic', 'hq')")
if isinstance(self.lags, int) and self.lags < 1:
raise ValueError(f"lags must be >= 1, got {self.lags}")
if self.resolved_error_dist.is_heavy_tailed and self.resolved_volatility.is_time_varying:
raise ValueError(
"Heavy-tailed observation errors are not yet supported with "
"time-varying volatility: the degrees of freedom and the "
"log-volatility innovation variance both absorb outliers, so "
"the two are only weakly identified jointly and NUTS mixes "
"poorly. Use volatility='constant' with error_dist='student_t', "
"or stochastic volatility with Gaussian errors."
)
return self
@property
def resolved_prior(self) -> Prior:
"""Resolve string prior shorthand to a Prior instance."""
if isinstance(self.prior, str):
return _PRIOR_REGISTRY[self.prior]()
return self.prior
@property
def resolved_volatility(self) -> PyMCVolatilityProcess:
"""Resolve string volatility shorthand to a PyMCVolatilityProcess instance."""
if isinstance(self.volatility, str):
return _VOLATILITY_REGISTRY[self.volatility]()
return self.volatility
@property
def resolved_error_dist(self) -> ErrorDistribution:
"""Resolve string error-distribution shorthand to an ErrorDistribution instance."""
if isinstance(self.error_dist, str):
return _ERROR_DIST_REGISTRY[self.error_dist]()
return self.error_dist
@staticmethod
def _default_sampler() -> Sampler:
"""Default sampler for VAR: cores=1 (macOS PyMC segfault), target_accept=0.8."""
from impulso.samplers import NUTSSampler
return NUTSSampler(cores=1, chains=4)
[docs]
def fit(
self,
data: VARData,
sampler: Sampler | None = None,
) -> "FittedVAR":
"""Estimate the Bayesian VAR model.
Args:
data: VARData instance.
sampler: Sampler protocol instance. Defaults to `_default_sampler()`
(`cores=1`, `chains=4`, `target_accept=0.8`). Pass an explicit
`NUTSSampler(cores=n)` to opt into parallel chains.
Returns:
FittedVAR with posterior draws.
"""
from impulso.fitted import FittedVAR
if sampler is None:
sampler = self._default_sampler()
model, n_lags = self._build_pymc_model(data)
# Sample
idata = sampler.sample(model)
return FittedVAR.model_construct(
idata=idata,
n_lags=n_lags,
data=data,
var_names=data.endog_names,
volatility=self.resolved_volatility,
error_dist=self.resolved_error_dist,
pymc_model=model,
)
[docs]
def prior_predictive(
self,
data: VARData,
*,
draws: int = 500,
random_seed: int | np.random.Generator | None = None,
) -> InferenceDataLike:
"""Simulate data from the prior, before seeing the likelihood.
Builds the same PyMC graph `fit` builds and calls
`pymc.sample_prior_predictive` on it, so the prior that gets
simulated is exactly the prior that gets sampled — no hand-rolled
second implementation to drift out of sync.
The simulated `obs` paths are **one-step-ahead given the observed
lags**: for each prior draw, `y_t = c + B x_t^obs (+ B_exog z_t) +
L_t eps_t` where `x_t^obs` stacks the *observed* lags of `data`.
The design matrices are baked into the graph, so this is the prior
predictive of the estimation-sample conditional means, not a
simulated path iterated from initial conditions. That is what
`arviz.plot_ppc(..., group="prior")` expects and what makes the
prior comparable to the data on the same time axis.
Note:
Under `volatility="sv"` the per-variable log-volatility priors
are seeded from the OLS residuals of `data` (see
`StochasticVolatility.build_pymc_latent`), so the "prior" is
mildly data-informed in its scale. The constant-volatility
default is not.
Note:
PyMC returns a single chain, so the `obs` variable has shape
`(1, draws, T - n_lags, n_vars)`.
Args:
data: VARData instance. Anchors the prior simulation on the real
lags (and, if present, the real exogenous regressors), and
fixes the lag order when `lags` is a selection criterion.
draws: Number of prior draws.
random_seed: Seed or Generator passed straight through to
`pymc.sample_prior_predictive`.
Returns:
InferenceData-schema container with `prior` (every latent), `prior_predictive`
(the simulated `obs`, dims `(chain, draw, time, var)`) and
`observed_data` (the realised `obs`) groups.
"""
import pymc as pm
model, _ = self._build_pymc_model(data)
with model:
return pm.sample_prior_predictive(draws=draws, random_seed=random_seed)
def _build_pymc_model(self, data: VARData) -> tuple[Any, int]:
"""Build the PyMC model graph for this specification.
Resolves the lag order (running `select_lag_order` when `lags` is a
criterion string), assembles the design matrices, and registers the
intercept, coefficient, exogenous, volatility and likelihood nodes.
The design matrices are baked into the graph as constants, so the
returned model is tied to `data`.
Shared by `fit` (which samples the graph) and `prior_predictive`
(which draws from it without conditioning on the observations). Every
prior lives here, including the scale-adaptive `B_exog` prior
(`_exog_prior_sigma`), so a prior-predictive check cannot describe a
different model from the one `fit` estimates.
Args:
data: VARData instance.
Returns:
Tuple of the built `pymc.Model` and the resolved lag order. The
model is typed `Any` so that importing `impulso.spec` does not
pull in PyMC — the same reason `FittedVAR.pymc_model` is.
"""
import pymc as pm
from impulso._lag_selection import select_lag_order
# Resolve lags
if isinstance(self.lags, str):
max_lags = self.max_lags or 12
ic = select_lag_order(data, max_lags=max_lags)
n_lags = getattr(ic, self.lags)
else:
n_lags = self.lags
# Build prior arrays
prior = self.resolved_prior
n_vars = data.endog.shape[1]
prior_params = prior.build_priors(n_vars=n_vars, n_lags=n_lags)
# Build data matrices
y = data.endog
Y = y[n_lags:]
X_parts = []
for lag in range(1, n_lags + 1):
X_parts.append(y[n_lags - lag : -lag])
X_lag = np.hstack(X_parts)
X_exog = data.exog[n_lags:] if data.exog is not None else None
# OLS residuals seed per-variable SV priors. Constant-volatility adapters
# ignore `data`; only stochastic adapters use it.
if X_exog is not None:
X_full = np.hstack([np.ones((Y.shape[0], 1)), X_lag, X_exog])
else:
X_full = np.hstack([np.ones((Y.shape[0], 1)), X_lag])
B_ols, *_ = np.linalg.lstsq(X_full, Y, rcond=None)
resid = Y - X_full @ B_ols
# Coordinates make the posterior self-describing: `B` comes back labelled
# by variable and by "L<lag>.<variable>" coefficient instead of positional
# `B_dim_0` / `B_dim_1`. Variable names come from `impulso._posterior` —
# the schema ConjugateVAR constructs against too, so both estimators
# agree. `coeff` is lag-major to mirror the X_lag hstack above.
coords: dict[str, object] = {
"var": data.endog_names,
"var1": data.endog_names,
"var2": data.endog_names,
"coeff": [f"L{lag}.{name}" for lag in range(1, n_lags + 1) for name in data.endog_names],
"time": data.index[n_lags:],
}
if data.exog_names is not None:
coords["exog"] = data.exog_names
# Build PyMC model
with pm.Model(coords=coords) as model:
# Intercept
intercept = pm.Normal(INTERCEPT, mu=0, sigma=1, dims="var")
# VAR coefficients with Minnesota prior
B = pm.Normal(
COEFFICIENTS,
mu=prior_params["B_mu"],
sigma=prior_params["B_sigma"],
dims=("var", "coeff"),
)
# Exogenous coefficients. The prior scales with the data so that it
# encodes the same belief regardless of the units the regressors
# happen to be measured in (#192).
if X_exog is not None:
B_exog = pm.Normal(
EXOG_COEFFICIENTS,
mu=0,
sigma=_exog_prior_sigma(y, X_exog, self.exog_prior_scale, data.exog_names),
dims=("var", "exog"),
)
mu = intercept + pm.math.dot(X_lag, B.T) + pm.math.dot(X_exog, B_exog.T)
else:
mu = intercept + pm.math.dot(X_lag, B.T)
# Volatility process: registers latent vars, returns L (Cholesky factor of Σ_t).
# For constant volatility, L is (n_vars, n_vars) and time-invariant.
# For stochastic volatility, L is (T, n_vars, n_vars) — per-t.
volatility = self.resolved_volatility
L = volatility.build_pymc_latent(n_vars=n_vars, T=Y.shape[0], data=resid)
# Sigma deterministic is only registered for time-invariant L —
# for SV, materialising (T, n, n) per draw is wasteful; users can
# reconstruct per-t Σ via `volatility.cholesky_at(posterior, t)`.
if L.ndim == 2:
pm.Deterministic("Sigma", pm.math.dot(L, L.T), dims=("var1", "var2"))
# Likelihood. The error-distribution seam owns which law is
# registered; PyMC handles batched chol natively either way (for
# 2D L every observation uses the same chol; for 3D L (T, n, n)
# observation t uses chol[t]). Under Student-t errors, L L' is the
# *scale* matrix rather than the covariance — see ADR-0007.
error_dist = self.resolved_error_dist
error_dist.build_likelihood("obs", mu=mu, chol=L, observed=Y, dims=("time", "var"))
return model, n_lags