clear all;close all
clc
addpath('occbin_21Nov/toolkit_files')


nperiods=100;
maxiter=1500;

modnam = 'model_childwork_noZLB';
modnamstar = 'model_childwork_ZLB'; 

constraint = 'Rs< -log(Rss)';
constraint_relax = 'Rtr> -log(Rss)';

param_save(0,0,0,3,0.25); % exogenous and constant QE


variables={'Y','X','Lamp','Lamb','LL','L','Lb','w','Cp','Cb','Cineq','Pi','Rs','rs','Rb','rb' ...
    'Rtr','Rre','Q','Omega','bFI','bCB','b','XFI','re','s','Theta','QE','pm'...      
    'A','pstar','x1','x2','vp','exr','lnY','lnCp','lnCb','lnCineq','lnL' ...
    'lnRs','lnrs','lnRb','lnrb','lnPi','lnQ','div','div_FI','surplus','transferb','laborinc'};



variables_IRF= {'lnY' 'Output' 1
                'lnPi' 'Inflation' 4
                'exr'   'Yield Spread' 4
                'lnRb'  'Nominal long yield' 4
                'lnrb'  'Real long yield' 4
                'lnRs'  'Nominal Deposit rate' 4
                'lnrs'  'Real Deposit rate' 4
                're'    'CB balance sheet' 1/100/1/4
                'lnCp'  'Consumption (savers)' 1
                'lnCb'  'Consumption (borrowers)' 1
                'lnCineq' 'Consumption inequality' 1
                'L' 'Labor supplied (savers)' 1
                'Lb' 'Labor supplied (borrowers)' 1
                'w'     'Real wage' 1
                'div'   'Real profits' 1
                };
% find the order of variables
for i=1:length(variables_IRF)
    eval(['p_' variables_IRF{i,1} '=find(strcmp(variables, variables_IRF{i,1} ));'])
end




%%
% case 1: no ZLB, MP, child work
% case 1: no ZLB, exogenous QE, child work
% case 2: ZLB, exogenous QE, child work
% case 4: ZLB, enodgenous QE, child work


% case 5: ZLB, exogenous QE, fixed transfer

baseline1= zeros(7,2);

scenario= baseline1;
baseline=baseline1;

baseline(1:6,1)= -25;
scenario(end,2)=-1;
    


%----------------- No ZLB -------------------------------
% CASE 1
% Credit shock for ZLB binding and MP shock
irfshock = char('et','er');

% First time we solve simulation only with baseline shocks
[sim_MP] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline,irfshock,nperiods,maxiter);

% Second time we solve simulation with baseline shocks and scenario
[sim_MP2] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline+scenario,irfshock,nperiods,maxiter);
irf1 = sim_MP2 - sim_MP; % no ZLB, no QE


% CASE 2
% Credit shock for ZLB binding and exogenous QE shock
param_save(0.6,0,0,3,0.25);

% First time we solve simulation only with baseline shocks
irfshock = char('et','eq');                    
shock_size = -4.5;

[sim_QE] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline,irfshock,nperiods,maxiter);

% Second time we solve simulation with baseline shocks and scenario
[sim_QE2] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline+scenario*shock_size,irfshock,nperiods,maxiter);
irf2 = sim_QE2 - sim_QE; % no ZLB, QE



%%
%********************* Impose ZLB **********************************
% CASE 3
irfshock = char('et','eq');
shock_size = -1.5;
% First time we solve simulation only with baseline shocks
[~, sim_zlb_exo] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline,irfshock,nperiods,maxiter);

% Second time we solve simulation with baseline shocks and scenario
[~, sim_zlb2_exo] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline+scenario*shock_size,irfshock,nperiods,maxiter);

irf3 = sim_zlb2_exo - sim_zlb_exo;

% CASE 4
% Endogenous QE
param_save(rhoR,0.25,5,3,0.25);


% First time we solve simulation only with baseline shocks
irfshock = char('et','eq');
shock_size = -4.0;

[~, sim_zlb_endo] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline,irfshock,nperiods,maxiter);

% Second time we solve simulation with baseline shocks and scenario
[~, sim_zlb2_endo] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline+scenario*shock_size,irfshock,nperiods,maxiter);

irf4 = sim_zlb2_endo - sim_zlb_endo;                    
                    
                    
%% PLOT
set(0, 'DefaultAxesFontSize',14);
set(0, 'DefaultTextFontSize', 14);
set(0, 'DefaultUicontrolFontSize',14);
co = [0    0.4470    0.7410
    %0.8500    0.3250    0.0980
    0.9290    0.6940    0.1250
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];
set(groot,'defaultAxesColorOrder',co)

% Display options
OptionGreycolor = 0;    % if want grey colorscheme instead of standard colors

vLinestyle = {'-','-.','--',':'};
FontsizeDefault = 12;
FontsizeAxis = 12;
FontSizeLegend = 12;
Fonttype = 'times';
LinewidthDefault = 1.6;
LinewidthAlt = 1;
ColorZeros = 'k';
StyleZeros = '-';

TitleSize = 14;
LabelSize = 12;
Font = 'times';

if OptionGreycolor == 0
vColors = {[4,30,150]/255,[1 0.5 0],[152,58,68]/255,[30,144,250]/255};
else
% vColors = repmat(linspace(0.2,0.6,4).',1,4);
% vColors = mat2cell(vColors,ones(1,size(vColors,1)),4);
vColors = {[0.2,0.2,0.2],[0.46,0.46,0.46],[0.6,0.6,0.6]};
end


for i = 1:length(variables_IRF)
    subplot(4,4,i)
    
    IRF_variable = eval(['sim_MP2(:,p_' variables_IRF{i,1} ') - sim_MP(:,p_' variables_IRF{i,1} ');']);
    IRF_variable = IRF_variable * 100*variables_IRF{i,3};
    plot(1:nperiods,IRF_variable,'Linewidth',1.5), hold on
    
    
    IRF_variable = eval(['sim_QE2(:,p_' variables_IRF{i,1} ') - sim_QE(:,p_' variables_IRF{i,1} ');']);
    IRF_variable = IRF_variable * 100*variables_IRF{i,3};
    plot(1:nperiods,IRF_variable,'Linewidth',1.5)
    
    IRF_variable = eval(['sim_zlb2_exo(:,p_' variables_IRF{i,1} ') - sim_zlb_exo(:,p_' variables_IRF{i,1} ');']);
    IRF_variable = IRF_variable * 100*variables_IRF{i,3};
    plot(1:nperiods,IRF_variable,'-.','Linewidth',1.5)
                   
    IRF_variable = eval(['sim_zlb2_endo(:,p_' variables_IRF{i,1} ') - sim_zlb_endo(:,p_' variables_IRF{i,1} ');']); 
    IRF_variable = IRF_variable * 100*variables_IRF{i,3};
    plot(1:nperiods,IRF_variable,'Linewidth',1.5)
    
    grid on
    box on                    
    title(variables_IRF{i,2})
    xlim([1 15])
 
end

hL=legend('no ZLB, MP','no ZLB, QE','ZLB,QE', 'ZLB, endogenous QE');%,'Location', 'SouthOutside','FontSize','FontSizeLegend') %,'ZLB (Endogenous QE)')
set(hL,'interpreter','latex','fontname','times','Location','southoutside','FontSize',FontSizeLegend)
newPosition = [0.3765    0.9542    0.2811    0.0373];
newUnits = 'normalized';
set(hL,'Position', newPosition,'Units', newUnits,'Orientation','horizontal');

