var Y_t C_t I_t K_t N_t w_t rk_t R_t mc_t  Pi_ast g1_t g2_t Pi_t z_t lambda_t  D_t    q_t  U_t ;  
 %   1   2   3   4   5   6    7   8    9      10     11   12   13  14      15    16    17   18
varexo eps_z eps_U;

parameters  sigma eta beta delta alpha theta omega rho_z rho_r delta_pi delta_y chi 
            Pi_bar D_bar N_bar rk_bar R_bar mc_bar K_bar I_bar rho_U U_bar sigma_U
            Y_bar C_bar lambda_bar w_bar g1_bar g2_bar z_bar hab sigma_z q_bar;

% set parameter values
sigma       = 1;        % Consumption preferences
eta         = 1;        % frisch elasticity
beta        = 0.99;     % discount factor
delta       = 0.025;    % discount rate 
alpha       = 1/3;     % production technology from RBC
theta       = 10;        % mark-up parameter
omega       = 0.75;     % Calvo parameter
rho_z       = 0.95;      % technology shock persistence
rho_r       = 0.5;      % degree of interest rate smoothing
rho_U       = 0.88;
sigma_U     = 0.002;
delta_pi    = 1.5;      % CB stance on inflation
delta_y     = 0.5;      % CB stance on deflation from ss output
phi         = 0;
sigma_z=0.01;

% chi is endogenously determined



% calculation of steady state
Pi_bar      = 1.006;
N_bar       = 0.33;
z_bar       = 1;
U_bar=0.01;

rk_bar      = 1/beta -1+delta;
R_bar       = Pi_bar/beta;
Pi_ast_bar  = ((Pi_bar^(1-theta)-omega)/(1-omega))^(1/(1-theta)); 
D_bar       = (1-omega)*(Pi_bar/Pi_ast_bar)^theta/(1-omega*Pi_bar^theta);
g1g2        = (Pi_ast_bar/Pi_bar) * (theta-1)/theta;
mc_bar      = g1g2* ((1-omega*beta*Pi_bar^theta)/(1-omega*beta*Pi_bar^(theta-1)));
K_bar       = (alpha*mc_bar/rk_bar)^(1/(1-alpha))*N_bar;
I_bar       = delta*K_bar;
Y_bar       = K_bar^alpha*N_bar^(1-alpha)/D_bar;
C_bar       = Y_bar - I_bar;
lambda_bar  = ((1-hab)*C_bar)^(-sigma);
w_bar       = (1-alpha)*mc_bar*Y_bar/N_bar;
chi         = w_bar*lambda_bar/N_bar^eta;
g1_bar      = lambda_bar*Y_bar*mc_bar/(1-omega*beta);
g2_bar      = lambda_bar*Y_bar/(1-omega*beta);
q_bar=lambda_bar;



model;

%eq 1 
exp(lambda_t) =exp(C_t)^(-sigma);

%eq 2 Labour Supply
exp(w_t) * exp(lambda_t) =chi*exp(N_t)^eta;

%eq 3 Capital Euler Equation
exp(q_t)=beta*(exp(lambda_t(+1))*(exp(rk_t(+1))-0.5*phi*(exp(I_t(+1))/exp(K_t)-delta)^2+phi*(exp(I_t(+1))/exp(K_t)-delta)*exp(I_t(+1))/exp(K_t))+exp(q_t(+1))*(1-delta));

%eq 4 foc w.r.t investment
exp(q_t)=exp(lambda_t)*(1+phi*(exp(I_t)/exp(K_t(-1))-delta));

%eq 5 Bonds Euler Equation
exp(lambda_t)*exp(Pi_t(+1)) = beta * exp(R_t) * exp(lambda_t(+1));


%------Production Sector---------------
%eq 6  Production Function
exp(Y_t)=exp(z_t)*exp(K_t(-1))^alpha * exp(N_t)^(1-alpha)/exp(D_t); 

%eq 7 Law of motion of capital
exp(K_t)=(1-delta)*exp(K_t(-1))+exp(I_t); 


