///// NK 2 country model  /////
//F. Gulcin Ozkan and D. Filiz Unsal (2012)//
 
var
C_HD 
C_MD 
C_D 
C_HF 
C_MF 
C_F 
p_HD 
p_MD 
p_D 
p_HF 
p_MF 
p_F 
mc_D 
mc_F 
p_XD 
p_XF 
I_D
I_HD 
I_MD 
I_F 
I_HF 
I_MF 
K_D 
K_F 
Q_D 
Q_F 
THETA_D 
THETA_F 
RHO_D 
RHO_F
om_D 
om_F 
z_D   
g_D
zph_D  
gph_D  
z_F  
g_F 
zph_F 
gph_F  
i_D  
i_F 
r_KD 
r_KF 
D_ED 
D_EF
C_ED 
C_EF 
C_EHD
C_EMD
C_EHF 
C_EMF
y_D 
y_HD 
y_XD 
y_F 
y_HF 
y_XF  
y_MD 
y_MF 
w_ED 
w_EF 
nw_D
nw_F
w_D  
w_F  
n_D 
n_F 
h_D 
h_F 
D_HD  
D_HF 
pi_D  
pi_F  
psi_Dt
r_D 
r_F 
PHI_D
PHI_F
s_D
s_F
;
 
varexo e_g_D e_g_F;
 
parameters  beta sigma gamma phi eta lamda delta omega psi_I psi_D psi_H psi_X psi_M epsilon_pi rho 
chi_D mu_D theta_D  
chi_F mu_F theta_F
alpha
alpha_F
xi S  S_F
PHI_D_ss 
PHI_F_ss 
D_HD_ss 
y_D_ss 
om_D_ss
om_F_ss
i_D_ss
i_F_ss
;

beta = 0.99;
sigma = 1;
gamma = 2;
phi = 2;
alpha = 0.315;
eta = 0.35;
lamda = 11;
delta = 0.025;
omega = 0.01;
alpha_F = 0.035;
psi_I = 12;
psi_D = 0.0075;
psi_H = 120;
psi_X = 120;
psi_M = 120;
epsilon_pi = 1.5;
rho = 0.5;
chi_D = 0.9;
mu_D = 0.2;
theta_D = 0.9873;
chi_F = 0.5;
mu_F = 0.12;
theta_F = 0.9966;
xi = 1 ;
PHI_D_ss = 0.02+7.9431e-06;
PHI_F_ss = 0.005-4.3497e-05;

/// MODEL EQUATIONS ///
model;
 
//DOMESTIC ECONOMY
// households
 
C_D = ((alpha)^(1/gamma)*C_HD^((gamma-1)/gamma) + (1-alpha)^(1/gamma)*C_MD^((gamma-1)/gamma))^(gamma/(gamma-1));
C_HD = (alpha)*((p_HD)/(p_D))^(-gamma)*C_D;
C_MD =  (1-alpha)*((p_MD)/(p_D))^(-gamma)*C_D;
p_D = ((alpha)*(p_HD)^(1-gamma) + (1-alpha)*(p_MD)^(1-gamma))^(1/(1-gamma));
h_D^phi/C_D^(-sigma) = w_D/p_D;
psi_Dt = psi_D/2*(exp((D_HD(+1)+D_ED)/(p_D*y_D) - (D_HD_ss + D_ED)/(p_D*y_D_ss)) - 1)^2;
C_D^(-sigma) = beta*(1+i_D)*(C_D(+1) )^(-sigma)*p_D/p_D(+1);
 
//production firms
y_D = n_D^(1-eta)*K_D^eta; 
w_D = ((1 - eta)*(1-omega)*y_D * mc_D)/h_D;
w_ED = ((1 - eta)*omega*y_D * mc_D);
r_D = (eta*y_D * mc_D)/K_D;
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))^(-phi)/(C_D )^(-phi) * p_D/p_D(+1);
 
//importing firms

