4. Reduced-form VAR Modeling

4.1. Description

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

With \(r_{t}=1,2,...,h\) and transition probabilities \(p_{r_{t},r_{t+1}}\left( I_{t}\right)\)

4.2. Quick-start examples

4.2.1. A constant-parameter VAR

Collecting and transforming data

% CLVMNACSCAB1GQNO : GDP Norway
% IR3TIB01NOQ156N : 3-month Interbank interest rate
% NORCPGRLE01IXOBQ : CPI excluding food and energy
% CCUSSP01NOQ650N : Spot Exchange rate : 1NOK = x USD
% POILWTIUSDQ : Global price of WTI Crude
xrange='1990Q1:2022Q3';
rawdb=fetch_fred({'NORCPGRLE01IXOBQ','IR3TIB01NOQ156N','CLVMNACSCAB1GQNO',...
    'CCUSSP01NOQ650N','POILWTIUSDQ'});
rawdb1.P=rawdb(1).series(xrange);
rawdb1.INTRATE=rawdb(2).series(xrange);
rawdb1.Y=rawdb(3).series(xrange);
rawdb1.EXRATE=rawdb(4).series(xrange);
rawdb1.POIL=rawdb(5).series(xrange);
rawdb=rawdb1; clear rawdb1

db=struct();
db.PAI=rawdb.P/lag(rawdb.P,1);
db.R=1+rawdb.INTRATE/100;
db.GROWTH=rawdb.Y/lag(rawdb.Y,1);
db.EXRATE=1/rawdb.EXRATE; % so that increase = depreciation of NOK
db.PAIOIL=rawdb.POIL/lag(rawdb.POIL,1);

Setting up the Reduced.form VAR

endog={'PAIOIL','GROWTH','PAI','R','EXRATE'};

exog={};

nlags=4;

const=true;

mdl = rfvar(endog, exog, nlags, const);

Estimating the VAR using classical techniques

data_range={db.GROWTH.start,db.GROWTH.finish};

mdlest=estimate(mdl,db,data_range);

Restrictions on the VAR : Domestic variables do not affect oil prices

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); %#ok<SAGROW>
%     end
% end

Estimate the restricted VAR

prior=[];

mdlest_restr=estimate(mdl,db,data_range,prior,linres);

Identification

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','+'
    % first set
    'PAIOIL{0}@demand',0
    'PAIOIL{0}@costpush',0
    'PAIOIL{0}@mp',0
    'PAIOIL{0}@forex',0
    % second set
    'GROWTH{0}@costpush',0
    'GROWTH{0}@mp',0
    'GROWTH{0}@forex',0
    % third set
    'PAI{0}@mp',0
    'PAI{0}@forex',0
    % fourth set
    '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);

disp(ident.restr)

disp(ident.unrestr)

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

% plots
%-------

titel='Structural shocks';

figure('name',titel);

for ii=1:numel(shock_names)

    thisname=shock_names{ii};

    subplot(3,2,ii)

    d=[sshocks.unrestr.(thisname),sshocks.restr.(thisname)];

    plot(d,'linewidth',2)

    xtickangle(45)

    title(thisname)

    if ii==1

        legend({'unrestr','restr'})

    end

end

% xrotate(45)

[~,h]=sup_label(titel,'t');

set(h,'fontsize',12)

Compare CHOLESKI irfs

cholShocks=[];

myirfs=irf([mdlest,mdlest_restr],cholShocks,40);

cholShocks=fieldnames(myirfs);

tex=struct();
tex.PAIOIL='Oil price inflation';
tex.GROWTH='GDP growth';
tex.PAI='Inflation';
tex.R = 'Interest rate';
tex.EXRATE='exchange rate';

for ishock=1:numel(cholShocks)

    shock=cholShocks{ishock};

    titel=['Impulse responses to a ',shock,' shock'];

    figure('name',titel);

    for ii=1:numel(endog)

        subplot(3,2,ii)

        plot(myirfs.(shock).(endog{ii}),'linewidth',2)

        title(tex.(endog{ii}))

        if ii==1

            legend({'unrestr','restr'})

        end

    end

    [~,h]=sup_label(titel,'t');

    set(h,'fontsize',12)

end

Compare irfs based on the identification scheme

params=[];

% Note the two models have the same identification !!!
% This is why we can run them together
myirfs=irf([mdlest,mdlest_restr],shock_names,40,params,Rfunc.unrestr);

