///// NK 2 country model  /////
//F. Gulcin Ozkan and D. Filiz Unsal (2012)//

var B_D C_D C_HD C_MD THETA_D D_ED D_HD h_D i_D I_D I_HD I_MD K_D mc_D nw_D p_D p_HD p_MD p_XD
Q_D r_D r_KD s_D w_D w_ED y_D y_HD y_XD psi_Dt z_D REX pi
C_F C_HF C_MF THETA_F D_EF h_F i_F I_F I_HF I_MF K_F mc_F nw_F p_F p_HF p_M p_XF 
Q_F r_F r_KF s_F w_F w_EF y_F y_HF g_D g_F om om_F n_D
n_F a_D C_EMF p_MF V_D V_F y_MD C_ED C_EHD C_EMD zph_D zph_F gph_D gph_F
z1 z2 ;

vaxero e

parameters n beta sigma gamma phi v eta lamda delta omega psi_I psi_D psi_H psi_X psi_M epsilon_pi epsilon_y rho 
chi_D phi_D mu_D theta_D
chi_F phi_F mu_F theta_F
alpha
xi S M S_F M_F
;
n = 0.1;
beta = 0.99;
sigma = 1;
gamma = 1;
phi = 2;
v = 0.35;
eta = 0.35;
lamda = 11;
delta = 0.025;
omega = 0.01;
psi_I = 12;
psi_D = 0.0075;
psi_H = 120;
psi_X = 120;
psi_M = 120;
epsilon_pi = 1.5;
epsilon_y = 0.5;
rho = 0.5;
chi_D = 0.3;
phi_D = 0.02;
mu_D = 0.2;
theta_D = 0.9933;
chi_F = 0.5;
phi_F = 0.005;
mu_F = 0.12;
theta_F = 0.9966;

alpha = (1-n)*v;
xi = 0.5 ;
S    = 0.207;
M    = -0.5*S^2;
S_F = 0.207;
M_F = -0.5*S_F^2
RHO_D = rho_D*RHO_D(-1) + xi*RHO_F+ e;

/// MODEL EQUATIONS ///
model(linear);

//DOMESTIC ECONOMY
// households
C_D = ((1-alpha)^(1/gamma)*C_HD^((gamma-1)/gamma) + alpha^(1/gamma)*C_MD^((gamma-1)/gamma))^(gamma/(gamma-1));
C_HD = (1 - alpha)*((p_HD)/(p_D))^(-gamma)*C_D;
C_MD =  alpha*((p_MD)/(p_D))^(-gamma)*C_D;
p_D = ((1-alpha)*(p_HD)^(1-gamma) + alpha*(p_MD)^(1-gamma))^(1/(1-gamma));
h_D^phi = w_D;
(C_D - h_D^(1+phi)/(1+phi))^(-sigma) = beta*(1+i_D)*(C_D(+1) - h_D(+1)^(1+phi)/(1+phi))^(-sigma)*p_D/p_D(+1);
(C_D - h_D^(1+phi)/(1+phi))^(-sigma) = beta*(1+i_F)*psi_Dt*(C_D(+1) - h_D(+1)^(1+phi)/(1+phi))^(-sigma)*p_D/p_D(+1)*s_D(+1)/s_D;
REX = s_D*p_F/p_D;
p_D*C_D + (1+i_D(-1))*B_D + (1+i_F(-1))*psi_Dt(-1)*s_D*D_HD = w_D*h_D + B_D(+1) + s_D*D_HD(+1) + pi;

//production firms


w_D = ((1 - eta)*(1-omega)*y_D * mc_D)/n_D;
w_ED = ((1 - eta)*omega*y_D * mc_D);
r_D = (eta*y_D * mc_D)/K_D;
mc_D =( r_D^eta*w_D^(1-eta))/(a_D*eta^eta*(1-eta)^(1-eta));
p_HD = lamda/(lamda - 1)*mc_D - psi_H/(lamda-1)*p_D/y_HD*p_HD/p_HD(-1)*(p_HD/p_HD(-1)-1) + psi_H/(lamda-1)* (THETA_D*p_D(+1)/y_HD*p_HD(+1)/p_HD*(p_HD(+1)/p_HD - 1));
s_D*p_XD = lamda/(lamda - 1)*mc_D - psi_X/(lamda-1)*p_D/y_XD*p_XD/p_XD(-1)*(p_XD/p_XD(-1)-1) + psi_X/(lamda-1) * (THETA_D*p_D(+1)/y_XD*p_XD(+1)/p_XD*(p_XD(+1)/p_XD - 1));
THETA_D = beta*(C_D(+1) - chi_D/(1 + phi)*h_D(+1)^(1+phi))^(-phi)/(C_D - chi_D/(1 + phi)*h_D^(1+phi))^(-phi) * p_D/p_D(+1);
y_XD = C_MF + C_EMF + I_MF + alpha*(p_MF/p_F)^(-gamma)*(psi_H/2*(p_HF/p_HF(-1)-1)^2 + psi_X/2*(p_XF/p_XF(-1)-1)^2 + psi_M/2*(p_MF/p_MF(-1)-1)^2 + V_F*r_KF/p_F*Q_F(-1)*K_F);

