% NKSOE Simple Model

%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Preamble
%----------------------------------------------------------------

// Endogenous variables (27)
var c la h w ppi ppi_S R R_ast; //8
var c_H c_F p_H p_F; //4
var y_H mc_H f1H f2H p_tilde_H ppi_H rer; //7
var c_H_ast d_ast tb m d_M s_tb s_d_ast vel; //8

// Exogenous state variables (5)
var z R_W y_ast ppi_ast e_m;

// Exogenous innovations (5)
varexo u_z u_R_W u_y_ast u_ppi_ast u_e_m;//

// Parameters 
parameters BETA SIGMA PSI VARPHI OMEGA ETA EPSILON_H THETA_H;
parameters PHI_d ETA_ast BAR_d PHI_m NU; 
parameters ALPHA_ppi ALPHA_y ALPHA_R ALPHA_ppiS;
parameters RHO_z RHO_R_W RHO_y_ast RHO_ppi_ast RHO_e_m;
parameters ppi_ss R_ss y_H_ss h_ss s_tb_ss p_H_ss ppi_S_ss s_d_ast_ss vel_ss;
parameters z_ss R_W_ss y_ast_ss ppi_ast_ss e_m_ss;
parameters GAMMA_ppi_H ppi_H_ss

// Calibrated parameters (Endogenous BETA ppi_ast_ss R_W_ss PSI y_ast_ss BAR_d NU)
GAMMA_ppi_H=1.3;
SIGMA=1;
VARPHI=1;
OMEGA=0.65;
ETA=1.3;
EPSILON_H=11;
THETA_H=0.7;
PHI_d=0.001;
ETA_ast=0.3;
ALPHA_R=0.7;
ALPHA_ppi=1.5;
ALPHA_y=0.2;
ALPHA_ppiS=0;
PHI_m=0.3;
RHO_z=0.95;
RHO_R_W=0.95;
RHO_y_ast=0;
RHO_ppi_ast=0.9;
RHO_e_m=0;

// Targeted steady state values
ppi_ss=1.018;
h_ss=1/3;
p_H_ss=1;
R_ss=1.019;
ppi_S_ss=0.99;
s_tb_ss=0.05;
vel_ss=3;
z_ss=1;
e_m_ss=1;

%----------------------------------------------------------------
% 2. Model (32=5+3+5+3+7+1+5+3 equations)
%----------------------------------------------------------------

model; 

// Households (5)
/*
la=c^(-SIGMA); //E1
la*w=PSI*h^VARPHI; //E2
1=R*BETA*la(+1)/la/ppi(+1); //E3
1=R_ast*BETA*la(+1)/la/ppi(+1)*ppi_S(+1); //E4
m=NU*y_H*R^(-PHI_m); //E5, this is an exogenous money demand (not derived from money in the utility)
*/
exp(la)=(exp(c))^(-SIGMA); //E1
exp(la)*exp(w)=PSI*exp(h)^VARPHI; //E2
1=exp(R)*BETA*exp(la(+1))/exp(la)/exp(ppi(+1)); //E3
1=exp(R_ast)*BETA*exp(la(+1))/exp(la)/exp(ppi(+1))*exp(ppi_S(+1)); //E4
exp(m)=NU*exp(y_H)*exp(R)^(-PHI_m); //E5

// Final consumption (3) 
/*
c=(OMEGA^(1/ETA)*c_H^(1-1/ETA)+(1-OMEGA)^(1/ETA)*c_F^(1-1/ETA))^(ETA/(ETA-1)); //E6
c_F=(1-OMEGA)*p_F^(-ETA)*c; //E7
c_H=OMEGA*p_H^(-ETA)*c; //E8
*/
exp(c)=(OMEGA^(1/ETA)*exp(c_H)^(1-1/ETA)+(1-OMEGA)^(1/ETA)*exp(c_F)^(1-1/ETA))^(ETA/(ETA-1)); //E6
exp(c_F)=(1-OMEGA)*exp(p_F)^(-ETA)*exp(c); //E7
exp(c_H)=OMEGA*exp(p_H)^(-ETA)*exp(c); //E8

