% Endogenous variables 
var Cih Cif Csh Csf kihh kihf kiff kifh ksh ksf Bih Bif Bsh Bsf q1 q2 Rh Rf Wh Wf Rkh Rkf Ah Af Zh Zf Yh Yf Jh Jf BETAi BETAs;

% Exogenous variables
varexo nu1 nu2 nu3 nu4;

parameters Epsilon ETA n RHO SIGMA OMEGA DELTA r ZETAi ZETAs sigmanu1 sigmanu2 sigmanu3 sigmanu4; 

%Parameters values
Epsilon = 0.36;     %Capital share in final goods production
ETA = 0.022;        %Discount function
n= 0.5;            %Proportion of investors
RHO = 0.9;         % Productivity shock persistence for investors in home country producing final good
SIGMA = 2;          % coefficient of relative risk aversion
ZETAi = 0.975;      %Investors discount function
ZETAs = 0.9825;     %Savers discountfunction
OMEGA = 0.1;        %Capital share in home production
DELTA = 0.0226;     %capital depreciation rate
r = 0.04;           %world interest rate
sigmanu1 = 0.02;
sigmanu2 = 0.02;
sigmanu3 = 0.0129;
sigmanu4 = 0.0129;
corr_nu1_nu2 =0.00;
corr_nu3_nu4 =0.00;

% steady state variables
BETAi_ss = 1/(1+r);
BETAs_ss = 1/(1+r);
Rh_ss = (1/BETAs_ss);
Rf_ss = (1/BETAs_ss);
Rkh_ss = (1/BETAs_ss);
Rkf_ss = (1/BETAs_ss);
kihh_ss = (1/(2*n))*(Rkh_ss/Epsilon)^(1/(Epsilon-1));
kihf_ss = (1/(2*n))*(Rkh_ss/Epsilon)^(1/(Epsilon-1));
kiff_ss = (1/(2*n))*(Rkh_ss/Epsilon)^(1/(Epsilon-1));
kifh_ss = (1/(2*n))*(Rkh_ss/Epsilon)^(1/(Epsilon-1));
kih_ss = n*(kihh_ss+kifh_ss);
kif_ss = n*(kiff_ss+kihf_ss);
ksh_ss = (1-kih_ss)/(1-n);
ksf_ss = (1-kih_ss)/(1-n);
Yh_ss =(kih_ss)^Epsilon;
Yf_ss =(kif_ss)^Epsilon;
Jh_ss = (ksh_ss)^OMEGA;
Jf_ss = (ksf_ss)^OMEGA;
Cih_ss = Yh_ss ;
Cif_ss = Yf_ss; 
Csh_ss = Jh_ss;
Csf_ss = Jf_ss; 
Wh_ss = (1-Epsilon)*(kih_ss)^Epsilon;
Wf_ss =  (1-Epsilon)*(kif_ss)^Epsilon;
Bih_ss = 1;
Bif_ss = 1;
Bsh_ss = 1;
Bsf_ss = 1;
q1_ss = 2;
q2_ss = 2;
Ah_ss = 1;
Af_ss = 1;
Zh_ss = 1;
Zf_ss = 1;


model;
%  discount factor 1
exp(BETAi)= ZETAi*(1+ exp(Cih))^(-ETA);
%discount factor 2
exp(BETAs)= ZETAs*(1+ exp(Csh))^(-ETA);
%equation 2  /budget constraint investment home country 3  FIRST SUBSCRIPT IS WHERE CONSUMED/USED  SECOND SUBSCRIPT IS WHERE PRODUCED
exp(Cih)+(exp(q1)*exp(kihh))+(exp(q2)*exp(kihf))= exp(Wh)+(exp(q1)+exp(Rkh))*exp(kihh(-1))+(exp(q2)+exp(Rkf))*exp(kihf(-1))+ exp(Bih)-exp(Rh(-1))*exp(Bih(-1));