for ishock=1:numel(shock_names)

    shock=shock_names{ishock};

    titel=['Impulse responses to a ',shock,' shock'];

    figure('name',titel);

    for ii=1:numel(endog)

        subplot(3,2,ii)

        plot(myirfs.(shock).(endog{ii}),'linewidth',2)

        title(tex.(endog{ii}))

        if ii==1

            legend({'unrestr','restr'})

        end

    end

    [~,h]=sup_label(titel,'t');

    set(h,'fontsize',12)

end

Variance decomposition

params=[];

vd=variance_decomposition(mdlest_restr,params,Rfunc.restr);

% plot decompositions
%---------------------

range='0:50'; % pick a range for the plots

figure('name','Variance Decomposition');

for iv=1:numel(endog)

    d=vd.conditional.(endog{iv});

    subplot(3,2,iv)

    plot_decomp(range,d)

    if iv==1

        legend(shock_names,'location','SE',...
            'Orientation','horizontal')

    end

    title(tex.(endog{iv}))

end

Historical decomposition

params=[];

hd=historical_decomposition(mdlest_restr,params,Rfunc.restr);

% plot decompositions
%----------------------

shock_only=true;

titel='Model with block exogeneity: Historical Decomposition';

figure('name',titel);

for iv=1:numel(endog)

    d=hd.(endog{iv});

    subplot(3,2,iv)

    if shock_only

        shock_tex=shock_names;

        d=d(shock_names);

    else

        shock_tex=d.varnames;

    end

    plot_decomp(d)

    if iv==1

        legend(shock_tex,'location','SE','Orientation','horizontal')

    end

    title(tex.(endog{iv}))

end

[~,h]=sup_label(titel,'t');

set(h,'fontsize',12)

Bootstrap

n=1000;

params=bootstrap(mdlest_restr,n);

Variance decomposition distribution

ci=[30,50,68,90];

vd=variance_decomposition(mdlest_restr,params,Rfunc.restr);

shock_tex=shock_names;

myrange='1:50';

for iv=1:numel(endog)

    vname=tex.(endog{iv});

    titel=['Variance Decomposition (in %) of ',vname];

    figure('name',titel);

    d=vd.conditional.(endog{iv});

    d.varnames=shock_names;

    d=pages2struct(d);

    contributors=fieldnames(d); % = shock_names

    for ii=1:numel(contributors)

        subplot(3,2,ii)

        % note we are multiplying by 100, this just by pure convenience
        %--------------------------------------------------------------
        out=fanchart(100*d.(contributors{ii})(myrange),ci);

        plot_fanchart(out)

        title(contributors{ii})

        axis tight

    end

    [~,h]=sup_label(titel,'t');

    set(h,'fontsize',12)

end

Historical decomposition distribution

ci=[30,50,68,90];

% compute decompositions

    hd=historical_decomposition(mdlest_restr,...
        params,Rfunc.restr);

% plot decompositions

for iv=1:numel(endog)

    vname=tex.(endog{iv});

    titel=['Model with restrictions : Historical Decomposition of ',vname];

    figure('name',titel);

    d=pages2struct(hd.(endog{iv}));

    contributors=fieldnames(d);

    for ii=1:numel(contributors)

        subplot(5,2,ii)

        out=fanchart(d.(contributors{ii}),ci);

        plot_fanchart(out)

        xtickangle(45)

        title(contributors{ii})

        axis tight

    end

    [~,h]=sup_label(titel,'t');

    set(h,'fontsize',12)

end

IRF distribution

    myirfs=irf(mdlest_restr,shock_names,...
        40,params,Rfunc.restr);

% IRFs plots
%--------------

for ishock=1:numel(shock_names)

    shock=shock_names{ishock};

    titel=['Model with restrictions : IRFs to a ',shock,' shock'];

    figure('name',titel);

    for ii=1:numel(endog)

        subplot(3,2,ii)

        d=myirfs.(shock).(endog{ii});

        out=fanchart(d,ci);

        plot_fanchart(out)

        title(tex.(endog{ii}))

        axis tight

    end

    [~,h]=sup_label(titel,'t');

    set(h,'fontsize',12)

end

Bayesian estimation

v=rfvar(endog,exog,nlags,const);

%set prior

var_prior=rfvar.prior_template();
% modify as needed

var_prior.type='sz';

prior=struct('var',var_prior);
% prior.nonvar
% prior.endogenous