// Domestic producers (5)
/*
p_H*mc_H*z=w; //E9
f1H=p_tilde_H^(1-EPSILON_H)*y_H*(EPSILON_H-1)/EPSILON_H+THETA_H*BETA*la(+1)/la/ppi(+1)*(p_tilde_H*ppi/p_tilde_H(+1))^(1-EPSILON_H)*ppi_H(+1)^(EPSILON_H)*f1H(+1); //E10
f2H=p_tilde_H^(-EPSILON_H)*y_H*mc_H+THETA_H*BETA*la(+1)/la/ppi(+1)*(p_tilde_H*ppi/p_tilde_H(+1))^(-EPSILON_H)*ppi_H(+1)^(EPSILON_H+1)*f2H(+1); //E11
f1H=f2H; //E12
1=THETA_H*(ppi_H/ppi(-1))^(EPSILON_H-1)+(1-THETA_H)*p_tilde_H^(1-EPSILON_H); //E13
*/
exp(p_H)*exp(mc_H)*exp(z)=exp(w); //E9
exp(f1H)=exp(p_tilde_H)^(1-EPSILON_H)*exp(y_H)*(EPSILON_H-1)/EPSILON_H+THETA_H*BETA*exp(la(+1))/exp(la)/exp(ppi(+1))*(exp(p_tilde_H)*exp(ppi)/exp(p_tilde_H(+1)))^(1-EPSILON_H)*exp(ppi_H(+1))^(EPSILON_H)*exp(f1H(+1)); //E10
exp(f2H)=exp(p_tilde_H)^(-EPSILON_H)*exp(y_H)*exp(mc_H)+THETA_H*BETA*exp(la(+1))/exp(la)/exp(ppi(+1))*(exp(p_tilde_H)*exp(ppi)/exp(p_tilde_H(+1)))^(-EPSILON_H) *exp(ppi_H(+1))^(EPSILON_H+1)*exp(f2H(+1)); //E11
exp(f1H)=exp(f2H); //E12
1=THETA_H*(exp(ppi_H)/exp(ppi(-1)))^(EPSILON_H-1)+(1-THETA_H)*exp(p_tilde_H)^(1-EPSILON_H); //E13

// Rest of the world (3)
/*
R_ast=R_W*exp(PHI_d*(d_ast-BAR_d)); //E14
p_F=rer; //E15
c_H_ast=(p_H/rer)^(-ETA_ast)*y_ast; //E16
*/
exp(R_ast)=exp(R_W)*exp(PHI_d*(d_ast-BAR_d)); //E14
exp(c_H_ast)=(exp(p_H)/exp(rer))^(-ETA_ast)*exp(y_ast); //E15
exp(p_F)=exp(rer); //E16

// Market clearing and aggregation (7)
/*
y_H=c_H+c_H_ast; //E17
y_H=z*h; //E18
p_H/p_H(-1)=ppi_H/ppi; //E19
rer/rer(-1)=ppi_S*ppi_ast/ppi; //E20
rer*d_ast+tb=rer*d_ast(-1)/ppi_ast*R_ast(-1); //E21
tb=p_H*c_H_ast-rer_F*c_F; //E22
d_M=m/m(-1)*ppi; //E23
*/
exp(y_H)=exp(c_H)+exp(c_H_ast); //E17
exp(y_H)=exp(z)*exp(h); //E18
exp(p_H)/exp(p_H(-1))=exp(ppi_H)/exp(ppi); //E19
exp(rer)/exp(rer(-1))=exp(ppi_S)*exp(ppi_ast)/exp(ppi); //E20
exp(rer)*d_ast+tb=exp(rer)*d_ast(-1)/exp(ppi_ast)*exp(R_ast(-1)); //E21
tb=exp(p_H)*exp(c_H_ast)-exp(rer)*exp(c_F); //E22
exp(d_M)=exp(m)/exp(m(-1))*exp(ppi); //E23

// Monetary policy (1)
/*
(ppi_S/ppi_S_ss)=(ppi_H/ppi_H_ss)^GAMMA_ppi_S; //E24
*/
exp(ppi_S)/ppi_S_ss)=(exp(ppi_H)/ppi_H_ss)^GAMMA_ppi_S; //E24

// Exogenous AR(1) processes (5)
/*
log(z/z_ss)=RHO_z*log(z(-1)/z_ss)+u_z; //E25
log(R_W/R_W_ss)=RHO_R_W*log(R_W(-1)/R_W_ss)+u_R_W; //E26
log(e_m/e_m_ss)=RHO_e_m*log(e_m(-1)/e_m_ss)+u_e_m; //E29
*/
log(exp(z)/z_ss)=RHO_z*log(exp(z(-1))/z_ss)+u_z; //E25
log(exp(R_W)/R_W_ss)=RHO_R_W*log(exp(R_W(-1))/R_W_ss)+u_R_W; //E26
log(exp(y_ast)/y_ast_ss)=RHO_y_ast*log(exp(y_ast(-1))/y_ast_ss)+u_y_ast; //E27
log(exp(ppi_ast)/ppi_ast_ss)=RHO_ppi_ast*log(exp(ppi_ast(-1))/ppi_ast_ss)+u_ppi_ast; //E28
log(exp(e_m)/e_m_ss)=RHO_e_m*log(exp(e_m(-1))/e_m_ss)+u_e_m; //E29

