cultivars.multivariate.large_dim.hierarchical#
Hierarchical shrinkage: the tightness is estimated, not asserted.
Every Bayesian VAR result depends on hyperparameters someone chose – the Minnesota tightness above all – and the standard practice of fixing them at folklore values means the reported bands condition on a number nobody defends. Giannone, Lenza and Primiceri (2015) close that gap: treat the hyperparameters \(\theta = (\lambda, \mu, \delta)\) as unknowns with their own hyperpriors and let the conjugate model’s closed-form marginal likelihood carry the data’s opinion about them,
where \(p(y \mid \theta)\) is exactly the number
BVAR already reports
and \(p(B, \Sigma \mid y, \theta)\) is its Normal-inverse-Wishart
posterior. The hyperparameter posterior is one density of at most three
dimensions, and everything downstream follows from standard machinery.
Both of the paper’s readings are offered, behind one surface. Empirical
Bayes (method="empirical") maximizes the hyperparameter posterior and
conditions on the mode: fast, and what most applied work does. The full
hierarchy (method="full", the default) runs an adaptive random-walk
Metropolis chain over \(\log\theta\) from the mode and draws
\((B, \Sigma)\) from the exact conditional posterior at every kept
state, so the retained draws marginalize over the tightness and the
bands stop pretending it was known. On stubbornly informative samples the
two agree; when they disagree, the disagreement is the finding.
Two commitments shape the surface. First, the hierarchy is the paper’s,
stated in its parameterization: the Minnesota tightness \(\lambda\)
and, optionally, the sum-of-coefficients and dummy-initial-observation
loosenesses \(\mu\) and \(\delta\) (larger is looser), under
Gamma hyperpriors with mode 0.2 and standard deviation 0.4 for
\(\lambda\) and mode 1 and standard deviation 1 for \(\mu\) and
\(\delta\); the lag decay, the deterministic looseness and the
persistence prior mean are held fixed as the paper holds them. Second,
what the reported evidence means is said on the result: its
log_marginal_likelihood is \(p(y \mid \hat\theta)\) at the mode,
conditional on it, not the integrated evidence over the hyperprior; and a
hyperparameter the data cannot update – \(\mu\) when the presample
mean is zero and the sum-of-coefficients rows vanish – comes back with
its hyperprior as its posterior rather than with a number that looks
estimated.
Layout. HierarchicalBVAR stores the fixed prior settings and
builds the base model on _BayesianVectorAutoRegressionModel in
_internals at the hyperprior modes; _prior_at maps a
hyperparameter vector to a composed
NormalInverseWishartPrior with the
dummy priors at reciprocal tightness, and fit turns the modes and
standard deviations into Gamma shapes and scales through
_gamma_from_mode from _core._samplers. The base’s
_fit_hierarchical finds the mode by Nelder-Mead in log space, runs
the Metropolis chain with a scalar step tuned toward a 0.3 acceptance
rate in blocks of fifty during burn-in and frozen afterwards, draws
(B, Sigma) through _draw_conjugate at each kept state, and packs
a _VectorHierarchicalFit. HierarchicalBVARResult extends
BVARResult, so the
credible intervals, impulse responses, predictive and stability share
are the conjugate family’s. The fixed-tightness model is
bayesian; the priors the
hierarchy ranges over are in priors.
References
Giannone, D., Lenza, M., & Primiceri, G. E. (2015). Prior selection for vector autoregressions. Review of Economics and Statistics, 97(2), 436-451.
Banbura, M., Giannone, D., & Reichlin, L. (2010). Large Bayesian vector auto regressions. Journal of Applied Econometrics, 25(1), 71-92.
Sims, C. A., & Zha, T. (1998). Bayesian methods for dynamic multivariate models. International Economic Review, 39(4), 949-968.
Example
The tightness the data prefer, under both readings, on a bivariate VAR(1) with persistence well below the random-walk prior mean:
>>> import numpy as np
>>> from cultivars.multivariate.large_dim.bayesian import BVAR
>>> from cultivars.bayes.priors import NormalInverseWishartPrior
>>> rng = np.random.default_rng(0)
>>> y = np.zeros((150, 2))
>>> for t in range(1, 150):
... y[t] = 0.6 * y[t - 1] + rng.standard_normal(2)
>>> y += 5.0
>>> model = HierarchicalBVAR(y, order=1)
>>> empirical = model.fit(method="empirical", n_draws=200, seed=0)
>>> full = model.fit(method="full", n_draws=300, n_burn=300, seed=0)
>>> bool(np.allclose(empirical.hyper_mode, full.hyper_mode))
True
>>> low, mean, high = full.hyperparameter_interval("lambda")
>>> bool(low < empirical.hyper_mode[0] < high)
True
>>> fixed = BVAR(y, order=1, prior=NormalInverseWishartPrior(tightness=0.2)).fit(
... n_draws=10, seed=0
... )
>>> bool(empirical.log_marginal_likelihood >= fixed.log_marginal_likelihood)
True
Classes#
Giannone-Lenza-Primiceri hierarchical shrinkage over the conjugate BVAR. |
|
A fitted hierarchical BVAR: the tightness carries a posterior too. |