//importing firms

p_MD =  lamda/(lamda - 1)*s_D*p_MF - psi_M/(lamda-1)*p_D/y_MD*p_MD/p_MD(-1)*(p_MD/p_MD(-1)-1) + psi_M/(lamda-1)* (THETA_D*p_D(+1)/y_MD*p_MD(+1)/p_MD*(p_MD(+1)/p_MD - 1));

//Unfinished captial producing firms

I_D = ((1-alpha)^(1/gamma)*I_HD^((gamma-1)/gamma) + alpha^(1/gamma)*I_MD^((gamma - 1)/gamma))^(gamma/(gamma-1));
I_HD = (1 - alpha)*(p_HD/p_D)^(-gamma)*I_D;
I_MD = alpha*(p_HD/p_D)^(-gamma)*I_D;
K_D = (I_D/K_D(-1) - psi_I/2*(I_D/K_D(-1) - delta)^2) * K_D(-1) + (1 - delta)*K_D(-1);
Q_D = ((1 - psi_I*(I_D/K_D - delta))^(-1))*p_D;

//entrepreneurs

r_KD*Q_D(-1)*K_D/s_D*g_D = (1+i_F(-1))*D_ED;
p_D*nw_D = theta_D*(r_KD*Q_D(-1)*K_D*(1 - V_D)- (1 - i_F(-1))*s_D*D_ED) + w_ED;
p_D*C_ED = (1-theta_D)*(r_KD*Q_D(-1)*K_D*(1-V_D) - (1+i_F(-1))*s_D*D_ED);
C_EHD = (1-alpha)*(p_HD/p_D)^(-gamma)*C_ED;
C_EMD = alpha*(p_HD/p_D)^(-gamma)*C_ED;
r_KD(+1) = (1+i_F)*(1+phi_D(+1));
1+ phi_D = zph_D/(g_D*zph_D - z_D*gph_D)*s_D/s_D(-1);
r_KD(+1) = r_D(+1)/Q_D + Q_D(+1)/Q_D * ((1-delta)+psi_I*(I_D(+1)/K_D(+1)-delta)*I_D(+1)/K_D(+1)-psi_I/2*(I_D(+1)/K_D(+1)-delta)^2);
 
z1 = (log(om)+0.5*S^2)/S;
z2 = (log(om/RHO_D)+0.5*S^2)/S;
z_D = 1 - (z1 - S) - om_ss*(1-normcdf(z1));
g_D = RHO_D*(om_ss/RHO_D*(1-normcdf(z2)) + (1-mu_D)*normcdf(z2 - S));
zph_D = -(1 - normcdf(z1));
gph_D = RHO_D*(1-normcdf(z2) - mu_D/(sqrt(2*3.14159)*S)*exp(-z2^2/2)); 
V_D = 1- z_D - g_D;


//monetary policy
1+i_D = (1+i_D_ss)*(pi_D)^(epsilon_pi)*(y_D/Y_D_ss)^(epsilon_y);
s_D/s_D(-1) = 1;

//general equilibrium
y_D = y_HD + y_XD;
y_HD = C_HD + C_EHD + I_HD+(1-alpha)*(p_HD/p_D)^(-gamma)*(psi_H/2*(p_HD/p_HD(-1)-1)^2 + psi_X/2*(p_XD/p_XD(-1)-1)^2 + psi_M/2*(p_MD/p_MD(-1)-1)^2 + V_D*r_KD/p_D*Q_D(-1)*K_D);
y_XF = C_MD + C_EMD + I_MD + alpha*(p_MD/p_D)^(-gamma)*(psi_H/2*(p_HD/p_HD(-1)-1)^2 + psi_X/2*(p_XD/p_XD(-1)-1)^2 + psi_M/2*(p_MD/p_MD(-1)-1)^2 + V_D*r_KD/p_D*Q_D(-1)*K_D);
n_D = h_D^(1-omega);
s_D*p_XD*y_XD - p_MF*y_MD = s_D*(1+i_F(-1))*(D_HD*psi_Dt(-1)+D_ED) - s_D*(D_HD(+1) + D_ED(+1));



//FOREIGN ECONOMY 

// households
 
C_HF = (1 - alpha_F)*((p_HF)/(p_F))^(-gamma)*C_F;
C_MF =  alpha_F*((p_MF)/(p_F))^(-gamma)*C_F;
p_F = ((1-alpha_F)*(p_HF)^(1-gamma) + alpha_F*(p_MF)^(1-gamma))^(1/(1-gamma));
h_F^phi = w_F;
(C_F - h_F^(1+phi)/(1+phi))^(-sigma) = beta*(1+i_F)*(C_F(+1) - h_F(+1)^(1+phi)/(1+phi))^(-sigma)*p_F/p_F(+1);
 
