%% 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 i R Nu Ynatural Ytilde;

% Pistar, V1, V2, S

varexo eps_A eps_Nu;

parameters beta sigma varphi alpha epsilon theta phi_Pi phi_Y rho_A rho_Nu 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;

nu = 1;

%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(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)));
// 8: Marginal Product of Labor
exp(MPN) = (1-alpha)*exp(A)*exp(N)^(-alpha);
// 9: Aggregate resource constraint
exp(Y)*(1-nu/2*(exp(Pi)-1))^2   =   exp(C);
// 10: Fisher equation
exp(i)  =   exp(R)*exp(Pi(+1));
// 11: Interest rate rule
exp(i)      =  ((1/beta)*exp(Pi)^phi_Pi*(exp(Ytilde))^(phi_Y))*exp(Nu);
// 12: Rotemberg Pricing
nu*(exp(Pi)-1) = (1-epsilon) + epsilon*(1/exp(Mu))+nu*beta*((exp(C(+1))/exp(C))^(-sigma))*(exp(Y(+1))/exp(Y))*(exp(Pi(+1))-1)*exp(Pi(+1));
// 13: TFP Shock
A  =   rho_A*A(-1)+eps_A;
// 14: 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);
%S       = log(S_ss);

i       = log(i_ss);
R       = log(R_ss);
Nu      = log(Nu_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 A Nu;
%stoch_simul(hp_filter = 1600, order = 2, irf=12, pruning) Y Ynatural Ytilde  N Pi i R S A Nu;