clc;clear;close all
cd('model_nonlinear')

Pis = 1; %% SS Trend inflation
betta=0.995; % Discount factor for Parent
Rss = Pis/betta; % SS short rate (2% at an annualized frequency)
Rres= Rss; % SS rate on reserves
sp = 2/400; % Targeted spread (2% annually)
Rbs = sp + Rss; %% SS long rate
bettab=betta * Rss/Rbs; % Discount factor for Child
epsip=11; % ss markup 10 percent
pms = (epsip-1)/epsip;  %% SS markup
ws = pms; %% SS wage
Ls = 1; % targeted SS labor
phip=0.75; %% Calvo price
pstars= (( 1-phip*Pis^(epsip-1))/(1-phip))^(1/(1-epsip)); %% SS reset price
vps = (1-phip)*pstars^(-epsip)/(1-phip*Pis^(epsip)) ; %% SS price dispersion
Ys = Ls/vps; % SS output
z= 0.33; % Targeted share of impatient types
Cbs = z * Ys; % SS consumption impatient
Cps = (1-z)*Ys; % SS counsumption patient
psi = ws / Cps; % 1.3636, Labor disutility scaling parameter (target L=1)
% psi = Cps^(-sigma)*ws/Ls^(chi);
sigm=1; % EIS
chi=1; % inverse Frisch elasticity

r = 0.1; % SS reserve/output ratio
res=r * Ys; % SS reserves
Omegas= betta* (Rbs-Rss);
kappa = 1-1/40; % tean year duration of the long bond
Qs = (Rbs-kappa)^(-1); %% SS long bond price
bs = Cbs / Qs; %% SS Total bond holdings
bCBs = res / Qs; % SS CB bond holdings 
bFIs = bs - bCBs; %% SS FI bond holdings
QEs= Qs * bCBs; % CB's balance sheet size 10% of output

Thetas = 5; % ss target of the risk-weighted leverage ratio
bar_XFIs = Qs * bFIs/Thetas; % Fixed equity of FI (0.0467) 
XFIs = bar_XFIs + kappa*Qs*bFIs; % Transfer to FI
Xbs = (Rbs* Qs*Pis^(-1)*bs); % Transfer to Children
ss= Qs*bFIs*(1-kappa*Pis^(-1)) + res - bar_XFIs; %% SS deposits/savingQEs= Qs * bCBs; % CB's balance sheet size 10% of output

% Checks
Ds = Ys - ws*Ls; % Production dividends
DFIs = Rbs*Qs*bFIs + Rres*res - Rss*ss; % FI dividends
Ts = - Rbs*Qs*bCBs + Rres*res; % Lump-sum transfer (should be negative bc it is a transfer)


% Parent's budgent constraint
check = Cps + ss - ws*Ls - Rss*ss-Ds - DFIs + Ts + XFIs + Xbs; % should be zero
check2 = Cbs + bs - Qs*bs*(1-kappa) - Xbs;
check3 = Cps + ss - (ws*Ls +  Ds)+ XFIs;

% Taylor rule
rhoR = 0.8;
phipi = 1.5;
phix = 0;
sr=-0.01; % (-) sign to generate positive MP shock (which is lowering interest rate)

% QE rule
rhoQE=0.8;
sq=0.01;
thetapi = -0;
thetax = -0;

% Shock process
rhoT=0.8;
rhoA=0.8;
sa=0.01;
st=0.01;

% save parameters
gamma = (1-phip)*(1-phip*betta)/phip; % Elasticity of inflation wrt real marginal cost (NKPC slope)
zeta = (chi*(1-z)+sigm)/(1-z); % Elasticity of real marginal cost wrt output gap (NKPC slope)
b1 = bFIs/bs;
b2 = bCBs/bs;


save param_nonlinear_exo  sigm chi Ls Ys Cps Cbs z betta bettab Rss Rres Rbs sp kappa Qs QEs bs bCBs bFIs b1 b2 ...
                        Thetas XFIs bar_XFIs Xbs Omegas res ss ...
                        epsip Pis vps pstars phip pms ws psi ...
                        rhoR phipi phix sr rhoQE  rhoT rhoA thetapi thetax sq sa st Ds DFIs Ts


% dynare linear_exo.mod nostrict
%% Normalize shocks to generate the same impact increase in output
xr = y_ea(1)/y_er(1);
xq = y_ea(1)/y_eq(1);
xt = y_ea(1)/y_et(1);
check1 = ystar_ea(1);
Yirf0= [y_ea xr*y_er xq*y_eq xt*y_et];
Piirf0= [pi_ea xr*pi_er xq*pi_eq xt*pi_et];
Xirf0 = [x_ea xr*x_er xq*x_eq xt*x_et];
Rirf0 =[rs_ea xr*rs_er xq*rs_eq xt*rs_et];
Exr0 = [exr_ea xr*exr_er xq*exr_eq xt*exr_et];
T_et = xt*theta_et;
QE_eq = xq* qe_eq;

T= size(y_ea,1);
t=1:T;

figure
plot(t,400*Exr0(:,2),'-k',t,400*Exr0(:,3),'--b','Linewidth',2)
grid on
xlabel('Horizon')
ylabel('Excess Return')
legend('MP Shock', 'QE Shock')
xlim([1 21])
xticks([1 6 11 16 21])
xticklabels({'0','5','10','15','20'})
saveas(gcf,'four_spreads','epsc') 









