*! version 1.0.1 21aug2026

Title

gvar bayes -- Bayesian estimation of the country models

Syntax

gvar bayes [, options]

options

Prior
prior(string) mn (default), ssvs, ng or hs.

Chain
draws(#) retained sweeps before thinning. Default 1000.

burnin(#) discarded sweeps. Default 1000.

thin(#) keep one sweep in #. Default 1.

seed(string) passed to set seed.

Minnesota, prior(mn)
lambda1(#) overall tightness. Default 0.1.

lambda2(#) cross-variable factor. Default 0.5.

lambda3(#) weakly exogenous factor. Default 0.5.

lambda4(#) deterministic factor. Default 100.

prmean(#) prior mean on the first own lag. Default 1.

Spike and slab, prior(ssvs)
tau0(#) spike sd, in OLS standard errors. Default 0.1.

tau1(#) slab sd, in OLS standard errors. Default 10.

pi(#) prior weight on the SPIKE. Default 0.5.

kappa0(#) spike sd for the covariance factor. Default 0.1.

kappa1(#) slab sd for the covariance factor. Default 10.

qij(#) prior weight on the spike in L. Default 0.5.

Normal-Gamma, prior(ng)
tautheta(#) shape of the local variances. Default 1/ln(M) per unit.

dlambda(#) gamma shape on the global variances. Default 0.01.

elambda(#) gamma rate on the global variances. Default 0.01.

nosampletau hold tau fixed instead of sampling it.

Horseshoe, prior(hs)
bgvarhs reproduce the source's endogenous zeta update. See Remarks.

Variances and reporting
sv stochastic volatility instead of a constant variance.

a1(#) inverse-gamma shape on sigma^2. Default 0.01. Ignored under sv.

b1(#) inverse-gamma rate on sigma^2. Default 0.01. Ignored under sv.

trim(#) discard draws whose largest companion eigenvalue modulus is at or above #. Default 1.05. See Remarks.

noeigentrim keep every draw. BGVAR's eigen=FALSE.

bgvarsv reproduce BGVAR's homoskedastic branch. See gvar_methods.

nosummary suppress the report.

The Horseshoe takes no tuning hyperparameters. That is the point of it: the
half-Cauchy scales are learnt rather than set.

gvar foreign must have run first. gvar bayes replaces
gvar estimate: it is an alternative way to obtain the
country-model coefficients, after which
gvar solve and the whole dynamic-analysis group proceed
unchanged.

Description

gvar bayes estimates every country VARX* by Markov chain Monte Carlo
rather than by maximum likelihood, and stores the retained draws. It follows
BVAR_linear.cpp in the BGVAR package.

Each unit is sampled independently, exactly as BGVAR does. The model is not
stacked until gvar solve runs, and the eigenvalue trim is
a global rule applied after stacking, because a country model can be
perfectly stable while the stacked system is not, and the reverse.

Options

Prior

prior(string) selects the prior. mn is the Minnesota prior of
helper.cpp get_Vminnesota; ssvs is stochastic search variable
selection; ng is the Normal-Gamma global-local shrinkage prior; and
hs is the Horseshoe. All four of BGVAR's priors are available. An
unrecognised string is refused by name rather than silently falling back to
the default.

Chain

draws(#), burnin(#) and thin(#) control the chain.
draws counts sweeps AFTER burn-in and BEFORE thinning, so
draws(1000) thin(2) keeps 500. A draws/thin below one is
refused rather than quietly keeping nothing.

seed(string) is passed straight to set seed. Set it. An MCMC
result that cannot be reproduced is not a result, and
gvar bayes makes no attempt to guess a seed for you.

Normal-Gamma

tautheta(#) is the shape of the local variances. Left empty it takes
BGVAR's per-entity default 1/ln(M), where M is that unit's number
of endogenous variables -- so it genuinely differs across units and cannot be
reported as one number.

The local variances are drawn from a generalised inverse Gaussian with shape
tautheta - 0.5. A negative shape is perfectly proper -- it is the
reciprocal branch of the GIG -- and 1/ln(M) itself falls below 0.5 for
any unit with eight or more endogenous variables, so values below 0.5 are
allowed. Only tautheta(0.5) exactly is refused: that puts the shape at
zero, where the sampler has no normalisable proposal. BGVAR's default rule
never produces it, since 1/ln(M) = 0.5 would need M = e^2.

nosampletau holds tau at tautheta() for the whole chain.
BGVAR samples it by default, by random-walk Metropolis on log tau
with the proposal standard deviation adapted over the first half of burn-in,
so nosampletau is the departure from the source, not the default.

Variances

bgvarsv reproduces BGVAR's homoskedastic branch, which writes a
variance into a slot every other line of the file treats as a log
variance. See gvar_methods and the remark below.

The sampler

Two loops per sweep, per unit:

1. the coefficients, equation by equation, as a GLS draw conditioned on
every other equation through L^-1;

2. the free elements of L, by regressing each equation's residual on the
earlier ones,

so Sigma = L D L' comes out of one sweep. L is unit lower triangular in
every draw: its diagonal is never sampled and its upper triangle never touched.

That direction is not a matter of taste. BVAR_linear.cpp:412 forms the
structural residuals as (Y - X A) L^-1', and :780 takes those as the series
whose variance is D. So D = L^-1 Sigma L^-1' and hence
Sigma = L D L'. utils.R:381 writes it out, and :809 confirms the stored
factor is L rather than its inverse. It is also what makes the sampler's own
whitening work: the coefficient draw multiplies by L^-1, which only yields
independent components if L^-1 Sigma L^-1' is diagonal.

Note that BVAR_linear.cpp also contains a joint version at line 413
which draws the coefficients and L together, conditioning only on the
preceding equations. It sits inside a comment block and is not what runs.
The two are different conditional distributions, both produce plausible
posteriors, and no convergence diagnostic distinguishes them; only reading the
file does.

The priors

All three priors are a rule for filling the prior variance matrix V,
not three different samplers. That is worth knowing before reading the
options: they share step 1 and step 2 above entirely.

Minnesota sets V once, deterministically:

own lag (l1/l)^2

cross lag (l1 l2/l)^2 s_i/s_j

weakly exog (l1 l3/(l+1))^2 s_i/s_j*

deterministic l4 s_i

The weakly exogenous block runs from lag zero -- the contemporaneous
star variables are regressors -- which is why its denominator is (l+1)
and not l.

SSVS refills V every sweep from a spike-and-slab, and reports a
posterior inclusion probability per coefficient. pi() is the prior
weight on the spike, so raising it keeps fewer variables.

Normal-Gamma refills V from a hierarchy with one global variance per
lag and one local variance per coefficient:

theta_ij ~ G(tau, tau lambda^2 / 2) local, drawn as a GIG

lambda^2 ~ G(d + tau k, e + tau sum(theta) prod / 2) global

lambda is a running product over lags, so shrinkage compounds with
lag length: lag 3 is shrunk by lambda1 lambda2 lambda3, not by
lambda3 alone.

The endogenous lags and the weakly exogenous block have separate
hierarchies -- columns 0 and 1 of lambda2_A in the source. Sharing one
would tie the shrinkage on a unit's own lags to the shrinkage on its
trade-weighted foreign variables, which is exactly what the two columns exist
to avoid.

Horseshoe is also global-local, but writes its half-Cauchy tails as a
pair of inverse-gamma steps:

lambda_j | nu_j, tau ~ IG(1, 1/nu_j + a_j^2/(2 tau))

nu_j | lambda_j ~ IG(1, 1 + 1/lambda_j)

tau | zeta, a ~ IG((n+1)/2, 1/zeta + sum(a^2/lambda)/2)

zeta | tau ~ IG(1, 1 + 1/tau)

so it needs no GIG draw at all and takes no tuning hyperparameters. It has
three hierarchies where Normal-Gamma has two: the endogenous lags, the
weakly exogenous block, and the free elements of L each carry their own
tau and zeta. Normal-Gamma handles L through a separate
lambda2_L instead, so the layout does not carry over between them.

Neither ng nor hs touches the deterministic rows. The constant and
any trend keep their flat prior variance, exactly as in the source -- the
hierarchy is not asked to shrink an intercept.

Remarks

Set a seed, and report it. Everything else in this help page is
secondary to that.

On a large system in levels, the default prior mean gives an unstable GVAR.
prmean() sets the prior mean on each variable's own first lag, and its
default of 1 is a random-walk prior -- which is the standard choice for
macro data in levels, and which pushes every one of those roots towards the unit
circle. Stack 26 country models and the assembled system tips over. Measured on
the shipped 26-unit demo with

. gvar bayes, prior(mn) prmean(p) lambda1(l) noeigentrim draws(60) burnin(60) seed(20260812)
. gvar solve

prmean lambda1 largest non-unit |eig| explosive

1 0.10 1.02737 3

0 0.10 0.99856 0

0 0.05 0.99422 0

0 0.02 0.71056 0

0 0.01 0.59987 0

Why "non-unit", and what the plain maximum would tell you. Two of the 272
companion roots sit on the unit circle at every one of those
configurations, so the plain max |eig| is exactly 1.00000 for all four
prmean(0) rows and cannot distinguish them. Those two are the dominant
block's: it is a three-variable VECM of cointegrating rank 1, so it carries
3 - 1 = 2 unit roots, and gvar bayes does not touch it -- the sampler fits
the country models as a VARX in levels, which is also why the 69 unit roots
of the ML fit (K minus the sum of the ranks) do not survive into these rows. The
column above therefore reports the largest modulus more than 1e-8 away from 1.
Roots above 1 are kept, which is why the prmean(1) row reads 1.02737
in both senses; a filter meant to remove unit roots must never hide an explosive
one. explosive counts roots beyond 1 + 1e-8.

So prmean(0) -- shrinking towards white noise rather than towards a random
walk -- is what makes this model usable, and lambda1(0.1) with it lands at
0.998: stable without crushing the dynamics, which lambda1 of 0.02 and 0.01
plainly do. This is why every example above passes prmean(0). The
prmean(1) row is genuinely explosive -- three roots outside the circle, not
a numerical artefact of roots sitting on it.

The default trim discards every draw of a large model in levels, and you
will meet this on the shipped demo. trim() keeps a draw only if the
largest companion eigenvalue modulus of the whole stacked GVAR is below
the cutoff. For the 26-unit demo -- 134 endogenous variables, 268 companion
eigenvalues -- those moduli run about 1.09 to 1.31, so BGVAR's 1.05 default
leaves nothing and the command stops.

That is not a sign of a bad fit. A GVAR estimated in levels has many roots
close to unity, and the maximum of several hundred of them sits a little
above one in almost every draw. BGVAR's 1.05 is calibrated on much smaller
models, and BGVAR stops in the same situation for the same reason -- so the
behaviour here is the source's, deliberately.

When it happens the command prints the minimum, median and maximum of the
moduli. Read them, then either raise trim() past them or pass
noeigentrim to keep every draw. Do not just lengthen the chain: that
changes how many draws there are, not where the moduli sit.

Trimming happens before the posterior mean is taken, so the mean is over
surviving draws only. And it is a global rule applied after stacking,
because a country model can be perfectly stable while the assembled system is
not, and the reverse -- trimming unit by unit would throw away draws the global
model would have kept.

What gvar bayes writes, and what it does not. It fills the reduced-form
country-model slots -- the lag coefficients, the weakly exogenous coefficients,
the intercept, the trend, and the residual covariance -- so
gvar solve, gvar irf,
gvar fevd, gvar hd and
gvar forecast all work from it unchanged. It does
not produce cointegrating vectors, because a VARX in levels imposes no
rank. gvar pp and gvar overid are defined on
those vectors and refuse rather than use whatever an earlier
gvar estimate left behind. gvar solve
likewise skips its K - sum(r) unit-root count, which has nothing to compare
against here.

The homoskedastic log-variance defect. BGVAR's homoskedastic branch
executes cur_sv.fill(sig2) into a slot that is a log variance
everywhere else in the file: it initialises to -3, all four of its uses
are exp(-0.5 Sv), the stochastic-volatility branch fills it from
update_fast_sv, and it is returned raw. The consequence is that the data
are weighted by exp(-sigma^2) rather than 1/sigma^2. The scaling
multiplies both sides of each equation, so a constant cancels from the
likelihood but not against the prior precision -- and for macro data in
logs sigma^2 is of order 0.01, so exp(-0.01) is essentially one and
the prior dominates far more than intended. This package writes
log sigma^2 by default; bgvarsv reproduces the source exactly, the
same way gvar tcdecomp handles the zeta/eta swap.

The GIG sampler behind prior(ng) is not transcribed from BGVAR.
do_rgig1.cpp handles only the degenerate ends and delegates the general
case to the GIGrvg package, a dependency rather than part of the source
tree. So it was written from the algorithm and gated on its own test against
the closed-form GIG moments, which Stata cannot evaluate directly -- there is
no besselk() -- so K_nu(x) is computed from its integral
representation and that quadrature is checked first against the exact
half-integer closed form.

Stochastic volatility. sv replaces the constant variance with

log(eps_t^2 + 1e-40) = h_t + log(chi^2_1)

h_t = mu + phi (h_ - mu) + sigma eta_t

using the 10-component normal mixture of Omori, Chib, Shephard and Nakajima
(2007) for log(chi^2_1), with BGVAR's priors: mu ~ N(0, 100^2),
phi a stretched Beta(25, 1.5) on (-1,1), and
sigma^2 ~ Gamma(1/2, 1/2). The last two are not conjugate, so both get
an independence Metropolis step -- which is what stochvol's
ProposalSigma2::INDEPENDENCE amounts to.

It is not bit-comparable to stochvol, and that is deliberate. BGVAR calls
stochvol::update_fast_sv, the Kastner and Fruehwirth-Schnatter (2014)
sampler with ancillarity-sufficiency interweaving. stochvol is a
dependency, not part of the BGVAR source, so the algorithm cannot be
transcribed -- the same situation as GIGrvg for prior(ng) and
coda for gvar bconv. What is reproduced is the
target: the same model, the same priors, the same mixture. ASIS is an
efficiency device -- it lowers autocorrelation, it does not change the posterior
-- so this sampler converges to the same distribution and may mix more slowly.
Run gvar bconv and look.

Measured on simulated data with a known volatility path, T = 500, true
(mu, phi, sigma) = (-1.00, 0.95, 0.25): the posterior mean path correlates
0.75 with the truth, and the parameters come back at
(-1.27, 0.86, 0.37). sigma is the hardest of the three and carries real
posterior uncertainty at that sample size. phi landing below 0.95 is the
Beta(25, 1.5) prior pulling towards its own mean of about 0.887, not an
error -- if you want a less opinionated phi, that prior is where to look.

What the stacked GVAR uses under sv. Sigma_t varies with t, so one
covariance has to be chosen for the stacking. The source takes the
elementwise median over time (utils.R:381 and :387, where the collapsed
margin is time -- not the mean, and not Sigma at one period). That is
reproduced. Note an elementwise median of positive-definite matrices need not
itself be positive definite; the package checks rather than assumes, so a
Cholesky-based decomposition may refuse and ask for shrink.

sv and bgvarsv are refused together. bgvarsv reproduces
BGVAR's homoskedastic branch and source defect #8, so it has nothing to say
about a stochastic-volatility fit.

The path is available from Mata as gvar_getbSV(i), a
T*M x draws matrix of log variances; colshape(v, T)' rebuilds one
draw's T x M surface. gvar_getbS() stays the per-equation summary.

The Horseshoe's three zeta updates are not the same expression in the
source. Two of them read

zeta = 1/R::rgamma(1, 1/(1 + 1 / tau));

and the endogenous one reads

zeta_A_endo = 1.0/R::rgamma(1, 1 + 1/(1 / tau_A_endo));

1 + 1/(1/tau) is 1 + tau, and it is not wrapped in 1/(...), so
that one block draws its auxiliary with scale 1 + tau where the other two
use 1/(1 + 1/tau) = tau/(1 + tau). At tau = 0.01 that is a rate of
0.99 against 101. Two of the three lines agree with each other and with the
published augmentation; one does not. This package uses the consistent form by
default and bgvarhs reproduces the source, the same arrangement as
bgvarsv for the log-variance defect.

Integer division in the Horseshoe's global shape. (n+1)/2,
(nstar+1)/2 and (v+1)/2 are integer divisions in the source, because
n, nstar and v are all declared int. With n = 50 the
shape is 25, not 25.5. Reproduced rather than silently promoted: it is what
runs, and the difference is too small to justify diverging over.

Examples

Setup, using the shipped 26-unit demo (see gvar_datasets):

. use gvar_demo26, clear

. gvar setup y Dp eq ep r lr, unit(country) time(quarter) global(poil pmat pmetal) gendog(poil=usa pmat=usa pmetal=usa) spec(gvar_demospec)

. gvar weights using gvar_flows, flow(trade) source(partner) destination(home) year(year) years(2009 2011) type(1) map(gvar_demoagg)

. gvar foreign

Minnesota prior, then solve and trace an impulse response as usual. Note
prmean(0) and noeigentrim -- see Remarks for why the defaults
do not work on this particular model:

. gvar bayes, prior(mn) prmean(0) noeigentrim draws(2000) burnin(2000) seed(20260811)

. gvar solve

. gvar irf, shock(usa:r) step(24) graph

Spike and slab, with a stronger prior weight on the spike:

. gvar bayes, prior(ssvs) pi(0.8) prmean(0) noeigentrim draws(2000) burnin(2000) seed(20260811)

Normal-Gamma, letting tautheta take its per-unit default and be
sampled:

. gvar bayes, prior(ng) prmean(0) noeigentrim draws(2000) burnin(2000) seed(20260811)

The same, with the shape pinned and the hierarchy made tighter:

. gvar bayes, prior(ng) tautheta(0.7) nosampletau prmean(0) noeigentrim draws(2000) burnin(2000) seed(20260811)

With stochastic volatility. Check the chains afterwards: the sampler here is not
the interweaved one, so mixing may be slower than stochvol's:

. gvar bayes, prior(mn) prmean(0) sv noeigentrim draws(2000) burnin(2000) seed(20260811)

. gvar bconv

Horseshoe, which takes no tuning hyperparameters at all:

. gvar bayes, prior(hs) prmean(0) noeigentrim draws(2000) burnin(2000) seed(20260811)

Horseshoe reproducing the source's endogenous zeta update, for
comparison:

. gvar bayes, prior(hs) bgvarhs prmean(0) noeigentrim draws(2000) burnin(2000) seed(20260811)

A historical decomposition needs a Cholesky factor of Sigma_zeta, which is
singular on this data under either estimator -- 136 variables, 134 quarters
-- so shrink is required. That is not a Bayesian matter; the ML fit needs
it too:

. gvar hd, variables(usa:y) shrink graph

Stored results

gvar bayes stores the following in r():

r(draws) retained draws per unit, after thinning

r(burnin) discarded sweeps

r(thin) the thinning interval

r(trim) the eigenvalue trim gvar solve will apply

r(prior) the prior, as typed

r(priorno) 1 Minnesota, 2 SSVS, 3 Normal-Gamma, 4 Horseshoe

The draws themselves are held in the model and reached from Mata:
gvar_getbA(i) the coefficients,
gvar_getbL(i) the covariance factor,
gvar_getbS(i) the residual variances, and
gvar_getbG(i) the prior-specific block -- inclusion
indicators under prior(ssvs), local prior variances under
prior(ng) and prior(hs), and unused under prior(mn).
gvar_getbprior() records which.

Source

BGVAR src/BVAR_linear.cpp, src/helper.cpp get_Vminnesota,
src/do_rgig1.cpp, R/utils.R.

Author

Dr Merwan Roudane

merwanroudane920@gmail.com

https://github.com/merwanroudane