%Written by Ahmed Jamal Pirzada, Department of Economics, University of Bristol (aj.pirzada@bristol.ac.uk)
%This code generates forecasts, conditional forecasts and ZLB constrained forecasts for any given .mod file. 
%The results produced for non-ZLB forecast are identical to those produced by Dynare. 
%The variables and their corresponding values are to be manually set in the conditional forecast section!
%The model variables to Plot must also be manually set in the PLOT section.

% clear; 
% close all; 
% clc 

%run (set =1 for estimation or =0 for using existing results)
run = 1;

%change the model name according to the .mod file
if run==1
dynare NGSusmodelffnews_4wp3.mod console
elseif run==0
load NGSusmodelffnews_4wp3_results.mat
end

%Set constraint=1 if you want Forecast under ZLB; set constraint=2 if you
%want conditional forecast as in the forecast paper
constraint = 0;
controlvar = 'robs'; %monetary policy variable 
controlshock = 'em'; %shock associated with the monetary policy (or interest rate) rule

conditionalvar1 = 'robs'; %first variable to condition on
var1shock = 'em'; %corresponding shock for var1
var1value = 0.13; %choose a value for first conditional variable
conditionalvar2 = 'spread'; %second variable to condition on
var2shock = 'efw'; %corresponding shock for var2
var2value = 1.4; %choose a value for second conditional variable
%-----------------------------------------------------------------------



%SETTING UP SIMULATION FOR THE FORECAST PERIOD
%number of periods over which to plot forecast
iperiod=20; 
%setting the initial period for simulation which is the last period of estimation
last= size(oo_.SmoothedVariables.Mean.('b'),1);

%set the last estimated value of shocks as their starting value for simulations
ex0=[];
    for shock_iter=1:M_.exo_nbr
    ex0=[ex0; 
    oo_.SmoothedShocks.Mean.(deblank(M_.exo_names(shock_iter,:)))(last)];
    end 
    
%set matrix of shocks over the forecast period
ex=zeros(iperiod,M_.exo_nbr);
%ex(2,:)=ex0;


%set the last estimated value of variables as their starting value for simulations
y0=[];

if M_.orig_endo_nbr==M_.endo_nbr

    for endo_iter=1:M_.endo_nbr
    y0 = [y0;
    oo_.SmoothedVariables.Mean.(deblank(M_.endo_names(endo_iter,:)))(last)];
    end;

else
    
    for endo_iter=1:M_.orig_endo_nbr
    y0 = [y0;
    oo_.SmoothedVariables.Mean.(deblank(M_.endo_names(endo_iter,:)))(last)];
    end;

    %deals with the auxilary variables (if there are any in the model)
    aux = M_.endo_nbr - M_.orig_endo_nbr;
    aux = zeros(aux,1);
    y0 = vertcat(y0, aux);
   
end

%make sure decision rules were updated
[oo_.dr,info,M_,options_] = resol(0,M_,options_,oo_);

dr = oo_.dr;
iorder=1;

y_=simult_(y0,dr,ex,iorder);
% %---------------------------------------------------------------------



%SIMULATION UNDER ZERO LOWER BOUND
%Turns off monetary policy response by finding a series of monetary policy shocks. 
%Decrease the threshold (i.e. 1.0e-9) even further for even smaller changes in interest rate. Will also need to add more zeros to 0.000000001 for this to work. Computation time increases.
if constraint == 1
    disp('STARTING CONSTRAINED SIMULATION')
    y_zlb = y_;
    z = y_zlb(strmatch(controlvar,M_.endo_names,'exact'),:); %Forecast for variable you want as zero. Here the varialbe is 'r'
    for j=1:iperiod
        if z(j) <= -1.0e-8
            while z(j) <= -1.0e-8
                ex(j,strmatch(controlshock,M_.exo_names,'exact')) = ex(j,strmatch(controlshock,M_.exo_names,'exact'))+0.00000001;      
                y_zlb = simult_(y0,oo_.dr,ex,1);
                z = y_zlb(strmatch(controlvar,M_.endo_names,'exact'),:); %IRFs for variable you want as zero. Here the varialbe is 'r'
            end 
         elseif z(j) > -1.0e-8
             y_zlb = simult_(y0,oo_.dr,ex,1);
         end
    end
    disp('END OF CONSTRAINED SIMULATION')
end
% %---------------------------------------------------------------------



%SIMULATION SIMILAR TO CONDITIONAL FORECAST
y_cond = y_;
ex_cond = ex;

if constraint == 2
    disp('STARTING SIMULATION FOR ROBS')
    z_robs = y_cond(strmatch(conditionalvar1,M_.endo_names,'exact'),:);
    if z_robs(2) >= var1value
        while z_robs(2) >= var1value
            ex_cond(1,strmatch(var1shock,M_.exo_names,'exact')) = ex_cond(1,strmatch(var1shock,M_.exo_names,'exact'))-0.001;  %you may have to change the sigh before '0.001' depending on the shock and variable you are conditioning
            y_cond = simult_(y0,oo_.dr,ex_cond,1);
            z_robs = y_cond(strmatch(conditionalvar1,M_.endo_names,'exact'),:);
        end
