% Simple New Keynesian Model

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var MULT Y C Lh L Kh Ks Pstar_P Wh_P Wstar_P W_P R Rk_P pi gammaw_1 gammaw_2 sw MC_P gammap_1 gammap_2 sp eps_b eps_a eps_R eps_l;
varexo eta_b eta_a eta_R eta_l;

parameters delta Ybar Rbar lambda_w lambda_p sigma_c lambda sigma_l theta_p theta_w PHI alpha beta rho_R phi_pi phi_y rho_a rho_b rho_r rho_l sigma_a sigma_b sigma_r sigma_l MULT_ss Y_ss C_ss Lh_ss L_ss Kh_ss Ks_ss Wh_P_ss Wstar_P_ss W_P_ss R_ss Rk_P_ss gammaw_1_ss gammaw_2_ss MC_P_ss gammap_1_ss gammap_2_ss r_ss premium_ss;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

load parameters;

set_param_value('delta' , delta);
set_param_value('lambda_w' , lambda_w);
set_param_value('lambda_p' , lambda_p);
set_param_value('sigma_c' , sigma_c);
set_param_value('lambda' , lambda);
set_param_value('sigma_l' , sigma_l);
set_param_value('theta_p' , theta_p);
set_param_value('theta_w' , theta_w);
set_param_value('PHI' , PHI);
set_param_value('alpha' , alpha);
set_param_value('beta' , beta);
set_param_value('rho_R' , rho_R);
set_param_value('phi_pi' , phi_pi);
set_param_value('phi_y' , phi_y);
set_param_value('rho_a' , rho_a);
set_param_value('rho_b' , rho_b);
set_param_value('rho_r' , rho_r);
set_param_value('sigma_a' , sigma_a);
set_param_value('sigma_b' , sigma_b);
set_param_value('sigma_r' , sigma_r);
set_param_value('Ybar' , Ybar);
set_param_value('Rbar' , Rbar);

set_param_value('MULT_ss',MULT_ss); 
set_param_value('Y_ss',Y_ss);
set_param_value('C_ss',C_ss);
set_param_value('Lh_ss',Lh_ss);
set_param_value('L_ss',L_ss);
set_param_value('Kh_ss',Kh_ss); 
set_param_value('Ks_ss',Ks_ss); 
set_param_value('Wh_P_ss',Wh_P_ss);
set_param_value('Wstar_P_ss',Wstar_P_ss);
set_param_value('W_P_ss',W_P_ss);
set_param_value('R_ss',R_ss);
set_param_value('Rk_P_ss',Rk_P_ss);
set_param_value('gammaw_1_ss',gammaw_1_ss); 
set_param_value('gammaw_2_ss',gammaw_2_ss);
set_param_value('MC_P_ss',MC_P_ss); 
set_param_value('gammap_1_ss',gammap_1_ss); 
set_param_value('gammap_2_ss',gammap_2_ss); 
sigma_l = 0.24;
rho_l = 0.12;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model;
% Households
MULT = eps_b/C; % FOC wrt C
Wh_P = (1/eps_b)*C*(Lh)^sigma_l;                                                               % intratemporal eq. in the labor market
R^(-1) = beta*(MULT(+1)/MULT)*(pi(+1))^(-1);                                  % intertemporal eq. in the bond market
MULT = beta*MULT(+1)*(1-delta+Rk_P(+1));                                            % intertemporal eq. in the capital market

% Unions 
(1/lambda_w)*(Wstar_P)*gammaw_1 = ((1+lambda_w)/lambda_w)*gammaw_2;
gammaw_1 = (W_P^((1+lambda_w)/lambda_w))*L+beta*theta_w*MULT(+1)/MULT*(pi(+1))^(1/lambda_w)*gammaw_1(+1);
gammaw_2 = (W_P^((1+lambda_w)/lambda_w))*Wh_P*L+beta*theta_w*MULT(+1)/MULT*(pi(+1))^((1+lambda_w)/lambda_w)*gammaw_2(+1);

% Firms
Ks/L = (alpha/(1-alpha))*W_P/Rk_P;
MC_P = ((alpha/(1-alpha))*W_P/Rk_P)^(-alpha)*(W_P/(1-alpha)*eps_a);
(1/lambda_p)*Pstar_P*gammap_1 = ((1+lambda_p)/lambda_p)*gammap_2;
gammap_1 = Y+beta*theta_p*MULT(+1)/MULT*(pi(+1))^(1/lambda_p)*gammap_1(+1);
gammap_2 = MC_P*Y+beta*theta_p*MULT(+1)/MULT*(pi(+1))^(1/lambda_p)*gammap_2(+1);

% Aggregate Price Dynamics
1 = (1-theta_p)*(Pstar_P)^(-1/lambda_p)+theta_p*pi^(1/lambda_p);
1 = (1-theta_w)*(Wstar_P*(1/W_P))^(-1/lambda_w)+theta_w*((1/pi)*W_P(-1)*(1/W_P))^(1/lambda_w);

% Central Bank
R/Rbar = ((R(-1)/Rbar)^rho_R)*((pi^phi_pi)*(Y/Ybar)^phi_y)^(1-rho_R)*eps_R;

% Market clearing and feasibility
Y = C+Kh-(1-delta)*Kh(-1);
sp*(C+Kh-(1-delta)*Kh(-1)) = eps_a*Ks^(alpha)*L^(1-alpha)-PHI;
sp = (1-theta_p)*Pstar_P^(-(1+lambda_p)/lambda_p)+theta_p*pi^((1+lambda_p)/lambda_p)*sp(-1);
Lh = L*sw;
sw = (1-theta_w)*(Wstar_P*W_P^(-1))^(-(1+lambda_w)/lambda_w)+theta_w*(pi^(-1)*W_P(-1)*W_P^(-1))^(-(1+lambda_w)/lambda_w)*sw(-1);
Ks = Kh(-1);

log(eps_b) = rho_b*log(eps_b(-1))+eta_b;
log(eps_a) = rho_a*log(eps_a(-1))+eta_a;
log(eps_R) = rho_r*log(eps_R(-1))+eta_R;
log(eps_l) = rho_l*log(eps_l(-1))+eta_l;

end;

initval;
MULT=MULT_ss; 
Y=Y_ss;
C=C_ss;
Lh=Lh_ss;
L=L_ss;
Kh=Kh_ss; 
Ks=Ks_ss; 
Pstar_P=1;
Wh_P=Wh_P_ss;
Wstar_P=Wstar_P_ss;
W_P=W_P_ss;
R=R_ss;
Rk_P=Rk_P_ss;
pi=1; 
gammaw_1=gammaw_1_ss; 
gammaw_2=gammaw_2_ss;
sw=1; 
MC_P=MC_P_ss; 
gammap_1=gammap_1_ss; 
gammap_2=gammap_2_ss; 
sp=1; 
eps_b=1; 
eps_a=1; 
eps_R=1;
eps_l=1;
end;
steady(solve_algo=4,maxit=10000);
check;

shocks;
var eta_b = sigma_b^2;
var eta_a = sigma_a^2;
var eta_R = sigma_r^2;
var eta_l = sigma_l^2;
end;

stoch_simul(hp_filter = 1600, order=1, irf=80);