%% Dynare mod file for solving a simple (non-linearized) NK model

clc;
close all;

%% Variable and Parameter Definition

var Q C Pi N WP Y A MPN Mu Pistar V1 V2 i R Nu S Ynatural Ytilde;

varexo eps_A eps_Nu;

parameters beta sigma varphi alpha epsilon theta phi_Pi phi_Y rho_A rho_Nu;

%parameters
beta    =0.99;
sigma   = 1;
varphi  = 1;
alpha   = 1/3; 
epsilon = 6;
theta   = 2/3;
phi_Pi  = 1.5;
phi_Y  = 0.5/4;
rho_A   = 0.9;
rho_Nu  = 0.3;

%steady state
Q_ss        = beta;
Nu_ss       = 1;
S_ss        = 1;
Pistar_ss   = 1;
Pi_ss       = 1;
A_ss        = 1;
i_ss        = (1/beta);
R_ss        = i_ss;
Mu_ss       = epsilon/(epsilon-1);
N_ss        = ((1-alpha)/Mu_ss)^(1/(sigma*(1-alpha)+varphi+alpha));
Y_ss        = N_ss^(1-alpha);
C_ss        = Y_ss;
WP_ss       = Y_ss^sigma*N_ss^varphi;
MPN_ss      = (1-alpha)*N_ss^(-alpha);
V1_ss       = Y_ss/(1-beta*theta);
V2_ss       = V1_ss/Mu_ss;
Ynatural_ss = Y_ss;
Ynatural_ss2 = ((1-alpha)*(epsilon-1/epsilon))^(1/(sigma+((varphi+alpha)/(1-alpha))));


%% The model (typed in non-linearly)

model;

// 1: Labor market equilibrium
(exp(N)^varphi)*(exp(C)^(-sigma))    =  exp(WP);  
 
// 2: Euler equation
exp(Q) = beta*((exp(C(+1))/exp(C))^(-sigma))*exp(Pi(+1))^(-1);

// 3: Arbitrage
exp(Q) = exp(i)^(-1);

// 4: Production Function
exp(S)*exp(Y)     =   exp(A)*(exp(N)^(1-alpha));

// 5: Natural Output
(1-alpha)*(1/(epsilon/(epsilon-1))) = (exp(Ynatural))^(sigma+((varphi+alpha)/(1-alpha))) * (exp(A)^(-(varphi+1)/(1-alpha)));

// 6: Output Gap
exp(Ytilde) = exp(Y)/exp(Ynatural);

// 7: Markup
(1-alpha)*(1/exp(Mu)) = exp(Y)^(sigma+((varphi+alpha)/(1-alpha))) * exp(A)^(-((varphi+1)/(1-alpha))) * exp(S)^((varphi+1)/(1-alpha));

// 8: Marginal Product of Labor
exp(MPN) = (1-alpha)*exp(A)/exp(S)*exp(N)^(-alpha);

// 9: Aggregate resource constraint
exp(Y)   =   exp(C);

// 10: F1: Optimal price choice
exp(Pistar)   =  epsilon/(epsilon-1)*exp(V2)/exp(V1);

// 11: F2: V1
exp(V1)       =   exp(Y)*exp(Pi)^(epsilon-1)+beta*theta*exp(Pi)^(epsilon-1)*(exp(C(+1))/exp(C))^(-sigma)*exp(V1(+1));
 
// 12: F3: V2
exp(V2)       =   (exp(Pi)^epsilon)*exp(Y)/exp(Mu)+beta*theta*exp(Pi)^epsilon*(exp(C(+1))/exp(C))^(-sigma)*exp(V2(+1));
 
// 13: Price Dynamics (with inflation indexation)
(exp(Pi))^(1-epsilon)     =   theta + (1-theta)*(exp(Pistar))^(1-epsilon);

// 14: Fisher equation
exp(i)  =   exp(R)*exp(Pi(+1));
 
// 15: Interest rate rule
exp(i)      =  ((1/beta)*exp(Pi)^phi_Pi*(exp(Ytilde))^(phi_Y))*exp(Nu);

// 16: Price dispersion
exp(S)    =  (1-theta)*exp(Pi)^(epsilon)*(exp(Pistar))^(-epsilon)+theta*exp(Pi)^epsilon*exp(S(-1));
 
// 17: TFP Shock
A  =   rho_A*A(-1)+eps_A;

// 18: Monetary Policy shock
Nu  =   rho_Nu*Nu(-1)+eps_Nu;

end;


initval;
Q       = log(Q_ss);
C       = log(C_ss);
Pi      = log(Pi_ss);
N       = log(N_ss);
WP      = log(WP_ss);
Y       = log(Y_ss);
A       = log(1);
MPN     = log(MPN_ss);
Mu      = log(Mu_ss);
Pistar  = log(Pistar_ss);
V1      = log(V1_ss);
V2      = log(V2_ss);
i       = log(i_ss);
R       = log(R_ss);
Nu      = log(Nu_ss);
S       = log(S_ss);
Ynatural = log(Ynatural_ss);
Ytilde   = log(1);
end;
 


shocks;
%var eps_A; stderr 1;
var eps_Nu; stderr 0.01;
end;


steady;

check;

stoch_simul(hp_filter = 1600, order = 1, irf=100,nograph);

%stoch_simul(hp_filter = 1600, order = 1, irf=12) Y Ynatural Ytilde  N Pi i R S A Nu;
%stoch_simul(hp_filter = 1600, order = 2, irf=12, pruning) Y Ynatural Ytilde  N Pi i R S A Nu;