%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Please cite as follows:
%Eric Sims and Jing Cynthia Wu "Evaluating Central Banks' Tool Kit: Past, 
%Present, and Future", Journal of Monetary Economics, forthcoming.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning off;

var 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 welf;

varexo eA eG et eB eb ef ec el em eD er;

parameters  alpha beta delta0 delta1 delta2 bh chi eta sp X kappaf kappab 
            sigma thetas Deltas Pis epsiw epsip phiw gammaw gammap 
            kappai Gs fcbs bcbs res rotemberg rhor phipi phiy rhoA rhom rhoc rhol 
            rhoG rhoB rhof rhob rhot rhoD sr sA sm sc sl sG sB sf sb st sD 
            bGs levs phis Res psif psib psi ;

load param;

verbatim;
for ii=1:length(M_.params)
    set_param_value(deblank(M_.param_names(ii,:)),eval(deblank(M_.param_names(ii,:))));
end
end;

model;

%%%%%%%%%%%%%% Financial intermediary %%%%%%%%%%%%%%
% (1) Return on private bond
exp(RF) = (1+kappaf*exp(QF))/exp(QF(-1));

% (2) Return on government bond
exp(RB) = (1 + kappab*exp(QB))/exp(QB(-1));

% (3) Extra discounting
Omega = 1 - sigma + sigma*a;

% (4) FOC private bonds
(exp(RF(+1)) - exp(Rd))*exp(Lam(+1))*Omega(+1)*exp(infl(+1))^(-1) = (lam/(1+lam))*exp(theta);

% (5) FOC government bonds
(exp(RB(+1)) - exp(Rd))*exp(Lam(+1))*Omega(+1)*exp(infl(+1))^(-1) = (lam/(1+lam))*exp(theta)*exp(Delta);

% (6) FOC on reserves
(exp(RE) - exp(Rd))*exp(Lam(+1))*Omega(+1)*exp(infl(+1))^(-1) = -omega/(1+lam);

% (7) Coefficient a in value function / Omega
a = exp(theta)*exp(phi);

% (8-10) Modified leverage term
phi_1 = (exp(Rd)*exp(infl(+1))^(-1)*exp(Lam(+1))*Omega(+1));
phi_2 = (exp(theta) - (exp(RF(+1)) - exp(Rd))*exp(infl(+1))^(-1)*exp(Lam(+1))*Omega(+1));
exp(phi) = phi_1/phi_2 - omega * re/exp(n)/exp(theta);

%%%%%%%%%%%%%% Household %%%%%%%%%%%%%%
% (11) marginal utility
exp(mu) = exp(psic)/(exp(C) - bh*exp(C(-1))) - beta*bh*exp(psic(+1))/(exp(C(+1)) - bh*exp(C));

% (12) FOC labor
exp(psic)*exp(psil)*chi*exp(L)^(eta) = exp(mu)*exp(mrs);

% (13) stochastic discount factor
exp(Lam) = beta*exp(mu)/exp(mu(-1));

% (14) FOC deposits
1 = exp(Lam(+1))*exp(Rd)*exp(infl(+1))^(-1);

%%%%%%%%%%%%%% Labor union %%%%%%%%%%%%%%
% (15) Reset wage
exp(wr) = (epsiw/(epsiw-1))*f1/f2;

% (16) Auxiliary wage-setting f1
f1 = exp(mrs)*exp(w)^(epsiw)*exp(Ld) + phiw*exp(Lam(+1))*exp(infl)^(-epsiw*gammaw)*exp(infl(+1))^(epsiw)*f1(+1);

% (17) Auxiliary wage-setting f2
f2 = exp(w)^(epsiw)*exp(Ld) + phiw*exp(Lam(+1))*exp(infl)^((1-epsiw)*gammaw)*exp(infl(+1))^(epsiw-1)*f2(+1);

%%%%%%%%%%%%%% Retail firm %%%%%%%%%%%%%%


%%%%%%%%%%%%%% Wholesale firm %%%%%%%%%%%%%%
% (21) Production
exp(Ym) = exp(A)*(exp(u)*exp(K(-1)))^(alpha)*exp(Ld)^(1-alpha);

% (22) Labor demand 
exp(w) = exp(pm)*(1-alpha)*exp(A)*(exp(u)*exp(K(-1)))^(alpha)*exp(Ld)^(-alpha);

