*! version 1.0.1 21aug2026

Title

gvar methods -- equations, the step-to-source map, and where the
reference implementations disagree

The model

N units, unit i carrying k_i domestic variables. The global
vector x_t stacks them all and has length K = sum k_i. Link
matrices W_i, of dimension (k_i + k*_i) x K, select each unit's own
variables and build its foreign block:

z_it = W_i x_t = (y_it', y*_it')'

Each foreign variable is a weighted average over the units that own the
matching domestic variable, with the weights renormalised to sum to one. A
weakly exogenous global variable, which no unit owns, gets an unweighted
indicator row instead.

Estimation

The country model is a VECMX* with y* treated as I(1) and weakly
exogenous. Reduced-rank ML follows Johansen, conditional on the foreign
block, with critical values from Pesaran, Shin and Smith (2000) indexed by
the deterministic case, n - r and k. Three deterministic cases are
supported, in the MacKinnon-Haug-Michelis numbering:

case 2 restricted intercept, no trend

case 3 unrestricted intercept, no trend

case 4 unrestricted intercept, restricted trend

In cases 2 and 4 the first row of beta carries the deterministic term,
which is why gvar overid lists _cons or
_trend as the leading row of beta.

Solving

Writing each country model in VARX* levels form and stacking through the link
matrices gives

G0 x_t = h0 + h1 t + sum_l H_l x(t-l) + zeta_t

with G0 = stack(A_i0 W_i) and H_l = stack(A_il W_i). Inverting:

x_t = d0 + d1 t + sum_l F_l x(t-l) + eta_t,
F_l = G0^-1 H_l, eta_t = G0^-1 zeta_t

gvar solve reports the eigenvalues of the companion matrix. A GVAR with
K variables and sum r_i cointegrating relations should have exactly
K - sum r_i unit roots; the command checks that and says so.

zeta and eta: which covariance goes where

This distinction causes more errors than anything else in the model, so it is
worth stating flatly.

zeta_t the stacked country-model residuals. Covariance
Sigma_zeta.

eta_t the reduced-form innovations, G0^-1 zeta_t. Covariance
Sigma_eta = G0^-1 Sigma_zeta G0^-1'.