p_MD =  lamda/(lamda - 1)*s_D*p_F - 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 = ((alpha)^(1/gamma)*(alpha*(p_HD/p_D)^(-gamma)*I_D)^((gamma-1)/gamma) + (1-alpha)^(1/gamma)*((1-alpha)*(p_MD/p_D)^(-gamma)*I_D)^((gamma - 1)/gamma))^(gamma/(gamma-1));
I_HD 	=	(alpha)*(p_HD/p_D)^(-gamma)*I_D;
I_MD 	=	(1-alpha)*(p_MD/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
 
C_EHD = (alpha)*(p_HD/p_D)^(-gamma)*C_ED;
C_EMD=	(1-alpha)*(p_MD/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(+1)/s_D);
r_KD*Q_D(-1)*K_D*g_D = (1+i_F(-1))*s_D*D_ED;
p_D*nw_D = theta_D*(r_KD*Q_D(-1)*K_D*z_D)+w_ED;
p_D*nw_D = theta_D*(r_KD*Q_D(-1)*K_D*(z_D+g_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*(z_D+g_D) - (1+i_F(-1))*s_D*D_ED);
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);
z_D = 1 - normcdf(((log(om_D)+0.5*S^2)/S) - S) - om_D_ss*(1-normcdf(((log(om_D)+0.5*S^2)/S)));
g_D = RHO_D*(om_D/RHO_D*(1-normcdf(((log(om_D/RHO_D)+0.5*S^2)/S))) + (1-mu_D)*normcdf(((log(om_D/RHO_D)+0.5*S^2)/S) - S));
zph_D = -(1 - normcdf(((log(om_D)+0.5*S^2)/S)));
gph_D = RHO_D*(1-normcdf(((log(om_D/RHO_D)+0.5*S^2)/S)) - mu_D/(sqrt(2*3.141592654)*S)*exp(-((log(om_D/RHO_D)+0.5*S^2)/S)^2/2)); 
om_F = om_D*RHO_D;
log(RHO_D) = rho*log(RHO_D(-1)) + xi*log(RHO_F) + e_g_D;

 //monetary policy
1+i_D = (1+i_D_ss)*(pi_D)^(epsilon_pi);
 
//general equilibrium
y_D = y_HD + y_XD;
y_HD = C_HD + C_EHD + I_HD + (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 + (1-z_D-g_D)*r_KD/p_D*Q_D(-1)*K_D);
y_MD = C_MD + C_EMD + I_MD + (1-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 + (1-z_D-g_D)*r_KD/p_D*Q_D(-1)*K_D);
n_D = h_D^(1-omega);
p_XD*y_XD - p_F*y_MD = (1+i_F(-1))*(D_HD*psi_Dt(-1)+D_ED) - (D_HD(+1) + D_ED(+1));
 
 

 
//FOREIGN ECONOMY 
 // households
 
C_F = ((alpha_F)^(1/gamma)*C_HF^((gamma-1)/gamma) + (1-alpha_F)^(1/gamma)*C_MF^((gamma-1)/gamma))^(gamma/(gamma-1));
C_HF = (alpha_F)*((p_HF)/(p_F))^(-gamma)*C_F;
C_MF = (1-alpha_F)*((p_MF)/(p_F))^(-gamma)*C_F;
p_F = ((alpha_F)*(p_HF)^(1-gamma) + (1-alpha_F)*(p_MF)^(1-gamma))^(1/(1-gamma));
h_F^phi/C_F^(-sigma) = w_F/p_F;
C_F^(-sigma) = beta*(1+i_F)*(C_F(+1))^(-sigma)*p_F/p_F(+1);

//production firms
y_F = n_F^(1-eta)*K_F^eta; 
w_F = ((1 - eta)*(1-omega)*y_F * mc_F)/h_F;
w_EF = ((1 - eta)*omega*y_F * mc_F);
r_F = (eta*y_F * mc_F)/K_F;
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))^(-phi)/(C_F )^(-phi) * p_F/p_F(+1);
 
//importing firms
s_F = 1/s_D; 
//Unfinished captial producing firms
 
