%  New Keynesian Model with 1 HH: Non-Stockholders

close all;

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

var C N	B_P	Lambda R	W_P	Pi	D_P	Y	Pstar_P	s	MC_P	Gamma1	Gamma2	eps_a	eps_r;
varexo eta_a eta_r;

parameters nu upsilon beta theta epsilon phi_p phi_y sigma_a sigma_r rho_a rho_r C_ss N_ss	B_P_ss	Lambda_ss	R_ss	W_P_ss	Pi_ss	D_P_ss	Y_ss	Pstar_P_ss	s_ss	MC_P_ss	Gamma1_ss	Gamma2_ss	eps_a_ss	eps_r_ss;


%----------------------------------------------------------------
% 2. Parameters
%----------------------------------------------------------------

load parameters;

set_param_value('nu',nu);
set_param_value('upsilon',upsilon);
set_param_value('beta',beta);
set_param_value('theta',theta);
set_param_value('epsilon',epsilon);
set_param_value('phi_p',phi_p);
set_param_value('phi_y',phi_y);
set_param_value('sigma_a',sigma_a);
set_param_value('sigma_r',sigma_r);
set_param_value('rho_a',rho_a);
set_param_value('rho_r',rho_r);

set_param_value('C_ss',C_ss);
set_param_value('N_ss',N_ss);
set_param_value('B_P_ss',B_P_ss);
set_param_value('Lambda_ss',Lambda_ss);
set_param_value('R_ss',R_ss);
set_param_value('W_P_ss',W_P_ss);
set_param_value('Pi_ss',Pi_ss);
set_param_value('D_P_ss',D_P_ss);
set_param_value('Y_ss',Y_ss);
set_param_value('Pstar_P_ss',Pstar_P_ss);
set_param_value('s_ss',s_ss);
set_param_value('MC_P_ss',MC_P_ss);
set_param_value('Gamma1_ss',Gamma1_ss);
set_param_value('Gamma2_ss',Gamma2_ss);
set_param_value('eps_a_ss',eps_a_ss);
set_param_value('eps_r_ss',eps_r_ss);

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

model;

%-----------
% Households
%-----------

% Stockholders
%-------------

W_P = C*N^upsilon;                                                          % (1) EE in labour market
1/R = Lambda(+1)/Pi(+1);                                                    % (2) EE in bond market
beta*C = Lambda(+1)*C(+1);                                                      % (3) Stochastic discount rate of stockholders

%-------
% Firms
%-------

% Production
%-----------

MC_P = W_P/eps_a;                                                                     % (4) Marginal cost

% Price
%------
Pstar_P = (epsilon/(epsilon-1))*Gamma1/Gamma2;                                         % (5) Price setting
Gamma1 = MC_P*Y + theta*Lambda(+1)*Gamma1(+1)*(Pi(+1))^(epsilon+1);      % (6)
Gamma2 = Y + theta*Lambda(+1)*Gamma2(+1)*(Pi(+1))^(epsilon);             % (7)

% Aggregate Profit
%------------------
D_P = Y*(theta*Pi^(epsilon-1) + (1-theta)*(Pstar_P)^(1-epsilon)) - Y*s*MC_P; % (8) Aggregate profit

%---------------------
% Aggregate Price Dynamic
%---------------------

Pstar_P = ((1-theta*Pi^(epsilon-1))/(1-theta))^(1/(1-epsilon));              % (9)            
s = theta*Pi^(epsilon) + (1-theta)*(Pstar_P)^(-epsilon);                     % (10) Price dispersion

%---------------------
% Market Clearing
%---------------------

Y = C;                                                                       % (11) Aggregate demand
s*Y = eps_a*N;                                                               % (12) Goods market clearing 
B_P = 0 ;                                                                    % (13) Bond market clearing
%C = W_P*N + D_P;                                                             % (13) Aggregte BC

%---------------------
% Monetary Policy Rule
%---------------------

R = R_ss*eps_r*Pi^phi_p;                                                       % (14) MP Rule              

%--------
% Shocks
%--------

log(eps_a) = rho_a*log(eps_a(-1))+eta_a;                                     % (15)
log(eps_r) = rho_r*log(eps_r(-1))+eta_r;                                     % (16)

end;

%----------------------------------------------------------------
% 4. Steady State
%----------------------------------------------------------------

initval;
C	=	C_ss	;
N           =	N_ss	;
B_P     =	B_P_ss	;
Lambda	=   Lambda_ss	;
R       =	R_ss	;
W_P     =	W_P_ss	;
Pi      =	Pi_ss	;
D_P     =	D_P_ss	;
Y       =	Y_ss	;
Pstar_P	=	Pstar_P_ss	;
s       =	s_ss	;
MC_P	=	MC_P_ss	;
Gamma1	=	Gamma1_ss	;
Gamma2	=	Gamma2_ss	;
eps_a	=	eps_a_ss	;
eps_r	=	eps_r_ss	;
end;

resid(1);
steady(solve_algo=3,maxit=100000);
check;

%----------------------------------------------------------------
% 5. Shocks
%----------------------------------------------------------------

shocks;
var eta_a = sigma_a^2;
var eta_r = sigma_r^2;
end;

%----------------------------------------------------------------
% 6. Calculation and Output
%----------------------------------------------------------------

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

output = Y_eta_g;