2. Reduced-form VAR Modeling
2.1. Description
with \(r_{t} = 1, 2, \dots, h\) and transition probabilities \(p_{r_{t}, r_{t+1}}(I_{t})\).
The modern toolbox builds the reduced-form VAR via the
rfvar_model factory; the option surface uses named arguments
matching the modern arguments blocks.
2.2. Quick-start example
A constant-parameter five-variable VAR on Norwegian quarterly data.
2.2.1. Collecting and transforming data
% CLVMNACSCAB1GQNO : GDP Norway
% IR3TIB01NOQ156N : 3-month interbank rate
% NORCPGRLE01IXOBQ : CPI excluding food and energy
% CCUSSP01NOQ650N : NOK / USD spot
% POILWTIUSDQ : Global price of WTI crude
xrange = '1990Q1:2022Q3';
rawdb = fetch_fred({'NORCPGRLE01IXOBQ','IR3TIB01NOQ156N','CLVMNACSCAB1GQNO', ...
'CCUSSP01NOQ650N','POILWTIUSDQ'});
raw.P = rawdb(1).series(xrange);
raw.INTRATE = rawdb(2).series(xrange);
raw.Y = rawdb(3).series(xrange);
raw.EXRATE = rawdb(4).series(xrange);
raw.POIL = rawdb(5).series(xrange);
db = struct();
db.PAI = raw.P / lag(raw.P, 1);
db.R = 1 + raw.INTRATE/100;
db.GROWTH = raw.Y / lag(raw.Y, 1);
db.EXRATE = 1 / raw.EXRATE; % rising => NOK depreciation
db.PAIOIL = raw.POIL / lag(raw.POIL, 1);
2.2.2. Setting up the reduced-form VAR
endog = {'PAIOIL','GROWTH','PAI','R','EXRATE'};
exog = {};
nlags = 4;
const = true;
mdl = rfvar_model(endog, ...
lag_length = nlags, ...
constant_term = const, ...
deterministic_vars = exog);
(The legacy constructor was rfvar(endog, exog, nlags, const)
positionally; the modern factory uses named arguments and a
slightly different name order.)
2.2.3. Estimating the VAR (classical / OLS)
mdlest = estimate(mdl, ...
data = db, ...
estim_start_date = date2serial(db.GROWTH.start), ...
estim_end_date = date2serial(db.GROWTH.finish));
Estimation dates are passed through date2serial – the modern
date validator does not accept char strings.
2.2.4. Restrictions on the VAR
“Domestic variables do not affect oil prices” via linear restrictions on the VAR coefficients:
linres = {
'b1(PAIOIL,PAI)=0'
'b1(PAIOIL,GROWTH)=0'
'b1(PAIOIL,R)=0'
'b1(PAIOIL,EXRATE)=0'
'b2(PAIOIL,PAI)=0'
'b2(PAIOIL,GROWTH)=0'
'b2(PAIOIL,R)=0'
'b2(PAIOIL,EXRATE)=0'
};
Or programmatically:
linres = cell(0, 1);
for ilag = 1:nlags
for iv = 2:numel(endog)
y = endog{iv};
linres{end+1, 1} = sprintf('b%0.0f(PAIOIL,%s)=0', ilag, y);
end
end
Estimate the restricted VAR:
mdlest_restr = estimate(mdl, ...
data = db, ...
estim_start_date = date2serial(db.GROWTH.start), ...
estim_end_date = date2serial(db.GROWTH.finish), ...
estim_linear_restrictions = linres);
2.3. Identification
A mixture of sign restrictions, contemporaneous zeros, and arbitrary restrictions on the impact matrix:
shock_names = {'oilp','demand','costpush','mp','forex'};
ident_restr1 = {
% normalization with sign restrictions
'PAIOIL{0}@oilp','+'
'GROWTH{0}@demand','+'
'PAI{0}@costpush','+'
'R{0}@mp','+'
'EXRATE{0}@forex','+'
% block: oil-price domestic neutrality
'PAIOIL{0}@demand',0
'PAIOIL{0}@costpush',0
'PAIOIL{0}@mp',0
'PAIOIL{0}@forex',0
% block: domestic-domestic
'GROWTH{0}@costpush',0
'GROWTH{0}@mp',0
'GROWTH{0}@forex',0
'PAI{0}@mp',0
'PAI{0}@forex',0
'R{0}@forex',0
};
agnostic = true;
max_trials = 6000;
Rfunc = struct(); ident = struct();
[Rfunc.unrestr, ident.unrestr] = ...
identification(mdlest, ident_restr1, shock_names, agnostic, max_trials);
[Rfunc.restr, ident.restr] = ...
identification(mdlest_restr, ident_restr1, shock_names, agnostic, max_trials);
With restrictions given as a cell array, each shock name is matched to
the restrictions that mention it, not to its position in the list:
the names can be listed in any order ({'demand','supply','mp'} as
well as {'demand','mp','supply'}), and each name labels the shock
its restrictions define, in irf, variance_decomposition,
historical_decomposition and structural_shocks alike, with or
without regime switching. Outputs that list several shocks list them
in the order the names were given; with fewer names than variables,
the remaining shocks keep their own names.
2.4. Recursive (Choleski) identification
identify(m, 'choleski', shock_names) (and
identification(m, 'choleski', shock_names)) imposes a recursive
structure: the impact matrix is the lower Cholesky factor of the
residual covariance. The recursion follows the order in which you
declared the variables, the list given to rfvar_model (or the
units and variables given to prfvar_model), even though RISE
stores the variables alphabetically. The first declared variable
responds on impact to its own shock only, the second to the first two
shocks, and so on. Shock k belongs to the k-th declared variable
and takes the k-th name:
m = rfvar_model({'ygap','pigap','Rgap'}, 'lag_length', 2);
me = estimate(m, 'data', db);
mi = identify(me, 'choleski', {'demand','supply','mp'});
r = irf(mi, 'irf_periods', 20); % r.mp: the shock of Rgap, which
% moves neither ygap nor pigap on impact
Without names, the shocks are called <variable>_SHOCK
(ygap_SHOCK, pigap_SHOCK, Rgap_SHOCK); with fewer names than
variables, the first names label the shocks of the first declared
variables and the others keep their default <variable>_SHOCK
names. A reduced-form VAR that is not identified at all uses the same
recursive factor: its shock e<k> is the shock of the k-th
declared variable.
To use another recursion order, pass 'ordering': the names of all
the variables in the order you want, or a permutation of
1:nvars that indexes the declared variables. The shock names then
follow that order:
mi = identify(me, 'choleski', {'mp','demand','supply'}, ...
'ordering', {'Rgap','ygap','pigap'});
% same as declaring rfvar_model({'Rgap','ygap','pigap'}, ...)
% and identifying with identify(me, 'choleski', {'mp','demand','supply'})
The declaration order is read from the model’s own equations (the
equation of the k-th declared variable is the k-th one, <var>{t}
= ... + e<k>{t}), so it is kept when the model is saved and loaded,
re-estimated, bootstrapped, or rebuilt by set.
This is the reduced-form counterpart of a recursive
structural VAR: an svar_model
whose \(A_{0}\) is lower triangular in the declaration order (the
zero restrictions a0(<earlier>,<later>)=0) has the same impulse
responses, its shock epsilon<k> matching the k-th Choleski shock.
Note
Change of behavior. Earlier versions ran the recursion in RISE’s
alphabetical order of the variables, so custom names given in the
order of the declared variables labeled the wrong shocks (with the
variables {'ygap','pigap','Rgap'}, the first name labeled the
shock of Rgap), and the default <variable>_SHOCK responses
followed the alphabetical recursion. Results of 'choleski' (and
of an unidentified reduced-form VAR) change unless the variables were
declared in alphabetical order. To reproduce the old results, pass
'ordering', sort(variables) and give the names in that order.
2.5. Structural shocks, IRFs and decompositions
Structural shocks:
params = [];
sshocks = struct();
sshocks.unrestr = structural_shocks(mdlest, params, Rfunc.unrestr, shock_names);
sshocks.restr = structural_shocks(mdlest_restr, params, Rfunc.restr, shock_names);
Cholesky IRFs (no identification scheme):
cholShocks = [];
myirfs = irf([mdlest, mdlest_restr], cholShocks, 40);
IRFs under the identification scheme:
params = [];
myirfs = irf([mdlest, mdlest_restr], shock_names, 40, params, Rfunc.unrestr);
Variance decomposition:
vd = variance_decomposition(mdlest_restr, params, Rfunc.restr);
Historical decomposition:
hd = historical_decomposition(mdlest_restr, params, Rfunc.restr);
The decomposition of a VAR is exact (the states are observed). Each
variable gets a ts with one column per contributor, in this order:
init: the initial conditions, i.e. the firstlag_lengthobservations propagated through the autoregressive coefficients;const: the constant term, accumulated from a zero start (when the VAR has'constant_term', true);one column per deterministic variable, under its name (the
deterministic_vars), accumulated from a zero start;one column per structural shock, under its name, accumulated from a zero start: the shock columns contain the contributions of the structural shocks only.
The columns add up to the data at every date: sum(double(d), 2),
with d = hd.(endog{1}), is the data of endog{1} from the first
date after the initial lags.
The constant and the initial conditions together make the deterministic
path of the VAR (what the data would be without any shock); plot the
shock columns alone (d(:, shock_names)) to see the story the shocks
tell. The same columns come out for svar_model, proxy_svar_model
and prfvar_model. With regime switching, the contributions of the
regimes are weighted in each period by the smoothed regime probabilities
of the filter, as for a switching DSGE model, and the output has one
page (not one page per regime). Grouping shocks (the groups
argument) or data with missing values send the computation to the
filter-based route of DSGE models, whose columns are the shocks (or
groups) and init.
2.6. Bootstrap and distributions
n = 1000;
params = bootstrap(mdlest_restr, n);
Variance-decomposition distribution:
ci = [30, 50, 68, 90];
vd = variance_decomposition(mdlest_restr, params, Rfunc.restr);
% fanchart over the draws via fanchart() + plot_fanchart()
Historical-decomposition distribution:
hd = historical_decomposition(mdlest_restr, params, Rfunc.restr);
IRF distribution:
myirfs = irf(mdlest_restr, shock_names, 40, params, Rfunc.restr);
In each case the output is multi-page (one page per bootstrap
draw); pass it to fanchart for central-tendency-plus-bands
plots (see Forecasting and simulation).
2.7. Bayesian estimation
VAR-coefficient priors are built with the var_priors
factories. The canonical informative prior is the Minnesota
family:
nlags = 4;
const = true;
exog = {};
var_prior = var_priors.minnesota(endog, const, exog, nlags, ...
tightness = 0.1, ...
lag_decay = 1.0, ...
ar_first_lag = 0.9);
The VAR prior is passed through estim_var_prior; structural /
non-VAR priors (if any) are passed through estim_priors. The
two options are independent and flat.
For a structural VAR these reduced-form priors reweight the prior of
\(A_{0}\); structural lag priors (generate_structural_prior,
Baumeister-Hamilton and Sims-Zha) do not: see
Structural VAR Modeling.
Then the unrestricted Bayesian estimate:
ve = estimate(mdl, ...
data = db, ...
estim_start_date = date2serial(db.GROWTH.start), ...
estim_end_date = date2serial(db.GROWTH.finish), ...
estim_var_prior = var_prior);
and the restricted Bayesian estimate:
ve_lr = estimate(mdl, ...
data = db, ...
estim_start_date = date2serial(db.GROWTH.start), ...
estim_end_date = date2serial(db.GROWTH.finish), ...
estim_var_prior = var_prior, ...
estim_linear_restrictions = linres);
How the mode is found. With constant parameters and one of the
conjugate templates (minnesota, inw, niw,
generate_niw_prior, generate_sims_zha_prior), estimate
does not search: the posterior mode of the likelihood (conditional on
the first lag_length observations) times the prior comes in closed
form, or by alternating two closed forms, the coefficients given the
residual covariance and the covariance given the coefficients:
priors whose coefficient covariance is \(\Sigma \otimes \Omega\) (
niw,generate_niw_prior, Sims-Zha) give the mode in one step; with Sims-Zha it is least squares on the data with the dummy observations appended;priors with a fixed coefficient covariance (
minnesota,inw) alternate the two steps until they agree to about 1e-12;linear restrictions (
estim_linear_restrictions, and the panel estimators of a panel VAR) are imposed on the coefficients in the coefficient step.
estimate prints Posterior mode in closed form, and the mode is
the one the numerical search would converge to (the search used to stop
short of it on larger systems, after minutes). It searches numerically,
as before, with regime switching, priors in estim_priors,
endogenous priors, nonlinear or general restrictions, restrictions on
the covariance, estim_mle, missing observations, data_demean or
a user-defined template, and when the closed-form mode is rejected by
the posterior kernel (an explosive VAR under solve_check_stability,
say). estim_var_closed_form = false asks for the search anyway.
A closed-form mode comes with the exact Hessian of the log posterior:
with estim_hessian_type = 'fd' (the default) or 'optimizer' it is
used instead of finite differences (hessian_source starts with
analytic), and the Laplace marginal data density is reported. On a
VAR with a hundred parameters the finite differences take minutes;
hessian(mdlest) still computes them.
params = struct();
params.ve = ve.estim_.sampler(1000);
params.ve_lr = ve_lr.estim_.sampler(1000);
2.8. Bayesian forecasting
myfkst = struct();
date_start = date2serial('2003Q1');
myfkst.ve = forecast(ve, db, date_start, params.ve);
myfkst.ve_lr = forecast(ve_lr, db, date_start, params.ve_lr);
The outputs are multi-page ts – one page per posterior
draw. fanchart + plot_fanchart give the central tendency
plus probability bands at ci = [30 50 68 90] percent.
2.9. Conditional forecasting
Condition on a path for the policy rate R:
date_start = date2serial('2003Q1');
nsteps = 12;
shock_uncertainty = false;
Rfunc = []; % no need for an identification scheme
conditions = struct();
conditions.R = {'2003Q1','2004Q4'};
myfkst.ve = forecast(ve, db, date_start, params.ve, ...
nsteps, shock_uncertainty, Rfunc, conditions);
conditions carries the conditioning ranges. The full surface
of conditional-forecasting options is in
Forecasting and simulation.
2.10. Markov-switching reduced-form VARs
A chain is declared with markov_chains and told which
parameters it governs. The shorthands expand to the reduced-form
parameter families:
Shorthand |
Expands to |
|---|---|
|
every |
|
every constant |
|
the diagonal of the residual covariance, |
|
the FULL residual covariance – variances
|
|
just the variances / just the covariances |
A reduced-form VAR carries one residual per equation, correlated
across equations. 's' therefore switches only the variances;
use 'covar' when the correlation structure should switch too:
markov_chains = struct( ...
'name' , 'vol', ...
'number_of_states' , 2, ...
'controlled_parameters' , {{'b','c','covar'}}, ...
'endogenous_probabilities', {{}}, ...
'probability_parameters', {{}});
mdl = rfvar_model(endog, lag_length = nlags, ...
constant_term = const, ...
markov_chains = markov_chains);
An entry that names nothing is an error rather than a silent
no-op, so a typo in controlled_parameters is reported at
construction.
Note
A covar_* you never set keeps its identity entry – unit
variance, zero covariance. Nobody should have to calibrate a
covariance matrix by hand just to get a VAR to solve, and
estimate seeds these from an OLS pre-fit in any case. Set
them only when you mean something other than uncorrelated unit
shocks.
A value you do supply must be a valid covariance: a negative
variance or an inconsistent pair stops the solve with return
code RFVAR_COV_NOT_PD. Unspecified defaults; specified-and-
wrong is reported.
The same machinery serves panel VARs. prfvar_model expands the
endogenous list to <unit>_<var> and builds an RFVAR on it, so
its covariance parameters are covar_sigma_<unit>_<var> and
covar_<unit_i>_<var_i>_<unit_j>_<var_j> and behave exactly as
above.
2.10.1. Priors under regime switching
An informative VAR prior is defined on the reduced form \((B, \Sigma_u)\). Under switching, each regime has its own reduced form, so the prior is evaluated at every regime and the contributions are summed – the usual reading of “the same prior applied in each regime”. Nothing extra is needed: pass the template as before and it spans the regimes:
var_prior = var_priors.minnesota(endog, const, exog, nlags, ...
tightness = 0.25);
ve = estimate(mdl, ...
data = db, ...
estim_var_prior = var_prior);
Two consequences are worth keeping in mind.
If two regimes carry an identical reduced form – because the switching is confined to the covariance, say – the density on that shared block is counted once per regime. When that duplication is unwanted, put the prior on the parameters directly through
estim_priorsinstead.estim_var_priorandestim_priorsare independent and additive. The first is the VAR-shaped template; the second carries priors on individually named parameters (including transition-probability parameters and any parameter appearing in anendogenous_probabilitiesexpression).
2.10.2. Threshold VARs: threshold_chains
A threshold VAR is a switching VAR whose transition probabilities are a
hard function of a lagged variable. You can write that chain by hand (see
the next section), but threshold_chains says it in one place and
encodes the four things the hand-written form gets wrong silently:
m = rfvar_model(endo, ...
lag_length = 1, ...
constant_term = true, ...
threshold_chains = struct( ...
'name' , 'sc', ...
'controlled_parameters', {{'b','c','covar'}}, ...
'threshold_variable' , 'ACR'));
That is all. It desugars to a plain two-state markov_chains entry with
sc_tp_1_2 = hard_threshold(ACR, sc_bar)
sc_tp_2_1 = 1-hard_threshold(ACR, sc_bar)
and one new estimable parameter, sc_bar, which is the threshold (pass
a number as threshold_value to fix it instead — see below). The model
then reports the chain it actually has: the sugar leaves no trace, so
describe_regimes, filter and everything else behave exactly as
they would for a hand-written chain.
Field |
Meaning |
|---|---|
|
chain name; creates |
|
as for |
|
an endogenous or deterministic name |
|
|
|
|
|
|
|
the |
|
2 (only two states are supported) |
State 1 is below the threshold, state 2 above.
Fixed values are literals, not parameters
threshold_value and slope each take 'estimate' or a number,
and the two do different things:
'estimate'creates a parameter —<name>_bar,<name>_gam— listed inprobability_parametersand free to be given a prior.a number is written into the transition expression as a literal. No parameter is created:
sc_tp_1_2 = hard_threshold(ACR, 0.25) sc_tp_2_1 = 1-hard_threshold(ACR, 0.25)
The asymmetry is deliberate. A fixed value that arrived as a parameter
would be one more name you must remember to set, and forgetting is
silent: the model builds, filters and returns a likelihood computed at
whatever the uncalibrated value happened to be. That is the class of
failure this sugar exists to remove, so a fixed threshold does not get a
name.
What the sugar encodes
Four things, each of which is invisible when wrong:
Which argument is the threshold. In
logistic(x,g,c)it is the third; nothing tells you that.The two probabilities must be complements. A threshold rule is memoryless in \(s_{t-1}\) — where you go depends only on \(z_{t-d}\), not on where you were. Two independent logistics give a smooth-transition Markov chain, which is a different model, and nothing warns you.
The timing.
delay = 1is written with no lag operator, because the transition matrix built at \(t\) governs the regime at \(t+1\) (see the note below). WritingACR{-1}yourself would be off by one period.The steepness.
transition = 'hard'useshard_threshold, a genuine indicator, so there is no slope to choose and no risk of leaving observations in the interior of a logistic — where they become genuine regime mixtures rather than threshold regimes.
For a smooth transition instead, transition = 'logistic' adds a slope
parameter <name>_gam and you are back to a smooth-transition VAR
(Auerbach–Gorodnichenko style) with the same declaration.
Note
hard_threshold(x,c) is the \(g \to \infty\) limit of
logistic(x,g,c), floored into \([\delta, 1-\delta]\) with
\(\delta = 10^{-10}\). It is named hard_threshold rather than
threshold because MATLAB’s Econometrics Toolbox ships a
@threshold class, and a class folder takes precedence over a
function on the path.
2.10.3. Time-varying transition probabilities
Transition probabilities may be written as expressions of the
model’s own variables, via endogenous_probabilities:
markov_chains = struct( ...
'name' , 'sc', ...
'number_of_states' , 2, ...
'controlled_parameters' , {{'b','c','covar'}}, ...
'endogenous_probabilities', {{ ...
'sc_tp_1_2 = logistic(ACR,gam,thr)', ...
'sc_tp_2_1 = 1-logistic(ACR,gam,thr)'}}, ...
'probability_parameters', {{'gam','thr'}});
Important
Timing. The transition matrix built from time-\(t\)
data governs the regime at \(t+1\). An expression written
on the contemporaneous variable ACR therefore implements
the lagged rule
Writing ACR{-1} would be off by one period. The same
convention holds in simulation, so filter and simulator agree.
With gam large the logistic approaches an indicator and the
chain becomes a threshold VAR: the regime is a deterministic
function of the data and the Hamilton filter collapses onto that
path. Keep gam steep enough that no observation lands in the
interior of the logistic – an observation that does is treated,
correctly, as a genuine mixture of the two regimes, which is no
longer a threshold model.