%equation 2' /budget constraint investment foreign country 4
exp(Cif)+ (exp(q2)*exp(kiff))+(exp(q1)*exp(kifh))= exp(Wf)+(exp(q2)+exp(Rkf))*exp(kiff(-1))+(exp(q1)+exp(Rkh))*exp(kifh(-1))+exp(Bif)-exp(Rf(-1))*exp(Bif(-1));

%equation 4 /FOC capital 1 home country  5 wrt kihh
exp(Cih)^(-SIGMA)= exp(BETAi) *exp(Cih(+1))^(-SIGMA)*(exp(q1(+1))+exp(Rkh(+1)))/exp(q1);

%equation 4' / FOC capital 1 foreign country 6 wrt kifh
exp(Cif)^(-SIGMA) = exp(BETAi) * exp(Cif(+1))^(-SIGMA)*(exp(q1(+1))+exp(Rkh(+1)))/exp(q1);

%equation 5 / FOC capital 2 home country 7 wrt kihf
exp(Cih)^(-SIGMA) = exp(BETAi) * exp(Cih(+1))^(-SIGMA)*(exp(q2(+1))+exp(Rkf(+1)))/exp(q2);

%equation 5'  / FOC capital 2 foreign country 8 wrt kiff
exp(Cif)^(-SIGMA) = exp(BETAi) * exp(Cif(+1))^(-SIGMA)*(exp(q2(+1))+exp(Rkf(+1)))/exp(q2);

%equation 6    / FOC bond home country 9
exp(Cih)^(-SIGMA) = exp(BETAi) * exp(Cih(+1))^(-SIGMA)*exp(Rh(+1));

%equation 6'   / FOC bond foreign country 10
 exp(Cif)^(-SIGMA) = exp(BETAi) * exp(Cif(+1))^(-SIGMA)*exp(Rf(+1));

%equation9  / budget constraint savers home country 11
exp(Csh)+ (exp(q1)*exp(ksh))= exp(Wh)+ exp(q1)*exp(ksh(-1))+ exp(Zh)*exp(ksh(-1))^(OMEGA)+ exp(Bsh)- (exp(Rh(-1))*exp(Bsh(-1)));

%equation 9'   / budget constraint saves foreign country 12
exp(Csf)+ exp(q2)*exp(ksf)= exp(Wf)+ exp(q2)*exp(ksf(-1))+exp(Zf)*exp(ksf(-1))^(OMEGA)+ exp(Bsf)- exp(Rf(-1))*exp(Bsf(-1));

%equation 11   / FOC capital home country 13
exp(Csh)^(-SIGMA) = exp(BETAs)*exp(Csh(+1))^(-SIGMA)*(exp(q1(+1))+ OMEGA*exp(Zh(+1))*exp(ksh)^(OMEGA-1))/exp(q1);


%equation 11'  / FOC capital foreign country 14
exp(Csf)^(-SIGMA) = exp(BETAs)*exp(Csf(+1))^(-SIGMA)*(exp(q2(+1))+ OMEGA*exp(Zf(+1))*exp(ksf)^(OMEGA-1))/exp(q2);

% Equation 12   / FOC bond home country (for savers)   15
exp(Csh)^(-SIGMA) = exp(BETAs)*exp(Csh(+1))^(-SIGMA)*exp(Rh);

% Equation 12' / FOC bond foreign country (for savers)  16
exp(Csf)^(-SIGMA) = exp(BETAs)*exp(Csf(+1))^(-SIGMA)* exp(Rf);

%Equation 15  / FOC with respect to L home country    17
exp(Wh)= exp(Ah)*(1-Epsilon)*(n*(exp(kihh(-1))+exp(kifh(-1))))^(Epsilon);

%Equation 15' / FOC with respect to L foreign country   18
exp(Wf)= exp(Af)*(1-Epsilon)*(n*(exp(kiff(-1))+exp(kihf(-1))))^(Epsilon);

%Equation 16  / FOC with respect to K home country   19
exp(Rkh)= exp(Ah)*(1-Epsilon)*(n*(exp(kihh(-1))+exp(kifh(-1))))^(Epsilon);