% (23) Utilization
exp(pk)*exp(M)*(delta1 + delta2*(exp(u) - 1)) = exp(pm)*alpha*exp(A)*(exp(u)*exp(K(-1)))^(alpha-1)*exp(Ld)^(1-alpha);

% (24) Capital FOC
exp(pk)*exp(M) = exp(Lam(+1))*(alpha*exp(pm(+1))*exp(A(+1))*(exp(u(+1))*exp(K))^(alpha-1)*exp(u(+1))*exp(Ld(+1))^(1-alpha) + (1-delta0-delta1*(exp(u(+1)) - 1)-(delta2/2)*(exp(u(+1))-1)^2)*exp(pk(+1))*exp(M(+1)));

% (25-26) Bonds FOC
exp(QF)*exp(M2) = exp(Lam(+1))*exp(infl(+1))^(-1)*(1+kappaf*exp(QF(+1))*exp(M2(+1)));

%(exp(M)-1)/(exp(M2)-1) = psi;
exp(M) - 1 = psi*(exp(M2) - 1);

% (27) Loan constraint
psi*exp(pk)*exp(Ih) = exp(QF)*(exp(fpr) - kappaf*exp(infl)^(-1)*exp(fpr(-1)));

% (28) Capital accumulation
exp(K) = exp(Ih) + (1-delta0-delta1*(exp(u)-1)-(delta2/2)*(exp(u)-1)^2)*exp(K(-1));

%%%%%%%%%%%%%% Capital producer %%%%%%%%%%%%%%
% (29) Production function
exp(Ih) = exp(mui)*(1 - (kappai/2)*(exp(I)/exp(I(-1)) - 1)^2)*exp(I);

% (30) Investment FOC
1 = exp(pk)*exp(mui)*(1 - (kappai/2)*(exp(I)/exp(I(-1)) - 1)^2 - kappai*(exp(I)/exp(I(-1))-1)*(exp(I)/exp(I(-1)))) + exp(Lam(+1))*exp(pk(+1))*exp(mui(+1))*kappai*(exp(I(+1))/exp(I) - 1)*(exp(I(+1))/exp(I))^2;

%%%%%%%%%%%%%% Government %%%%%%%%%%%%%%
% (31) Budget constraint
exp(G) + exp(bG(-1))*exp(infl)^(-1) = T + Tcb + exp(QB)*(exp(bG) - kappab*exp(infl)^(-1)*exp(bG(-1)));

%%%%%%%%%%%%%% Central bank %%%%%%%%%%%%%%
% (32) Balance sheet condition
exp(QF)*fcb + exp(QB)*bcb = re;

% (33) Transfer to fiscal authority
Tcb = (1+kappaf*exp(QF))*exp(infl)^(-1)*fcb(-1) + (1+kappab*exp(QB))*exp(infl)^(-1)*bcb(-1) - exp(RE(-1))*exp(infl)^(-1)*re(-1);

% (34-36) Taylor rule
RTR = (1-rhor)*log(Res) + rhor*RTR(-1) + (1-rhor)*(phipi*(infl - log(Pis)) + phiy*(Y - Y(-1))) + sr*er;
RE  = 0;
Rd = RE;

%%%%%%%%%%%%%% Aggregate conditions %%%%%%%%%%%%%%
% (37) Wholesale vs retail output
exp(Ym) = exp(Y);


% (39) Price evolution


%% ROTEMBERG/ PRICE EVOLUTION


%(infl-1)*infl = 1/rotemberg*(epsip*pm*Y)+1-epsip)+beta*Lam*Y(+1)/Y*infl(+1)*(infl(+1)-1);

epsip-1=epsip*exp(pm)-rotemberg*exp(infl)*exp(infl)-1+exp(Lam)*rotemberg*exp(infl(+1))*exp((infl(+1))-1);


% (40) Labor demand vs labor supply
exp(L) = exp(Ld)*vw;

% (41) Wage dispersion
vw = (1-phiw)*(exp(wr)/exp(w))^(-epsiw) + phiw*exp(infl(-1))^(-gammaw*epsiw)*exp(infl)^(epsiw)*exp(w)^(epsiw)*exp(w(-1))^(-epsiw)*vw(-1);