//production firms
n_F = h_F^(1-omega);
w_F = ((1 - eta)*(1-omega)*y_F * mc_F)/n_F;
w_EF = ((1 - eta)*omega*y_F * mc_F);
r_F = (eta*y_F * mc_F)/k_F;
mc_F =( r_F^eta*w_F^(1-eta))/(a_F*eta^eta*(1-eta)^(1-eta));
p_HF = lamda/(lamda - 1)*mc_F - psi_H/(lamda-1)*p_F/y_HF*p_HF/p_HF(-1)*(p_HF/p_HF(-1)-1) + psi_H/(lamda-1)* (THETA_F*p_F(+1)/y_HF*p_HF(+1)/p_HF*(p_HF(+1)/p_HF - 1));
s_F*p_XF = lamda/(lamda - 1)*mc_F - psi_X/(lamda-1)*p_F/y_XF*p_XF/p_XF(-1)*(p_XF/p_XF(-1)-1) + psi_X/(lamda-1) * (THETA_F*p_F(+1)/y_XF*p_XF(+1)/p_XF*(p_XF(+1)/p_XF - 1));
THETA_F = beta*(C_F(+1) - chi_F/(1 + phi)*h_F(+1)^(1+phi))^(-phi)/(C_F - chi_F/(1 + phi)*h_F^(1+phi))^(-phi) * p_F/p_F(+1);
 
//importing firms
 
p_MD =  lamda/(lamda - 1)*p_XD - psi_M/(lamda-1)*p_D/y_MD*p_MD/p_MD(-1)*(p_MD/p_MD(-1)-1) + psi_D/(lamda-1)* (THETA_F*p_D(+1)/y_MD*p_MD(+1)/p_MD*(p_MD(+1)/p_MD - 1));
 
//Unfinished captial producing firms
 
I_HF = (1 – alpha_F_F)*(p_HF/p_F)^(-gamma)*I_F;
I_MF = alpha_F_F*(p_HF/p_F)^(-gamma)*I_F;
K_F(+1) = (I_F/K_F - psi_I/2*(I_F/K_F - delta)^2) * K_F + (1 - delta)*K_F;
Q_F/p_F = (1 - psi_I*(I_F/K_F - delta))^(-1);
 
//entrepreneurs
 
 
r_KF*Q_D(-1)*K_F*g_F = (1+i_F(-1))*D_EF;
p_F*nw_F = theta_F*(r_KF*Q_F(-1)*K_F*(1 - V_F)- (1 - i_F(-1))*s_F*D_F) + w_EF;
p_F*C_EF = (1-theta_F)*(r_KF*Q_F(-1)*K_F*(1-V_F) - (1+i_F(-1))*s_F*D_F);
C_EHF = (1-alpha_F)*(p_HF/p_F)^(-gamma)*C_EF;
C_EMF = alpha_F*(p_HF/p_F)^(-gamma)*C_EF;
r_KF(+1) = (1+i_F)*(1+phi_F(+1));
1+ phi_F = zph_F/(g_F*zph_F - z_F*gph_F);
z1_F = (log(om_F)+0.5*S^2)/S;
z2_F = (log(om/RHO_D)+0.5*S^2)/S;
z_F = 1 - (z1_F - S_F) - om_F_ss*(1-normcdf(z1_F));
g_F = RHO_F*(om_F_ss/RHO_F*(1-normcdf(z2_F)) + (1-mu_D)*normcdf(z2_F - S_F));
zph_F = -(1 - normcdf(z1_F));
gph_F = RHO_F*(1-normcdf(z2_F) - mu_D/(sqrt(2*3.14159)*S_F)*exp(-z2_F^2/2)); 
V_F = 1- z_F - g_F;
r_KF(+1) = R_F(+1)/Q_F + Q_F(+1)/Q_F * ((1-delta)+psi_I*(I_F(+1)/K_F(+1)-delta)*I_F(+1)/K_F(+1)-psi_I/2*(I_D(+1)/K_F(+1)-delta)^2);

   
//monetary policy
1+i_F = (1+i_F_ss)*(pi_F)^(epsilon_pi)*(y_F/Y_F_ss)^(epsilon_y);
 
//general equilibrium
s_F = 1/s_D;
y_F = y_HF + y_XF;
y_HF = C_HF + C_EHF + I_HF+(1-alpha_F_F)*(p_HF/p_F)^(-gamma)*(psi_H/2*(p_HF/p_HF(-1)-1)^2 + psi_X/2*(p_XF/p_XF(-1)-1)^2 + psi_M/2*(p_MF/p_MF(-1)-1)^2 + V_F*r_KF/p_F*Q_F(-1)*K_F);
y_XF = C_MF + C_EMF + I_MF + alpha_F_F*(p_MF/p_F)^(-gamma)*(psi_H/2*(p_HF/p_HF(-1)-1)^2 + psi_X/2*(p_XF/p_XF(-1)-1)^2 + psi_M/2*(p_MF/p_MF(-1)-1)^2 + V_F*r_KF/p_F*Q_F(-1)*K_F);

end;

steady(solve_algo=3); check;


shocks;
    var e;
    stderr 0.01; 
end;

stoch_simul(order=2,irf=40,nograph);