close all;

var yh yf ch cf gh gf kh kf ih if th tf ah af hh hf qh qf;
varexo e_ah e_af e_th e_tf;

parameters beta alpha delta eta Lambda phi psi CoY IoY GoY TH TF rho_ah rho_af rho_th rho_tf;

beta=0.99;
delta = 0.025;
phi = 0.25;
psi = 0.5;
alpha = 0.34;
eta = 1/3;
Lambda = (1+0)*(1+0);
GoY = 0.2;
GoK = GoY*(1/alpha*(1/beta - (1-delta)));
IoK = Lambda - (1-delta);
IoY = IoK*(1/(1/alpha*(1/beta - (1-delta))));
CoK = 1/alpha*(1/beta - (1-delta)) - IoK - GoY*(1/alpha*(1/beta - (1-delta)));
CoY = CoK*(1/(1/alpha*(1/beta - (1-delta))));
TH = GoY;
TF = TH;
rho_ah = 0.9;
rho_af = 0.9;
rho_th = 0.9;
rho_tf = 0.9;

model(linear);
    % Euler home and foreign
	ch(+1)-ch=(beta*alpha)*(yh(+1)-kh + (1-delta));
    cf(+1)-cf=beta*alpha*(yf(+1)-kf + (1-delta));
	% Hours
	yh - hh - ch = eta*hh;
    yf - hf - cf = eta*hf;
	% Capital law of motion home and foreign
	Lambda*kh = delta*ih + (1-delta)*kh(-1);
    Lambda*kf = delta*if + (1-delta)*kf(-1);
	% Production functions
	yh = ah + alpha*kh + (1-alpha)*hh;
    yf = af + alpha*kf + (1-alpha)*hf;
    % Resources constraint
	psi*(yh-CoY*ch-IoY*ih-GoY*gh) = (1-psi)*(yf-CoY*cf-IoY*if-GoY*gf);
    % Tobin's Q
    qh = phi*(ih-kh);
    qf = phi*(if-kf);
    % Governments
    gh = 1/TH*th + yh;
    gf = 1/TF*tf + yf;
    % Exogenous shocks
	ah = rho_ah*ah(-1) + e_ah;
	af = rho_af*af(-1) + e_af;
	th = rho_th*th(-1) + e_th;
	tf = rho_tf*tf(-1) + e_tf;
    % Perfect risk sharing
    -ch + ch(+1) = -cf + cf(+1);
end;

initval;
yh = 0;
yf = 0;
ch = 0;
cf = 0;
ih = 0;
if = 0;
gh = 0;
gf = 0;
ah = 0;
af = 0;
th = 0;
tf = 0;
qh = 0;
qf = 0;
kh = 0;
kf = 0;
end;

steady;
check;
resid(1);


shocks;
var e_af; stderr 1;
end;


stoch_simul(order = 1,irf=30) yh yf ch cf kh kf;