%Equation 16'   / FOC with respect to K foreign country  20
exp(Rkf)=exp(Af)*(1-Epsilon)*(n*(exp(kiff(-1))+exp(kihf(-1))))^(Epsilon);

%equation 17    / AR(1) home investors  21
Ah= RHO*Ah(-1)+ nu1;
  
% equation 18  / AR(1) foreign investors  22
Af= RHO*Af(-1)+ nu2;

%The stochastic process for home goods productivity in the home and foreign countries
% equation 19   / AR(1) Home savers 23
Zh= RHO*Zh(-1)+ nu3;

% equation20  / AR(1) foreign savers  24
Zf= RHO*Zf(-1)+ nu4;

%Equation 21    / bond market home country  25
%n*exp(Bih)+(1-n)*exp(Bsh)=0;

%Equation 22 / fixed asset equilibrium home country  26
n*exp(kihh)+n*exp(kifh)+(1-n)*exp(ksh)=1;

%Equation 21'   / bond market foreign country  27
n*exp(Bif)+(1-n)*exp(Bsf)=0;

%Equation 22'   / fixed asset equilibrium foreign country   28
n*exp(kiff) + n*exp(kihf)+ (1-n)*exp(ksf)=1;

% Production functions for home and foreign investors
% equation 13    / production function home investors  29
exp(Yh)= exp(Ah)*(n*(exp(kihh(-1))+exp(kifh(-1))))^Epsilon;
%Kh =kihh+kifh
% equation 13'    / production function foreign investors    30
exp(Yf)= exp(Af)*(n*(exp(kiff(-1))+exp(kihf(-1))))^Epsilon;
%Kf = kiff+kihf
% Production functions for home and foreign savers

%equation 14     / production function home savers 31
exp(Jh)= exp(Zh)*(exp(ksh(-1)))^(OMEGA);

% equation 14'    / production function foreign savers 32
exp(Jf)= exp(Zf)*(exp(ksf(-1)))^(OMEGA);

% equation 23    / world market clearing condition   33
n*exp(Cih)+ n*exp(Cif)+(1-n)*(exp(Csh)+exp(Csf))= exp(Yh)+exp(Yf)+(1-n)*exp(Jh)+(1-n)*exp(Jf);

end;

initval;
BETAi = log(BETAi_ss);
BETAs = log(BETAs_ss);
Rh    = log(Rh_ss);
Rf    = log(Rf_ss);
Rkh   = log(Rkh_ss);
Rkf   = log(Rkf_ss);
kihh  = log(kihh_ss);
kihf  = log(kihf_ss);
kiff  = log(kiff_ss);
kifh  = log(kifh_ss);
ksh   = log(ksh_ss);
ksf   = log(ksf_ss);
Yh    = log(Yh_ss);
Yf    = log(Yf_ss);
Jh    = log(Jh_ss);
Jf    = log(Jf_ss);
Cih   = log(Cih_ss) ;
Cif   = log(Cif_ss); 
Csh   = log(Csh_ss);
Csf   = log(Csf_ss);
Wh    = log(Wh_ss);
Wf    = log(Wf_ss);
Bih   = log(Bih_ss);
Bif   = log(Bif_ss);
Bsh   = log(Bsh_ss);
Bsf   = log(Bsf_ss);
q1    = log(q1_ss);
q2    = log(q2_ss);
Ah    = log(Ah_ss);
Af    = log(Af_ss);
Zh    = log(Zh_ss);
Zf    = log(Zf_ss);
nu1   =0;
nu2   =0;
nu3   =0;
nu4   =0;
end;
shocks;
var nu1 = sigmanu1^2;
var nu2 =sigmanu2^2;
var nu3 =sigmanu3^2;
var nu4 = sigmanu4^2;
end;
resid(1);
steady(solve_algo = 4,maxit=1000); 

stoch_simul(hp_filter=1600,order=1,irf=80);