BVAR.fit#

BVAR.fit(*, n_draws=1000, seed=None)[source]#

Update the posterior exactly and draw from it.

Runs the closed-form Normal-inverse-Wishart update – the design and the prior’s dummy rows stacked, \(K = \Omega^{-1} + X^\top X\), \(\bar B\), \(\bar S\), \(\bar\nu\) – and the closed-form log marginal likelihood, then draws n_draws independent pairs \((B^{(s)}, \Sigma^{(s)})\) from the posterior: \(\Sigma^{(s)} \sim \mathcal{IW}(\bar S, \bar\nu)\) by Bartlett decomposition, and \(B^{(s)} = \bar B + L^{-\top} Z (\Sigma^{(s)})^{1/2\top}\) with \(L L^\top = K\) factored once, \(Z\) a standard normal (w, k) matrix, and \((\Sigma^{(s)})^{1/2}\) that draw’s lower Cholesky factor. The update formulas are given in full on BVARResult.

Because the posterior is exact, n_draws sets Monte Carlo precision only. The posterior mean beta_mean, the covariance mean sigma_u, the log_marginal_likelihood, and every summary-table number are identical across draw counts and seeds; what tightens with more draws is the bands of credible_interval(), irf_bands(), and forecast(), and the resolution of stable_share. A thousand draws puts a 68% band’s endpoints within about a twentieth of a posterior standard deviation; the cost is one (w, w) Cholesky factorization plus, per draw, an inverse-Wishart draw and two (w, k) matrix products, so a hundred-variable system at four lags takes seconds, not minutes.

Parameters:

  • n_draws (int) – Independent posterior draws of \((B, \Sigma)\) to retain, all of them; there is no burn-in or thinning to subtract. Becomes n_kept.

  • seed (int | Generator | None) – Seed or generator for the draws. An integer gives a fresh generator and a reproducible result; a numpy.random.Generator is used in place and advanced; None draws from operating-system entropy and is not reproducible. The posterior itself does not depend on it.

Return type:

BVARResult

Returns:

The fitted BVARResult – exact posterior-mean blocks, the (n_draws, w, k) and (n_draws, k, k) draws, the log marginal likelihood, and the sample accounting.

Raises:

  • SpecificationError – If n_draws is not positive; if the prior is improper, leaving some coefficient variance infinite or non-positive; or if the prior’s variances do not factor as \(\sigma_i^2 \cdot \omega_j\), which is the case for a MinnesotaPrior with cross_equation != 1 and for IndependentNormalWishartPrior.

  • DimensionError – If the sample has fewer than order + 2 rows, which the prior’s univariate residual scales need.

  • NumericalError – If a posterior scale matrix is not positive definite, which a proper prior rules out except under numerically degenerate data.

Example

>>> import numpy as np
>>> from cultivars.multivariate.large_dim.bayesian import BVAR
>>> rng = np.random.default_rng(0)
>>> y = np.zeros((120, 2))
>>> for t in range(1, 120):
...     y[t] = 0.7 * y[t - 1] + rng.standard_normal(2)
>>> model = BVAR(y, order=1, names=("growth", "inflation"))
>>> res = model.fit(n_draws=200, seed=0)
>>> res.n_kept, res.beta_draws.shape
(200, (200, 3, 2))

The posterior is exact, so the draw count changes only the Monte Carlo resolution of the bands, never the posterior mean or the marginal likelihood:

>>> few = model.fit(n_draws=10, seed=1)
>>> np.allclose(few.beta_mean, res.beta_mean)
True
>>> few.log_marginal_likelihood == res.log_marginal_likelihood
True

An integer seed reproduces the draws; a generator is advanced:

>>> np.allclose(model.fit(n_draws=50, seed=3).beta_draws,
...             model.fit(n_draws=50, seed=3).beta_draws)
True
>>> gen = np.random.default_rng(3)
>>> a = model.fit(n_draws=50, seed=gen).beta_draws
>>> b = model.fit(n_draws=50, seed=gen).beta_draws
>>> np.allclose(a, b)
False
>>> model.fit(n_draws=0)
Traceback (most recent call last):
    ...
cultivars.exceptions.SpecificationError: n_draws must be positive; got 0.

See also

  • BVARResult – what the fit returns, with the update in closed form.

  • prior_replications() and prior_predictive_check() – inspect the prior on the scale of the data before fitting.

  • HierarchicalBVAR – when the tightness should be estimated rather than set.

On this page