clear all
close all
addpath c:\dynare\5.0\matlab


% set parameters
beta = 0.995;
kappa = 1-40^(-1);
spF = 0.03;
spB = 0.01;
psi = 0.81;
epsip = 11;
epsiw = 11;
alpha = 1/3;
delta0 = 0.025;
g = 0.2;
bb = 0.8; 
bcs = 0.1; % central bank balance sheet / Y
bcbGs = 0.9; % percentage of balance sheet holding government bonds
by = 0.5; % steady state government debt / GDP
levs = 5;
sigma = 0.95;

% solve for steady state
Rds = beta^(-1);
RFs = (1+spF)^(1/4)*Rds;
RBs = (1+spB)^(1/4)*Rds;
Qs = (RFs - kappa)^(-1);
QBs = (RBs - kappa)^(-1);
M1s = beta/(Qs*(1-beta*kappa));
M2s = 1 + (M1s-1)*psi;
pws = (epsip-1)/epsip;
Ks = (alpha*pws/(M2s*(1/beta - (1-delta0))))^(1/(1-alpha));
Ys = Ks^(alpha);
Is = delta0*Ks;
ws = (1-alpha)*pws*Ks^(alpha);
delta1 = (alpha*pws*Ks^(alpha-1)/M2s);
Cs = (1-g)*Ys - Is;
mus = (1/Cs)*(1-beta*bb)/(1-bb);
mrss = ((epsiw-1)/epsiw)*ws;
chi = mus*mrss;
fws = psi*Is/(Qs*(1-kappa));
res = bcs*Ys;
bcbs = bcbGs*res/QBs;
fcbs = (res - QBs*bcbs)/Qs;
fs = fws - fcbs;
bGs = by*Ys/QBs;
bs = bGs - bcbs;
ns = (Qs*fs + QBs*bs + res)/levs;
ds = Qs*fs + QBs*bs + res - ns;
Delta = (RBs - Rds)/(RFs - Rds);
phis = (Qs*fs + Delta*QBs*bs)/ns;
thetas = (1 - sigma + phis*beta*(1-sigma)*(RFs - Rds))/((1-sigma)*phis - beta*sigma*phis^2*(RFs - Rds));
X = ns - sigma*( (RFs-Rds)*Qs*fs + (RBs - Rds)*QBs*bs + Rds*ns);
Omegas = 1 - sigma + sigma*phis*thetas;
lams = (thetas/(beta*(RFs - Rds)*(1-sigma + sigma*phis*thetas)) - 1)^(-1);

% parameters unrelated to steady state
phip = 0.75;
phiw = 0.75;
psik = 2; % capital adjustment cost
delta2 = 0.01; % squared term capital adjustment cost
eta = 1; % Frisch elasticity

% Taylor rule
rhor = 0.8;
phipi = 1.5;
phiy = 0.15;
sr = 0.0025;

% Shock processes
rhof = 0.97;
rhob = 0.97;
sf = 0.01;
sb = 0.01;
rhoA = 0.9;
sA = 0.01;
rhot = 0.95;
st = 0.01;
rhoG = 0.9;
sG = 0.01;
rhoB = 0.9;
sB = 0.9;

addpath('occbin_20140630\toolkit_files')
%addpath('functions')
nperiods=29;
maxiter=50;

modnam = 'sw_2020';
modnamstar = 'sw_2020_zlb'; 

constraint = 'Rre<-0.005';
constraint_relax = 'Rtr>-0.005';

save param_sw beta kappa spF spB psi epsip epsiw alpha delta0 g bb bcs bcbGs by levs sigma Rds RFs RBs Qs QBs M1s M2s pws Ks Ys Is ws delta1 Cs mus mrss chi fws res bcbs fcbs fs bGs bs ns ds Delta phis thetas X Omegas lams phip phiw psik delta2 eta rhor phipi phiy sr rhof rhob sf sb rhoA sA rhot st rhoG sG rhoB sB

baseline1= zeros(10,3);

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

baseline=baseline1;
baseline(1:9,1) = 25.5;



% Do private QE shocks
irfshock = char('et','ef');