%eq 8 Labour demand
exp(w_t)=(1-alpha)*exp(mc_t)*exp(Y_t)/exp(N_t);

%eq 9 Factor Price ratio (implies capital demand)
exp(rk_t) = alpha * exp(mc_t) * exp(Y_t) / exp(K_t(-1));


%--------Price Setting equations---------------------------------------

%eq 10 Nonlinear NK-Phillips-Curve (recursive) A
exp(Pi_ast)=exp(Pi_t)*(theta/(theta-1)) * exp(g1_t) / exp(g2_t);


%eq 11 Nonlinear NK-Phillips-Curve (recursive) B
exp(g1_t) = exp(lambda_t) * exp(Y_t) * exp(mc_t) + omega*beta * exp(Pi_t(+1))^theta * exp(g1_t(+1));


%eq 12 Nonlinear NK-Phillips-Curve (recursive) C
exp(g2_t) = exp(lambda_t) * exp(Y_t) + omega*beta * exp(Pi_t(+1))^(theta-1) * exp(g2_t(+1));


%eq 13  Price Dispersion
exp(D_t)=(1-omega)* (exp(Pi_t)/exp(Pi_ast))^theta+omega*exp(Pi_t)^theta*exp(D_t(-1));

%eq 14  Price evolution
exp(Pi_t)^(1-theta) = (1-omega)*exp(Pi_ast)^(1-theta) +omega;


%--------Closing the model-------------------------

%eq15 Goods Market Clearing
exp(Y_t)=exp(I_t)+exp(C_t)+phi/2*(exp(I_t)/exp(K_t(-1))-delta)^2*exp(K_t(-1));



%eq 16 Taylor Rule
(exp(R_t)/R_bar) = (exp(R_t(-1))/R_bar)^rho_r *  ((exp(Pi_t)/Pi_bar)^delta_pi * (exp(Y_t)/Y_bar)^delta_y )^(1-rho_r);


%eq 17 TFP process
z_t=rho_z*z_t(-1)+U_t*eps_z;

%eq 18 Uncertainty process
U_t=(1-rho_U)*U_bar+rho_U*U_t(-1)+sigma_U*eps_U;

end;


initval;
Y_t         = log(Y_bar);  
C_t         = log(C_bar);
I_t         = log(I_bar);
K_t         = log(K_bar);
N_t         = log(N_bar);
w_t         = log(w_bar);
rk_t        = log(rk_bar);
R_t         = log(R_bar);
mc_t        = log(mc_bar);
Pi_ast      = log(Pi_ast_bar);
g1_t        = log(g1_bar);
g2_t        = log(g2_bar);
D_t         = log(D_bar);
Pi_t        = log(Pi_bar);
z_t         = log(z_bar);
lambda_t    = log(lambda_bar);
q_t         = log(lambda_bar);
U_t         = U_bar;
end;


resid;

steady;

check;

shocks;
var eps_z=1;
var eps_U=1;

end;




//---------------------------------------------------------------------
// 6. Run Dynare with code from Gohde-Meyer and Han
//---------------------------------------------------------------------

stoch_simul(irf=20,order=3,nograph,noprint);
options_.irf=30;
my_nonlinear_MA_options.plot_simulations=0;
my_nonlinear_MA_options.shock_scale=1;
my_nonlinear_MA_options.plot_irf=0;
[MA_]=nonlinear_MA(M_,oo_,options_,my_nonlinear_MA_options);

//---------------------------------------------------------------------

var Y_t C_t I_t K_t N_t w_t rk_t R_t mc_t  Pi_ast g1_t g2_t Pi_t z_t lambda_t  D_t    q_t  U_t ;  
 %   1   2   3   4   5   6    7   8    9      10     11   12   13  14      15    16    17   18

irfs=MA_.irf.total([1:8,13]+18,:)*100; 
ttitle={'Y_t', 'C_t' ,'I_t', 'K_t', 'N_t', 'w_t', 'rk_t', 'R_t', '\Pi_t'};
figure
for i=1:length(ttitle)
subplot(3,3,i)
plot(irfs(i,:));
title(ttitle(i));
end