// Definitions (3)
/*
s_tb=tb/(p_H*y_H); //E30
s_d_ast=rer*d_ast/(p_H*y_H); //E31
vel=p_H*y_H/m; //E32
*/
s_tb=tb/(exp(p_H)*exp(y_H)); //E30
s_d_ast=exp(rer)*d_ast/(exp(p_H)*exp(y_H)); //E31
exp(vel)=exp(p_H)*exp(y_H)/exp(m); //E32

end; 

%----------------------------------------------------------------
% 3. Steady State
%----------------------------------------------------------------

steady_state_model; 

// Computing the steady state and endogenous parameters
BETA=ppi_ss/R_ss;
ppi_ast_ss=ppi_ss/ppi_S_ss;
R_ast_ss=R_ss/ppi_S_ss;
R_W_ss=R_ast_ss;
ppi_H_ss=ppi_ss;
p_tilde_H_ss=1;
mc_H_ss=(EPSILON_H-1)/EPSILON_H;
y_H_ss=z_ss*h_ss;
f1H_ss=p_tilde_H_ss^(-EPSILON_H)/(1-THETA_H*BETA)*y_H_ss*mc_H_ss;
f2H_ss=f1H_ss;
p_F_ss=((1-OMEGA*p_H_ss^(1-ETA))/(1-OMEGA))^(1/(1-ETA));
rer_ss=p_F_ss;
tb_ss=s_tb_ss*p_H_ss*y_H_ss;
c_ss=p_H_ss*y_H_ss-tb_ss;
c_F_ss=(1-OMEGA)*rer_ss^(-ETA)*c_ss;
c_H_ss=OMEGA*p_H_ss^(-ETA)*c_ss;
c_H_ast_ss=y_H_ss-c_H_ss;
y_ast_ss=c_H_ast_ss*(p_H_ss/rer_ss)^(ETA_ast);
d_ast_ss=(tb_ss)*(rer_ss*(R_ast_ss/ppi_ast_ss-1))^(-1);
BAR_d=d_ast_ss;
w_ss=mc_H_ss*z_ss*p_H_ss;
la_ss=c_ss^(-SIGMA);
PSI=w_ss*la_ss*h_ss^(-VARPHI);
m_ss=p_H_ss*y_H_ss/vel_ss;
NU=m_ss/y_H_ss*R_ss^PHI_m;
d_M_ss=ppi_ss;
s_d_ast_ss=rer_ss*d_ast_ss/(p_H_ss*y_H_ss);

// Initial values for numerical solver
c=log(c_ss);
la=log(la_ss);
h=log(h_ss);
w=log(w_ss);
ppi=log(ppi_ss);
ppi_S=log(ppi_S_ss);
R=log(R_ss);
R_ast=log(R_ast_ss);
c_H=log(c_H_ss);
c_F=log(c_F_ss);
p_H=log(p_H_ss);
p_F=log(p_F_ss);
y_H=log(y_H_ss);
mc_H=log(mc_H_ss);
f1H=log(f1H_ss);
f2H=log(f2H_ss);
p_tilde_H=log(p_tilde_H_ss);
ppi_H=log(ppi_H_ss);
rer=log(rer_ss);
c_H_ast=log(c_H_ast_ss);
d_ast=d_ast_ss;
tb=tb_ss;
m=log(m_ss);
d_M=log(d_M_ss);
z=log(z_ss);
R_W=log(R_W_ss);
y_ast=log(y_ast_ss);
ppi_ast=log(ppi_ast_ss);
e_m=log(e_m_ss);
s_tb=s_tb_ss;
s_d_ast=s_d_ast_ss;
vel=log(vel_ss);

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

steady;
check;

shocks;
var u_z; stderr 0.01;
var u_R_W; stderr 0.01;
//var u_y_ast; stderr 0;
//var u_ppi_ast; stderr 0;
var u_e_m; stderr 0.01;
end;

stoch_simul(order=1, periods=0, irf=40, nograph);
save NKSOE_simple_ej1g1.mat