%%%%%%%%%%%%% No ZLB
%%%%%%%%%%%%%%%%%%
% 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);
irf0 = sim_MP2 - sim_MP; % IRF with no Taylor rule

%%%%%%%%%%%%%%%% Impose ZLB
%%%%%%%%%%%%%%%%%%%%%%

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

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



irf1 = sim_zlb2 - sim_zlb;

T = size(irf1(10:end,1),1);
t = 1:T;

% variable positions
pY = 49; % note this is log
pC = 50; 
pI = 51;
pL = 52;
pPi = 53;
pRd = 54;


figure
subplot(2,3,1)
plot(t,irf0(10:end,pY),'-k',t,irf1(10:end,pY),'--k','Linewidth',1.5)
title('Output')

subplot(2,3,2)
plot(t,irf0(10:end,pC),'-k',t,irf1(10:end,pC),'--k','Linewidth',1.5)
title('Consumption')

subplot(2,3,3)
plot(t,irf0(10:end,pI),'-k',t,irf1(10:end,pI),'--k','Linewidth',1.5)
title('Investment')

subplot(2,3,4)
plot(t,irf0(10:end,pL),'-k',t,irf1(10:end,pL),'--k','Linewidth',1.5)
title('Labor')

subplot(2,3,5)
plot(t,irf0(10:end,pPi),'-k',t,irf1(10:end,pPi),'--k','Linewidth',1.5)
title('\pi')

subplot(2,3,6)
plot(t,irf0(10:end,pRd),'-k',t,irf1(10:end,pRd),'--k','Linewidth',1.5)
title('r^{d}')




% Do government spending shocks
irfshock = char('et','eG');

%%%%%%%%%%%%% No ZLB
%%%%%%%%%%%%%%%%%%
% 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);
irf0 = sim_MP2 - sim_MP; % IRF with no Taylor rule

%%%%%%%%%%%%%%%% Impose ZLB
%%%%%%%%%%%%%%%%%%%%%%

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

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



irf1 = sim_zlb2 - sim_zlb;

T = size(irf1(10:end,1),1);
t = 1:T;

% variable positions
pY = 49; % note this is log
pC = 50; 
pI = 51;
pL = 52;
pPi = 53;
pRd = 54;


figure
subplot(2,3,1)
plot(t,irf0(10:end,pY),'-k',t,irf1(10:end,pY),'--k','Linewidth',1.5)
title('Output')

subplot(2,3,2)
plot(t,irf0(10:end,pC),'-k',t,irf1(10:end,pC),'--k','Linewidth',1.5)
title('Consumption')

subplot(2,3,3)
plot(t,irf0(10:end,pI),'-k',t,irf1(10:end,pI),'--k','Linewidth',1.5)
title('Investment')

subplot(2,3,4)
plot(t,irf0(10:end,pL),'-k',t,irf1(10:end,pL),'--k','Linewidth',1.5)
title('Labor')

subplot(2,3,5)
plot(t,irf0(10:end,pPi),'-k',t,irf1(10:end,pPi),'--k','Linewidth',1.5)
title('\pi')

subplot(2,3,6)
plot(t,irf0(10:end,pRd),'-k',t,irf1(10:end,pRd),'--k','Linewidth',1.5)
title('r^{d}')


% Do credit shocks
irfshock = char('et','et2'); % note extra shock that is the same to rig this to work

%%%%%%%%%%%%% No ZLB
%%%%%%%%%%%%%%%%%%
% 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);
irf0 = sim_MP2 - sim_MP; % IRF with no Taylor rule

%%%%%%%%%%%%%%%% Impose ZLB
%%%%%%%%%%%%%%%%%%%%%%

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

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



irf1 = sim_zlb2 - sim_zlb;

T = size(irf1(10:end,1),1);
t = 1:T;

% variable positions
pY = 49; % note this is log
pC = 50; 
pI = 51;
pL = 52;
pPi = 53;
pRd = 54;


figure
subplot(2,3,1)
plot(t,irf0(10:end,pY),'-k',t,irf1(10:end,pY),'--k','Linewidth',1.5)
title('Output')