I_F =  ((alpha_F)^(1/gamma)*I_HF^((gamma-1)/gamma) + (1-alpha_F)^(1/gamma)*I_MF^((gamma - 1)/gamma))^(gamma/(gamma-1));
I_HF = (alpha_F)*(p_HF/p_F)^(-gamma)*I_F;
I_MF = (1-alpha_F)*(p_MF/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_F(-1)*K_F*g_F = (1+i_F(-1))*s_F*D_EF;
C_EHF = (alpha_F)*(p_HF/p_F)^(-gamma)*C_EF;
C_EMF = (1-alpha_F)*(p_MF/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)*(s_F(+1)/s_F);
p_F*nw_F = theta_F*(r_KF*Q_F(-1)*K_F*z_F)+w_EF;
p_F*nw_F = theta_F*(r_KF*Q_F(-1)*K_F*(z_F+g_F)- (1 +i_D(-1))*s_F*D_EF) + w_EF;
p_F*C_EF = (1-theta_F)*(r_KF*Q_F(-1)*K_F*(z_F+g_F) - (1+i_D(-1))*s_F*D_EF);

z_F = 1 - normcdf(((log(om_F)+0.5*S_F^2)/S_F) - S_F) - om_F_ss*(1-normcdf(((log(om_F)+0.5*S_F^2)/S_F)));
g_F = RHO_F*(om_F/RHO_F*(1-normcdf(((log(om_F/RHO_F)+0.5*S_F^2)/S_F))) + (1-mu_F)*normcdf(((log(om_F/RHO_F)+0.5*S_F^2)/S_F) - S_F));
zph_F = -(1 - normcdf(((log(om_F)+0.5*S_F^2)/S_F)));
gph_F = RHO_F*(1-normcdf(((log(om_F/RHO_F)+0.5*S_F^2)/S_F)) - mu_F/(sqrt(2*3.141592654)*S_F)*exp(-((log(om_F/RHO_F)+0.5*S_F^2)/S_F)^2/2)); 
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_F(+1)/K_F(+1)-delta)^2);
log(RHO_F) = rho*log(RHO_F(-1)) + xi*log(RHO_D) + e_g_F;
om_D = om_F*RHO_F;

//monetary policy
1+i_F = (1+i_F_ss)*(pi_F)^(epsilon_pi);
 
//general equilibrium
y_F = y_HF + y_XF;
y_HF = C_HF + C_EHF + I_HF + (alpha_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 + (1-z_F-g_F)*r_KF/p_F*Q_F(-1)*K_F);
y_MF = C_MF + C_EMF + I_MF + (1-alpha_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 + (1-z_F-g_F)*r_KF/p_F*Q_F(-1)*K_F);
n_F = h_F^(1-omega);
p_XF*y_XF - p_D*y_MF = (1+i_D(-1))*(D_HF*psi_Dt(-1)+D_EF) - (D_HF(+1) + D_EF(+1));
 
end;
 
initval;
C_HD 	=	0;
C_MD 	=	0;
C_D 	=	0;
C_HF 	=	0;
C_MF 	=	0;
C_F 	=	0;
p_HD 	=	0;
p_MD 	=	0;
p_D 	=	0;
p_HF 	=	0;
p_MF 	=	0;
p_F 	=	0;
mc_D 	=	0;
mc_F 	=	0;
p_XD 	=	0;
p_XF 	=	0;
I_D 	=	0;
I_HD    =   0;
I_MD    =   0;
I_F 	=	0;
I_HF 	=	0;
I_MF 	=	0;
K_D 	=	0;
K_F 	=	0;
Q_D 	=	0;
Q_F 	=	0;
THETA_D 	=	0;
THETA_F 	=	0;
RHO_D 	=	0;
RHO_F	=	0;
om_D 	=	0;
om_F 	=	0;
z_D   	=	0;
g_D	=	0;
zph_D  	=	0;
gph_D  	=	0;
z_F  	=	0;
g_F 	=	0;
zph_F 	=	0;
gph_F  	=	0;
i_D  	=	0;
i_F 	=	0;
r_KD 	=	0;
r_KF 	=	0;
D_ED 	=	0;
D_EF	=	0;
C_ED 	=	0;
C_EF 	=	0;
C_EHD	=	0;
C_EMD	=	0;
C_EHF 	=	0;
C_EMF	=	0;
y_D 	=	0;
y_HD 	=	0;
y_XD 	=	0;
y_F 	=	0;
y_HF 	=	0;
y_XF  	=	0;
y_MD 	=	0;
y_MF 	=	0;
w_ED 	=	0;
w_EF 	=	0;
nw_D	=	0;
nw_F	=	0;
w_D  	=	0;
w_F  	=	0;
n_D 	=	0;
n_F 	=	0;
h_D 	=	0;
h_F 	=	0;
D_HD  	=	0;
D_HF 	=	0;
pi_D  	=	0;
pi_F  	=	0;
psi_Dt	=	0;
r_D 	=	0;
r_F 	=	0;
PHI_D	=	0;
PHI_F	=	0;
s_D     =   0;
s_F     =   0;

end;
steady(solve_algo=3); check;
 
 
shocks;
    var e_g_D;
    stderr 0.01; 
    var e_g_F;
    stderr 0.01;
end;
 
stoch_simul(order=2,irf=4000,nograph);