% unrestricted model
%--------------------

ve=estimate(v,db,data_range,prior);%,restrictions

% restricted model
%------------------

ve_lr=estimate(v,db,data_range,prior,linres);

Posterior sampling of parameters

params=struct();
params.ve=ve.estim_.sampler(1000);
params.ve_lr=ve_lr.estim_.sampler(1000);

Bayesian forecasting

myfkst=struct();
date_start='2003Q1';
myfkst.ve=forecast(ve,db,date_start,params.ve);

myfkst.ve_lr=forecast(ve_lr,db,date_start,params.ve_lr);

% set environment

ci=[30,50,68,90];

modelnames={'ve','ve_lr'};

% Forecast plots

for jj=1:numel(modelnames)

    modname=modelnames{jj};

    figure('name',['model (',modname,') Forecasts of Norwegian Data']);

    for ii=1:numel(endog)

        subplot(3,2,ii)

        d=myfkst.(modname).(endog{ii});

        out=fanchart(d,ci);

        plot_fanchart(out,[244, 122, 66]/244)

        title(tex.(endog{ii}))

        axis tight

    end

    xrotate(45)

end

Conditional forecasting

myfkst=struct();
% date_start=[];
date_start='2003Q1';
nsteps=12;
shock_uncertainty=false;
Rfunc=[]; % No need for identification

conditions=struct();
conditions.R={'2003Q1','2004Q4'}; % range over which we want to condition

myfkst.ve=forecast(ve,db,date_start,params.ve,nsteps,...
    shock_uncertainty,Rfunc,conditions);

myfkst.ve_lr=forecast(ve_lr,db,date_start,params.ve_lr,nsteps,...
    shock_uncertainty,Rfunc,conditions);

% set environment

ci=[30,50,68,90];

% Forecast plots

for jj=1:numel(modelnames)

    modname=modelnames{jj};

    figure('name',['model (',modname,') Forecasts of Norwegian Data']);

    for ii=1:numel(endog)

        subplot(3,2,ii)

        d=myfkst.(modname).(endog{ii});

        out=fanchart(d,ci);

        plot_fanchart(out,[244, 122, 66]/255)

        title(tex.(endog{ii}))

        axis tight

    end

    xrotate(45)

end

4.2.2. Adding regime switching with constant transition probabilities

%% adding constant transition probabilities

mc=struct();
mc.name='policy';
mc.number_of_states=2;
mc.controlled_parameters={'b(R,:)'};
mc.endogenous_probabilities=[];
mc.probability_parameters=[];

mdl = rfvar(endog, exog, nlags, const,mc);

%% Set priors and Estimate the VAR

switch_prior=struct();
switch_prior.policy_tp_1_2={0.5,0.1,0.3,'beta'};
switch_prior.policy_tp_2_1={0.5,0.1,0.3,'beta'};

prior=struct();
prior.nonvar=switch_prior;

% ___Optional prior on the VAR__
% var_prior=svar.prior_template();
% var_prior.type='sz';
% prior.var=var_prior;

data_range={db.GROWTH.start,db.GROWTH.finish};

mdlest=estimate(mdl,db,data_range,prior);

4.2.3. Adding regime switching with time~varying transition probabilities

%% adding time-varying transition probabilities

mc=struct();
mc.name='policy';
mc.number_of_states=2;
mc.controlled_parameters={'b(R,:)'};
mc.endogenous_probabilities={
    'policy_tp_1_2=1/(1+exp(a12*(PAI-c12)))'
    'policy_tp_2_1=1/(1+exp(a21*(PAI-c21)))'
    };
mc.probability_parameters={'a12','a21','c12','c21'};

mdl = rfvar(endog, exog, nlags, const,mc);

%% Set priors and estimate the RFVAR

switch_prior=struct();
switch_prior.a12={0,0,1,'normal'};
switch_prior.a21={0,0,1,'normal'};
switch_prior.c12={1.05^.25,1.05^.25,0.05,'gamma'};
switch_prior.c21={1.05^.25,1.05^.25,0.05,'gamma'};

prior=struct();
prior.nonvar=switch_prior;

% ___Optional prior on the VAR__
% var_prior=svar.prior_template();
% var_prior.type='sz';
% prior.var=var_prior;

data_range={db.GROWTH.start,db.GROWTH.finish};

mdlest=estimate(mdl,db,data_range,prior);