Impulse responses, variance decompositions and persistence profiles
take Sigma_zeta.
The formula is
Phi_h G0^-1 Sigma_zeta e_j / sqrt(e_j' Sigma_zeta e_j),
and irf.m forms G\Sigma_u internally. Supplying
Sigma_eta applies G0^-1 twice and cross-mixes every shock.

Historical decompositions and forecast error bands take Sigma_eta,
because they work from the reduced form directly.

The Beveridge-Nelson trend/cycle decomposition takes Sigma_eta too - or
rather it cumulates eta, not zeta. See
below.

The singular covariance, and what to do about it

Sigma_zeta is K x K but is built from T observations, so its
rank is at most T. In any GVAR with more variables than periods it is
singular by construction. In the shipped 26-unit demo K = 136,
T = 134 and the rank is 133.

Nothing that needs a Cholesky factor can therefore be computed unaided: no
orthogonalised IRF, no structural GIRF over a large leading block, no
orthogonal FEVD, no historical decomposition, and no bootstrap draw in
orthogonalised space. Generalized responses need no factor and are what the
literature reports for systems this size.

The remedies, both from the Toolbox:

shrink shrink the correlation matrix towards the identity with
the intensity that minimises expected quadratic loss, then rescale by the
original standard deviations, so the variances are untouched. In the demo
this gives lambda* = 0.375 and restores full rank.

lambda(#) set the intensity by hand.

vcov(blockdiag) zero every cross-unit covariance, which imposes
that shocks are uncorrelated across units.

shuffle for the bootstrap only: resample whole date columns
instead of orthogonalised scalars. This needs no factor at all and preserves
the empirical cross-section correlation exactly.

For a model with K > T, shuffle is the better bootstrap. Measured
on the demo over 100 kept replications, it discarded 13 draws as unstable
against 59 for shrinkdraw, and gave narrower bands (mean width 0.0094
against 0.0108).

Dynamic analysis

Generalized IRF (Pesaran and Shin 1998). Invariant to the ordering of
the variables, which is why it is the default: with 136 variables no Cholesky
ordering is defensible. The shocks are correlated by construction, so the
responses do not decompose into orthogonal contributions.

Orthogonalised IRF. Cholesky on the whole system. Depends on the order
of x_t; see gvar describe, order. Note that in a
GVAR the impact matrix is G0^-1 P, not P, so it is not
triangular even though P is.

Structural GIRF. Orthogonalises only a leading block. Which variables
lead is the identifying assumption, set with first() and
vorder(). The block size n0 is derived from that choice, never
typed.

Generalized FEVD does not sum to one across shocks, because the shocks
are correlated. That is a property of the estimator, not a defect; the
reported TOTAL column shows how far from one it falls. Use type(oirf)
for shares that do sum to one.

Persistence profiles (Pesaran and Shin 1996) start at one by
construction and must converge to zero. A relation that does not settle back
is not a long-run relation, whatever the rank test said.

The bootstrap

reps() on gvar irf, gvar pp and
gvar stability runs the full model-level bootstrap of
bootstrap_GVAR.m. Every step of the point estimate is replicated:

1. recentre the country-model residuals by row means

2. draw, either by resampling K x T orthogonalised scalars or, with
shuffle, whole date columns

3. regenerate the global vector from the reduced form, holding the first
pmax observations at their actual values

4. rebuild each unit's domestic and foreign blocks through W_i

5. re-estimate all country models at their own (p, q, case, rank)

6. re-stack, re-solve, discard the replication if unstable, capped at
2B

7. recompute the statistic

The number discarded is reported, and matters: a high discard rate means the
bands are conditioned on a non-random subset of draws.

Quantiles follow the source: (0.05, 0.50, 0.95) for the dynamic objects,
(0.90, 0.95, 0.99) for the stability battery, which is one-sided. Note that a
95th percentile from B draws is the 0.95B-th order statistic; these
statistics are right-skewed, so a small B biases the critical value
down and inflates rejections. Use reps(200) or more.

The Bayesian branch

gvar bayes replaces
gvar estimate: it samples every country VARX* by MCMC
instead of fitting it by reduced-rank ML, then writes the posterior mean into
the same country-model slots, so gvar solve and the whole
dynamic-analysis group proceed unchanged. It follows BVAR_linear.cpp in
BGVAR.

The sampler is two loops per sweep, not one. The coefficients are drawn
equation by equation as a GLS draw conditioned on every other equation
through L^-1; then the free elements of L are drawn by regressing each
equation's residual on the earlier ones. BVAR_linear.cpp also contains a
joint version at line 413 that draws the two together, conditioning only on the
preceding equations -- it sits inside a comment block and is not what runs.
Two different conditional distributions, both producing plausible posteriors, and
no convergence diagnostic distinguishes them.

Sigma = L D L', and the direction is fixed by the source rather than
chosen. 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'. utils.R:381 writes the product out and
:809 confirms the stored factor is L and not its inverse. It is also what
makes the whitening work: the coefficient draw multiplies by L^-1, which
only yields independent components if L^-1 Sigma L^-1' is diagonal.

The four priors are four fill-rules for V, not four samplers. They share
both loops above entirely; all that changes is how the prior variance matrix is
refilled each sweep. Minnesota sets it once and deterministically; SSVS refills
it from a spike and slab; Normal-Gamma and Horseshoe refill it from a
global-local hierarchy. See gvar bayes for each.

Three pieces are not in the BGVAR source tree and had to be written from
their published algorithms, then gated on their own tests rather than against
another implementation:

the generalised inverse Gaussian draw for prior(ng) -- do_rgig1.cpp
delegates to the GIGrvg package. Gated against the closed-form GIG moments,
E[X] = sqrt(chi/psi) K_(l+1)(w)/K_l(w).

Geweke's diagnostic for gvar bconv -- conv.diag
delegates to coda. Gated on the spectral density at zero against its
closed form for an AR(1), on the size of the test under i.i.d. chains, and
on its power against a drifting chain and a random walk.

the stochastic-volatility sampler -- BGVAR calls
stochvol::update_fast_sv. The target is reproduced exactly, the
algorithm is not: the source interweaves (Kastner and
Fruehwirth-Schnatter 2014), which lowers autocorrelation without changing the
posterior, so this sampler converges to the same distribution and may mix more
slowly. Gated on the mixture's moments against log(chi^2_1), on recovery of
a known volatility path, and on not inventing variation in
constant-variance data.

Stability is screened per draw, globally. The eigenvalue trim stacks each
retained draw, takes the largest companion eigenvalue modulus, and discards
draws at or above the cutoff -- 1.05 by default, as BGVAR's
.gvar.stacking.wrapper does, stopping outright below 10 stable draws. It is
a global rule because a country model can be perfectly stable while the
assembled system is not, and the reverse. On the shipped 26-unit demo those
moduli run about 1.09 to 1.31, so the default discards everything and the
command says so with the numbers; see gvar bayes.

What the Bayesian branch does not produce is a cointegrating rank. A VARX
in levels imposes none, so beta and alpha do not exist for it.
gvar pp and gvar overid are defined on those
vectors and refuse rather than use whatever an earlier
gvar estimate left behind, and
gvar solve skips its K - sum(r) unit-root count, which
has nothing to compare against. The reduced-form tools -- irf,
fevd, hd,
forecast -- work from either estimator.

Step to source map

step

link matrices create_linkmatrices.m

weight matrices build_wmat.m, weightmat.m,
update_matrix.m

foreign variables create_foreignvariables.m

unit root tests adf.m, ws.m, unitroot_tests.m

lag selection select_varxlag.m, AIC_SBC.m

cointegration cointegration_test.m, get_rank.m

VECMX* ML mlcoint.m; restricted, mlcoint_r.m

VARX* recovery vecx2varx.m

stacking and solving solve_GVAR.m

weak exogeneity test_weakexogeneity.m,
select_lags_we.m

contemporaneous effects contmpcoeff.m

average correlations avgcorrs.m, corrmat.m

serial correlation F Ftest_rsc.m

multivariate diagnostics GVARX .jb.multi, .pt.multi,
.bgserial, .arch.multi

structural stability kraplob.m, nyblom.m,
schow.m

over-identifying restrictions overid_restr.m,
mlcoint_r.m

Granger causality GVARX .grangerGVAR, vars
causality()

impulse responses irf.m, phi.m

variance decomposition fevd.m

persistence profiles pprofile.m

reordering for SGIRF reorder_GVAR.m

covariance handling transform_varcov.m,
ShrinkageCorrLstar.m

forecasting forecast_GVAR.m, con_forecast_GVAR.m

trend/cycle TCdecomp.m, TC_trend_restr.m

historical decomposition BGVAR hd.R

connectedness BGVAR conn; Diebold and Yilmaz (2014)

bootstrap bootstrap_GVAR.m,
bootstrap_GVAR_ss.m

Where the sources are wrong

Eight defects were found in the reference implementations while porting them.
Where reproducing the source's behaviour is useful -- because published results
depend on it -- the source's version is available behind an option and the
corrected version is the default. Two of the eight are not offered that way, and
the entries say why.

1. gvar.m passes the wrong residual to TCdecomp.
TCdecomp.m documents its third argument as the reduced-form residual and
builds its multiplier from H0\H(:,:,j), but gvar.m:3106 calls it
with zeta, the structural residual. Cumulating zeta against a
reduced-form long-run multiplier leaves C(1)(I - G0) sum eta_s in the
"cycle", which is I(1). Measured on the demo, ADF over all 136 cycle series
finds 136 of 136 stationary using eta against 103 of 136 using
zeta. Default residuals(eta); residuals(zeta) reproduces the
Toolbox.

2. .bgserial omits the error-correction terms.
GVARX's Breusch-Godfrey test takes ylagged <- x$datamat[, -(1:K)], every
regressor of the fitted model. For a VECMX* that includes the r_i
error-correction terms. Passing only the short-run block inflates the
statistic in proportion to r_i / k_i: empirical size of the
Edgerton-Shukur F moved from 0.108 to 0.066 on average, and from 0.248 to
0.116 for Australia, whose rank is 5 of 6.

3. con_forecast_GVAR.m indexes Omega beyond its own dimensions.
It sizes omega_Hbar_tilda at K x H_bar but then indexes it with the
forecast horizon, so it only works when H <= H_bar. Truncating the update
at H_bar would be wrong, because the conditioning still moves later
horizons through the cross-horizon covariance block. Omega is built to
max(H, H_bar) here: same formula, wider index range.

4. hd.R feeds the trend coefficient in unscaled.
HDtrend_big[,nn] <- TT + Fcomp %*% HDtrend_big[,nn-1] supplies a
constant TT every period, but the model's deterministic input at period
nn is d1 * trend_nn. The error grows linearly and is propagated by
the companion matrix. hd.R also starts its initial-condition block one
application of Fcomp short. Both are hidden by its trailing residual
slice, which is defined as data minus the sum of everything else. Corrected,
the leftover falls from 48.7 to 9.4e-13 on a series whose level is 4.8, and
the largest contributor to US output changes. bgvar reproduces
hd.R.

5. BVAR_linear.cpp writes a variance into a slot that holds a LOG variance.
The homoskedastic branch executes cur_sv.fill(sig2) where Sv_draw is a
log variance everywhere else in the file: it initialises to -3, all four of
its uses are exp(-0.5 Sv), the SV branch fills it from
update_fast_sv, and it is returned raw. The commented-out log() on
the next line was stranded by a switch to unsafe_col. 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. bgvarsv reproduces it.

6. The Horseshoe's endogenous zeta update is not the same expression as the
other two. Two of the three read
zeta = 1/rgamma(1, 1/(1 + 1/tau)) and the endogenous one reads
zeta_A_endo = 1.0/rgamma(1, 1 + 1/(1/tau_A_endo)). Since
1 + 1/(1/tau) is 1 + tau and is not wrapped in 1/(...), that block
draws its auxiliary with scale 1 + tau where the others use
tau/(1 + tau) -- a rate of 0.99 against 101 at tau = 0.01. Two of the
three lines agree with each other and with the published augmentation; one does
not. bgvarhs reproduces it.

7. predict.R freezes the trend and starts the state a period early.
.get_companion carries the deterministic block forward with
diag(nd), so a trend stays at its last in-sample value and the drift
contribution is a1*T at every horizon instead of a1*(T+h): a trended
forecast is flat in its trend component. Separately, the state is initialised at
Xn[bigT,], which .mlag fills with y_(T-1)...y_(T-p), so the
first companion multiply returns the fitted value at T and horizon
h is really T+h-1. gvar bforecast reuses the
same routine gvar forecast uses, which starts from the
observed terminal values and advances the trend, so both are fixed by
construction. Not offered as an option: a forecast whose first horizon is an
in-sample fit is not a defensible convention to reproduce.

8. Two source paths disagree about whether Bmu is a variance or an sd.
utils.R:752 passes sd = Bmu to the SV prior while
BVAR_linear.cpp:252 passes sqrt(Bmu) as the sd, i.e. treats Bmu
as a variance. With the default Bmu = 100^2 the two readings differ by a
factor of 10^4 in the prior variance. The C++ path is the one that runs for
BGVAR's sampler, so Bmu is a variance here. Recorded rather than corrected,
because there is nothing to correct -- the source is inconsistent with itself,
and the prior is diffuse enough either way that it changes little in practice.

Weak exogeneity: a known deviation

gvar wetest implements test_weakexogeneity.m and reproduces its
formula, its restriction count and its degrees of freedom, but not the
F values printed on the exogeneity_test sheet of the shipped full
demo. This is stated here rather than buried because it is the one published
quantity in the Toolbox's demo output that this package does not reproduce.

What agrees, against that sheet, for the 26-country demo:

marginal-model orders (p*, q*) 25 of 26 units

degrees of freedom 24 of 26 units

5% critical values the same 24 units, exactly

the 5% verdict per test 198 of 206, 96.1%

F values to 1e-6 relative 0 of 206

So the deviation is systematic in the F values while the dimensions are
right: bra:lr is 3.808566 here against a published 7.940544 at degrees of
freedom that agree exactly (118). The two degrees-of-freedom misses are
chl, which follows from its order miss and so is not independent, and
usa, which is off by one in the other direction and is.

Read the pattern, not the level. 96.1% of the tests reach the same
conclusion at 5% as the published sheet, and the rejection rate is what
Dees, di Mauro, Pesaran and Smith (2007) interpret. Individual F values
should not be quoted against the Toolbox's.

Why this is not fixable from the reference. Computing the test in Python
from the Toolbox's own published beta (ECMS_CVs) and its own
published data (countrydata.xls), with the documented formula, also gives
3.808566. Formula, data, beta, dimensions and restriction count all agree,
so the error-correction terms actually fed to the published run are not
beta'Z1 for the beta that run published. That is a property of the
reference output. Note also that all 26 log-likelihoods match to 2.3e-12, which
pins the cointegrating space, and the F test is invariant to any nonsingular
transformation of the error-correction terms -- so a different beta basis
cannot be the explanation either.

Ruled out, so as not to be re-attempted. The Toolbox specifies the
marginal model's foreign block separately from the country model's weakly
exogenous block (fvflag_we / gvflag_we, gvar.m:1680-1714), and its
on-screen note at that pause offers "include the foreign variable, eps, in all
country models" for DdPS(2007). That facility is available here as
weforeign() on gvar wetest, and it defaults to off --
because it is not what the published sheet was run with. Solving
degfr = T - (1 + r + ls*k + ln*ks_we) against the sheet's own critical
values gives 24 of 26 units with ks_we = ks and none at all with
ks_we = ks + 1; bra settles it alone, at a published 118 against 117
with the extra column. Also ruled out: the star variables (2.7e-15), the weight
matrix (1.7e-16), the F formula, and the critical-value function.

_test62.do keeps every number above under test, so the two move together.

Checks that look convincing and are not

Recorded because each cost real time during development.

A decomposition whose last slice is defined residually cannot be
validated by checking that it sums to the data. That identity holds by
construction however wrong the other blocks are. Check the size of the
residual slice.

A persistence profile starts at one for any covariance matrix, because
the same Sigma sits in its numerator and denominator at h = 0. It cannot
detect a wrong covariance argument. An orthogonalised FEVD summing to one
can.

Comparing the size of two candidate cycles cannot separate them. In
TCdecomp the cycle is the OLS residual on [1, trend], hence
orthogonal to the trend by construction under either residual. Only a
stationarity test discriminates.

A discrepancy that grows with the horizon cannot come from an
initialisation error. Read the shape of the error before forming a
hypothesis about its cause.

svmat defaults to float. Validating a Mata computation against Stata's
regress by round-tripping through svmat truncates to seven
significant digits and looks like a real discrepancy. Use svmat double.

A check that reads the producer's output cannot detect that the consumer
never reads it. gvar bayes sampled correctly and stored its draws, and
41 checks confirmed the sampler by reading those draws directly -- while nothing
copied them into the slots gvar solve reads, so every impulse
response after it was the ML model. At least one check per feature has to travel
the documented user-facing path.

Every property of Sigma except Sigma. The country covariance was built as
inv(L) D inv(L)' where the source says L D L' -- a 31.5% error, live in
three commands. Checks existed for its shape, its symmetry, the ranges of what
depended on it, and how it responded to the prior. All four are invariant to
inverting L. A quantity needs at least one check against an independent
construction of the same quantity -- the same number arrived at another way, not
a property it happens to have.

A suite that samples only the easy region certifies nothing about the hard
one. The GIG sampler passed 17 of 17 while broken for omega below about
1e-4, because the smallest omega any check used was 0.63 -- and the
Normal-Gamma prior routinely visits 1e-3 and below. The range a test covers
is part of what it asserts, and nothing in a green result says what that range
was.

In Stata and Mata a missing value compares GREATER than any number, so
missing > x is true and a comparison involving a missing silently passes.
This produced a green check whose reported quantity printed as . -- caught
only because the check printed the number as well as the verdict. Note the
asymmetry: missing < x is false and fails correctly, so half of all
comparisons are safe and the pattern "we use comparisons everywhere and they
work" holds right up to the one that does not. Guard > and >= with
x < ., and never initialise an extremum to -., which is missing.

rc == 0 means a command ran, not that its answer means anything.
irf, fevd, hd and forecast all returned 0 on a model
gvar solve had just declared unstable with six roots above
unity. Usability needs its own assertion.

Narrowing a false claim until it passes is fitting the test to the output.
An assertion that pD falls as the prior tightens failed; restricting it to a
"sane" range failed again; the claim was simply false -- pD is
concentration times curvature and therefore U-shaped. The third version asserted
something true instead. Two related habits: when a fix lands, grep for what the
documentation still claims -- the Sigma correction left the wrong formula in
four places including the command's own printed output -- and when a class of
defect appears, grep the package for it rather than fixing the instance, which is
how one new display as error inside a summary block turned up sixteen
others.

Author

Dr Merwan Roudane

merwanroudane920@gmail.com

https://github.com/merwanroudane