% (42) Wage evolution
exp(w)^(1-epsiw) = (1-phiw)*exp(wr)^(1-epsiw) + phiw*exp(infl(-1))^(gammaw*(1-epsiw))*exp(infl)^(epsiw-1)*exp(w(-1))^(1-epsiw);

% (43) Financial intermediary balance sheet condition
exp(QF)*exp(f) + exp(QB)*exp(b) + re = exp(d) + exp(n);

% (44) Binding constraint
exp(QF)*exp(f) + exp(Delta)*exp(QB)*exp(b) = exp(phi)*exp(n);

% (45) Evolution of net worth
exp(n) = sigma*exp(infl)^(-1)*((exp(RF) - exp(Rd(-1)))*exp(QF(-1))*exp(f(-1)) + (exp(RB) - exp(Rd(-1)))*exp(QB(-1))*exp(b(-1)) + (exp(RE(-1)) - exp(Rd(-1))) * re(-1) +  exp(Rd(-1))*exp(n(-1))) + X; 

% (46) Private bond market-clearing
exp(fpr) = exp(f) + fcb;

% (47) Government bond market-clearing
exp(bG) = exp(b) + bcb;

% (48) Resource constraint
exp(Y) = exp(C) + exp(I) + exp(G);

% (49) Yield on private bond
exp(RLF) = exp(QF)^(-1) + kappaf;

% (50) Yield on government bond
exp(RLB) = exp(QB)^(-1) + kappab;

% (51) Hypothetical expectations hypothesis bond
exp(QEH) = (1+kappab*exp(QEH(+1)))/exp(Rd);

% (52) Yield on expectations hypothesis bond
exp(RLEH) = exp(QEH)^(-1) + kappab;

%%%%%%%%%%%%%% Exogenous processes %%%%%%%%%%%%%%
% (53) Productivity
A = rhoA*A(-1) + sA*eA;

% (54) Investment
mui = rhom*mui(-1) + sm*em;

% (55) Intertemporal preference
psic = rhoc*psic(-1) + sc*ec;

% (56) Intratemporal preference
psil = rhol*psil(-1) + sl*el;

% (57) Government spending
G = (1-rhoG)*log(Gs) + rhoG*G(-1) + sG*eG;

% (58) Government debt
bG = (1-rhoB)*log(bGs) + rhoB*bG(-1) + sB*eB;

% (59) Central bank private bond holdings
fcb = (1-rhof)*fcbs + rhof*fcb(-1) + psif*(1-rhof)*(phipi*(infl - log(Pis)) + phiy*(Y - Y(-1))) + sf*ef;

% (60) Central bank government bond holdings
bcb = (1-rhob)*bcbs + rhob*bcb(-1) + psib*(1-rhob)*(phipi*(infl - log(Pis)) + phiy*(Y - Y(-1))) + sb*eb;

% (61) theta
theta = (1-rhot)*log(thetas) + rhot*theta(-1) + st*et;

% (62) Delta
Delta = (1-rhoD)*log(Deltas) + rhoD*Delta(-1) + sD*eD;

%%%%%%%%%%%%%% speards %%%%%%%%%%%%%%
% (63) returnf
returnf = exp(RF) - exp(Rd(-1));

% (64) returnb
returnb = exp(RB) - exp(Rd(-1));

% (65) returnfb
returnfb = exp(RF) - exp(RB);

% (66) spread
spreadf = RLF - Rd;

% (67)
spreadb = RLB - Rd;

%%%%%%%%%%%%%% real rates %%%%%%%%%%%%%%
% (68)
rlb_real = RLB - infl(+1);

% (69)
rlf_real = RLF - infl(+1);

% (70)
YY = exp(Y);

% (71) 
GG = exp(G);

% (72)
dY = (Y - Y(-1));

%%%%%%%%%%%%%% welfare %%%%%%%%%%%%%%
% (73) Welfare
welf = log(exp(C) - bh*exp(C(-1))) - psi*exp(L)^(1+eta)/(1+eta) + beta*welf(+1);

end;

options_.noprint=1;

steady;

shocks;
var eA = 1;
var eG = 1;
var et = 1;
var eB = 1;
var eb = 1;
var ef = 1;
var ec = 1;
var el = 1;
var em = 1;
var eD = 1;
var er = 1;
end;

stoch_simul(order=1,irf=20,nograph,ar=0,qz_zero_threshold=1e-8);