2. Reduced-form VAR Modeling

2.1. Description

\[y_{t} = C(r_{t}) x_{t} + B_{1}(r_{t}) y_{t-1} + \cdots + B_{p}(r_{t}) y_{t-p} + u_{t}\]

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);

2.4. 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);

2.5. 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.6. 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_priors; structural / non-VAR priors (if any) are passed through estim_priors. The two options are independent and flat.

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_priors = 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_priors          = var_prior, ...
    estim_linear_restrictions = linres);

2.7. Posterior sampling

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.