%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Thie file replicates  results in Figure 1 of 
% Eric Sims and Jing Cynthia Wu "Evaluating Central Banks' Tool Kit: Past, 
% Present, and Future", Journal of Monetary Economics, forthcoming.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear 
close all

clc
addpath c:\dynare\4.6.4\matlab
%addpath c:\dynare\4.5.6\matlab
addpath('Occbin_update-master\toolkit_files')
%addpath('occbin_20140630\toolkit_files')
addpath('functions')
cd('models')

scale_QE = -7;

nperiods=20;
maxiter=30;

modnam = 'model';
modnamstar = 'model_zlb'; 

constraint = 'RTR<-log(Res)';
constraint_relax = 'RTR>-log(Res)';



param_save(0,0)
param_ss


variabes = {'RF' 'QF' 'RB' 'QB' 'RE' 'Rd' 'infl' 'Lam' 'Omega' 'lam' 'a' 'phi' 'mu' 'C' 'L' 'mrs' 'wr' 'f1' 'f2' 'w' 'Ld' ... 
   'pm' 'Y' 'Ym' 'u' 'K' 'pk' 'M' 'Ih' 'fpr' 'I' 'G' 'bG' 'T' 'Tcb' 'fcb' 'bcb' 're'  'vw' 'f' 'b' 'd' 'n' ...
    'RLF' 'RLB' 'QEH' 'RLEH' 'A' 'mui' 'psic' 'psil' 'theta' 'Delta' 'returnf' 'returnb' 'returnfb' ...
    'phi_1' 'phi_2' 'omega' 'spreadf' 'spreadb' 'RTR' 'rlf_real' 'rlb_real' 'M2' 'YY' 'GG' 'dY'};


variables_IRF = {'Y' 'Output' 1
    'I' 'Investment' 1
    'C' 'Consumption' 1
    'infl' 'Inflation' 4
    'RTR' 'TR rate' 4
    'RE' 'Policy rate' 4
    'Rd' 'Deposit rate' 4
    're' 'CB Balance sheet' 1/100/Ys/4
    'rlf_real' 'Real long yield' 4
    'RLF' 'Nominal long yield' 4
    'spreadf' 'Spread' 4
    };

% find the order of variables
for i = 1:length(variables_IRF)
    eval(['p_' variables_IRF{i,1} '= find(strcmp(variabes, variables_IRF{i,1} ));'])
end






baseline1= zeros(7,3);

scenario= baseline1;
scenario(end,2) = -1;

baseline=baseline1;
baseline(1:6,1) = 1.5;


%===============================MP=========================================
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);

%==================================QE======================================
irfshock = char('et','ef'); % Shocks we look at: preference and policy shock
shock_size = scale_QE;

% First time we solve simulation only with baseline shocks
[~, 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);


%===================================FG=====================================
irfshock = char('et','er');
shock_size = 2.2;
% First time we solve simulation only with baseline shocks
[~, sim_FG] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  baseline,irfshock,nperiods,maxiter);

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


%===================================NIRP===================================
irfshock = char('et','er');
shock_size = 2.4;
modnam = 'model';
modnamstar = 'model_nirp'; % for NIRP -- only deposit rate constrained

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

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

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)


for i = 1:length(variables_IRF)
    subplot(3,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_NIRP2(:,p_' variables_IRF{i,1} ') - sim_NIRP(:,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_FG2(:,p_' variables_IRF{i,1} ') - sim_FG(:,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})
end

legend('MP','QE','NIRP','FG')
cd('..')