subplot(2,3,2)
plot(t,irf0(10:end,pC),'-k',t,irf1(10:end,pC),'--k','Linewidth',1.5)
title('Consumption')

subplot(2,3,3)
plot(t,irf0(10:end,pI),'-k',t,irf1(10:end,pI),'--k','Linewidth',1.5)
title('Investment')

subplot(2,3,4)
plot(t,irf0(10:end,pL),'-k',t,irf1(10:end,pL),'--k','Linewidth',1.5)
title('Labor')

subplot(2,3,5)
plot(t,irf0(10:end,pPi),'-k',t,irf1(10:end,pPi),'--k','Linewidth',1.5)
title('\pi')

subplot(2,3,6)
plot(t,irf0(10:end,pRd),'-k',t,irf1(10:end,pRd),'--k','Linewidth',1.5)
title('r^{d}')

% Do productivity shocks
irfshock = char('et','eA');

%%%%%%%%%%%%% No ZLB
%%%%%%%%%%%%%%%%%%
% 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);
irf0 = sim_MP2 - sim_MP; % IRF with no Taylor rule

%%%%%%%%%%%%%%%% Impose ZLB
%%%%%%%%%%%%%%%%%%%%%%

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

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



irf1 = sim_zlb2 - sim_zlb;

T = size(irf1(10:end,1),1);
t = 1:T;

% variable positions
pY = 49; % note this is log
pC = 50; 
pI = 51;
pL = 52;
pPi = 53;
pRd = 54;


figure
subplot(2,3,1)
plot(t,irf0(10:end,pY),'-k',t,irf1(10:end,pY),'--k','Linewidth',1.5)
title('Output')

subplot(2,3,2)
plot(t,irf0(10:end,pC),'-k',t,irf1(10:end,pC),'--k','Linewidth',1.5)
title('Consumption')

subplot(2,3,3)
plot(t,irf0(10:end,pI),'-k',t,irf1(10:end,pI),'--k','Linewidth',1.5)
title('Investment')

subplot(2,3,4)
plot(t,irf0(10:end,pL),'-k',t,irf1(10:end,pL),'--k','Linewidth',1.5)
title('Labor')

subplot(2,3,5)
plot(t,irf0(10:end,pPi),'-k',t,irf1(10:end,pPi),'--k','Linewidth',1.5)
title('\pi')

subplot(2,3,6)
plot(t,irf0(10:end,pRd),'-k',t,irf1(10:end,pRd),'--k','Linewidth',1.5)
title('r^{d}')

% Do public QE shocks
irfshock = char('et','eb');

%%%%%%%%%%%%% No ZLB
%%%%%%%%%%%%%%%%%%
% 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);
irf0 = sim_MP2 - sim_MP; % IRF with no Taylor rule

%%%%%%%%%%%%%%%% Impose ZLB
%%%%%%%%%%%%%%%%%%%%%%

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

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



irf1 = sim_zlb2 - sim_zlb;

T = size(irf1(10:end,1),1);
t = 1:T;

% variable positions
pY = 49; % note this is log
pC = 50; 
pI = 51;
pL = 52;
pPi = 53;
pRd = 54;


figure
subplot(2,3,1)
plot(t,irf0(10:end,pY),'-k',t,irf1(10:end,pY),'--k','Linewidth',1.5)
title('Output')

subplot(2,3,2)
plot(t,irf0(10:end,pC),'-k',t,irf1(10:end,pC),'--k','Linewidth',1.5)
title('Consumption')

subplot(2,3,3)
plot(t,irf0(10:end,pI),'-k',t,irf1(10:end,pI),'--k','Linewidth',1.5)
title('Investment')

subplot(2,3,4)
plot(t,irf0(10:end,pL),'-k',t,irf1(10:end,pL),'--k','Linewidth',1.5)
title('Labor')

subplot(2,3,5)
plot(t,irf0(10:end,pPi),'-k',t,irf1(10:end,pPi),'--k','Linewidth',1.5)
title('\pi')

subplot(2,3,6)
plot(t,irf0(10:end,pRd),'-k',t,irf1(10:end,pRd),'--k','Linewidth',1.5)
title('r^{d}')
