Heavy-Tailed Observation Errors#

Macroeconomic samples contain observations that a Gaussian VAR cannot accommodate: 2020Q2, the 2008 collapse, a devaluation, a data revision. Under Gaussian errors a single such quarter pulls the coefficient estimates and inflates the estimated covariance for the whole sample. Student-t observation errors downweight it automatically, inside the model, with no dummies to pick by hand.

When to reach for it#

The residuals from a Gaussian fit have a few very large values rather than a uniformly wide spread.
You would otherwise be adding outlier dummies, and would rather not choose the dates yourself.
The estimated innovation covariance looks implausibly large relative to the bulk of the sample.

If instead the volatility level drifts over the sample — quiet decades and turbulent ones — that is a stochastic-volatility problem, not a heavy-tails problem. The two cannot currently be combined (see Limitations).

Fitting#

The string shorthand infers the degrees of freedom from the data:

from impulso import VAR

fitted = VAR(lags=4, error_dist="student_t").fit(data)

Fixing them instead is the robust choice on short samples, where nu is only weakly identified:

from impulso import VAR, StudentT

fitted = VAR(lags=4, error_dist=StudentT(nu=5.0)).fit(data)

nu must be strictly greater than 2 — below that the t has infinite variance and forecast bands, variance decompositions and the innovation covariance all stop being defined. Values around 4–6 are aggressively robust; above roughly 30 the fit is indistinguishable from Gaussian.

Reading the posterior for nu#

The degrees of freedom land in the posterior as nu under both parameterisations, so the same code reads them either way:

import arviz as az

az.summary(fitted.idata, var_names=["nu"])

Under inference the free parameter is nu_excess and nu = 2 + nu_excess. The shift means the prior has zero density at the boundary, so the sampler is never dragged toward the infinite-variance edge. A posterior median below about 10 says the data genuinely want heavy tails; a median that has drifted up toward the prior mean (22 by default) says the sample carries little information about the tail, and a fixed nu is the more honest specification.

If you see divergences

nu and the innovation scale trade off against each other, which can make the posterior geometry awkward on short samples. Raise target_accept to 0.9, or fix nu rather than inferring it.

What changes, and what does not#

Under the t, Ω = L Lᵀ is the scale matrix, not the covariance. The distinction propagates predictably:

Quantity	Effect
`fitted.sigma()`	Returns Ω unchanged — the scale matrix, no longer the covariance
`fitted.innovation_covariance()`	New accessor; returns `nu/(nu−2)·Ω`, the actual second moment
`fitted.forecast()`	Innovations follow the t, so bands have fatter tails at the same interquartile width
`identified.impulse_response()`	A “unit shock” is one scale unit = `sqrt((nu−2)/nu)` unconditional sd
`identified.fevd()`	Exactly unchanged — shares are ratios, and the scale cancels
`identified.historical_decomposition()`	Exactly unchanged — shocks are backed out and re-propagated
`identified.counterfactual()`	Exactly unchanged for zero edits; non-zero `ShockPath` values are in scale units
`fitted.conditional_forecast()`	Raises `NotImplementedError`
`identified.structural_scenario()`	Raises `NotImplementedError`

To restore the one-standard-deviation shock convention for impulse responses, scale them by sqrt(nu/(nu−2)):

import numpy as np

nu = fitted.idata.posterior["nu"].values          # (chain, draw)
irf = identified.impulse_response(horizon=20)
factor = np.sqrt(nu / (nu - 2.0))[:, :, None, None, None]
sd_units = irf.idata.posterior_predictive["irf"].values * factor

Is it worth it? Comparing against the Gaussian fit#

Both fits carry a pointwise log-likelihood, so ArviZ compares them directly:

import arviz as az

from impulso import VAR, NUTSSampler

sampler = NUTSSampler(nuts_sampler="pymc")
gaussian = VAR(lags=4).fit(data, sampler=sampler)
student = VAR(lags=4, error_dist="student_t").fit(data, sampler=sampler)

az.compare({"gaussian": gaussian.idata, "student_t": student.idata})

Sampler backend

nuts_sampler="pymc" is required here. The nutpie backend ignores the idata_kwargs that request the log_likelihood group, so a nutpie-sampled InferenceData has nothing for az.compare to read. This is a pre-existing quirk of the nutpie integration, not specific to Student-t errors.

Limitations#

Two combinations are refused outright rather than approximated.

Stochastic volatility. VAR(volatility="sv", error_dist="student_t") raises at construction. The degrees of freedom and the log-volatility innovation variance both absorb outliers, so the two are weakly identified jointly and NUTS mixes poorly. Pick the mechanism that matches the problem: drifting volatility level → stochastic volatility with Gaussian errors; isolated extreme observations → constant volatility with Student-t errors.

Conditional forecasts and structural scenarios. Both raise NotImplementedError under heavy tails. The Waggoner–Zha constrained draw and the ADPRR three-way partition are Gaussian conditional-law results, and the plausibility statistic’s χ² reference assumes Gaussian shocks. Under a t the conditional law has updated degrees of freedom and a Mahalanobis-inflated scale, and the plausibility reference becomes an F/Hotelling statistic. Returning a half-valid scenario would be worse than the error, because nothing in the output would look wrong. Unconditional density forecasts via forecast() and in-sample counterfactual() are both fully valid under the t.

See ADR-0007 for the full derivation of which quantities change and which are exactly invariant.