%     elseif z_robs(1) <= .13
%         while z_robs(1) <= 0.13  
%             ex_cond(1,strmatch('em',M_.exo_names,'exact')) = ex_cond(1,strmatch('em',M_.exo_names,'exact'))+0.001;
%             y_cond = simult_(y0,oo_.dr,ex_cond,1);
%             z_robs = y_cond(strmatch('robs',M_.endo_names,'exact'),:);
%         end
    end
    disp('END OF FIRST SIMULATION')

    disp('STARTING SIMULATION FOR SPREAD')
    z_spread = y_cond(strmatch(conditionalvar2,M_.endo_names,'exact'),:);
    if z_spread(2) <= var2value
        while z_spread(2) <= var2value
            ex_cond(1,strmatch(var2shock,M_.exo_names,'exact')) = ex_cond(1,strmatch(var2shock,M_.exo_names,'exact'))+0.001;  %you may have to change the sigh before '0.001' depending on the shock and variable you are conditioning
            y_cond = simult_(y0,oo_.dr,ex_cond,1);
            z_spread = y_cond(strmatch(conditionalvar2,M_.endo_names,'exact'),:);
        end
%     elseif z_spread(1) <= 1.4
%         while z_spread(1) <= 1.4  
%             ex_cond(1,strmatch('efw',M_.exo_names,'exact')) = ex_cond(1,strmatch('efw',M_.exo_names,'exact'))+0.001;
%             y_cond = simult_(y0,oo_.dr,ex_cond,1);
%             z_spread = y_cond(strmatch('spread',M_.endo_names,'exact'),:);
%         end
    end
    disp('END OF SECOND SIMULATION')
    
    disp('IMPOSING ZLB')
    y_cond_zlb = y_cond;
    z_zlb = y_cond_zlb(strmatch(controlvar,M_.endo_names,'exact'),:); %Forecast for variable you want as zero. Here the varialbe is 'r'
    for j=1:iperiod
        if z_zlb(j+1) <= -1.0e-8
            while z_zlb(j+1) <= -1.0e-8
                ex_cond(j,strmatch(controlshock,M_.exo_names,'exact')) = ex_cond(j,strmatch(controlshock,M_.exo_names,'exact'))+0.001;      
                y_cond_zlb = simult_(y0,oo_.dr,ex_cond,1);
                z_zlb = y_cond_zlb(strmatch(controlvar,M_.endo_names,'exact'),:); %IRFs for variable you want as zero. Here the varialbe is 'r'
            end 
%          elseif z_zlb(j) > -1.0e-8
%              y_cond_zlb = simult_(y0,oo_.dr,ex,1);
         end
    end
    disp('END OF ZLB')
end
% %---------------------------------------------------------------------



%PLOT

plot_vars={'dy';'dc';'dinve';'pinfobs';'robs';'spread'}; %choose the variables you want to plot the IRFs for
plot_vars_heading={'OutputGr';'ConsGr';'InvGr';'Infl';'Fed Rate';'Spread'}; %names of chosen variables

IRF_matrix=[];

if constraint == 0
    figure
    for i=1:6
        subplot(3,2,i)

        %Creates a matrix with IRFs for variables of interest.
        x = y_(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:); %selecting the IRF for the i'th variable in plot_vars
        IRF_matrix= [IRF_matrix x']; %creates a matrix with IRFs of variables in plot_vars. For copy-pasting to excel

        %Creates the subplots.
        plot_=y_(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:);
        plot(1:iperiod,plot_(1,1:iperiod),'r','LineWidth',1.2);
        title(plot_vars_heading(i,:))
        axis tight
        hold on
        plot((1:iperiod),zeros(iperiod,1),'k')
    end
elseif constraint == 1
    figure
    for i=1:6
        subplot(3,2,i)

        %Creates a matrix with IRFs for variables of interest.
        x_zlb = y_zlb(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:); %selecting the IRF for the i'th variable in plot_vars
        IRF_matrix_zlb= [IRF_matrix x_zlb']; %creates a matrix with IRFs of variables in plot_vars. For copy-pasting to excel

        %Creates the subplots.
        plot_=y_(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:);
        plot_zlb=y_zlb(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:);
        plot(1:iperiod,plot_(1,1:iperiod),'r',1:iperiod,plot_zlb(1,1:iperiod),'--b','LineWidth',1.2);
        if i==5
        legend('Baseline','No Policy','location','northwest')
        legend boxoff
        end
        title(plot_vars_heading(i,:))
        axis tight
        hold on
%        plot((1:iperiod),zeros(iperiod,1),'k')
    end
elseif constraint == 2
        figure
    for i=1:6
        subplot(3,2,i)

        %Creates a matrix with IRFs for variables of interest.
        x_cond = y_cond(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:); %selecting the IRF for the i'th variable in plot_vars
        IRF_matrix_zlb= [IRF_matrix x_cond']; %creates a matrix with IRFs of variables in plot_vars. For copy-pasting to excel

        %Creates the subplots.
        plot_=y_(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:);
        plot_cond=y_cond(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:);
        plot_cond_zlb=y_cond_zlb(strmatch(plot_vars(i,:),M_.endo_names,'exact'),:);
        plot(1:iperiod,plot_(1,1:iperiod),'r',1:iperiod,plot_cond(1,1:iperiod),'--b',1:iperiod,plot_cond_zlb(1,1:iperiod),'-.k','LineWidth',1.2);
        if i==5
        legend('Baseline','Conditional','With ZLB','location','northwest')
        legend boxoff
        end
        title(plot_vars_heading(i,:))
        axis tight
        hold on
%        plot((1:iperiod),zeros(iperiod,1),'k')
    end
end