% DYNARE CODE
%
%

%----------------------------------------------------------------
% 0. Let's start
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var rpn c ph ct cn ch cf cstar cnstar ctstar chstar cfstar pt pf pn phstar pfstar ptstar p pstar pnstar yn ynstar yt ytstar tot rer z bb bbb;
varexo epsilont epsilontstar epsilonn epsilonnstar e_z;

parameters at ah ahstar phi ytbar rho rhon rhot sigmat rho_z;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

at = 0.2;
ah  = 0.85;
ahstar = 0.15;

phi = 1.5;
ytbar = 0.5;

rho = 1.5;

rhon = 0.99;
rhot = 0.99;
sigma = 1.5;
sigmat = 0.0075;

rho_z=0.8;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 

ch = ah*(ph/pt)^(1-phi)*yt(-1);
cf = (1-ah)*(pf/pt)^(-phi)*(ph/pt)*yt(-1);
cn = (1-at)*(pn/p)^(-rho)*(yt(-1)*ph+yn(-1)*pn)/p;
chstar = ahstar*(phstar/ptstar)^(-phi)*(ytstar(-1)*pf/ptstar);
cfstar = (1-ahstar)*(pfstar/ptstar)^(1-phi)*ytstar(-1);
cnstar = (1-at)*(pnstar/pstar)^(-rho)*(ytstar(-1)*pf+ynstar(-1)*pnstar)/pstar;

c=(at^(1/rho)*ct^(1-1/rho)+(1-at)^(1/rho)*cn^(1-1/rho))^(rho/(rho-1));
ct=(ah^(1/phi)*ch^(1-1/phi)+(1-ah)^(1/phi)*cf^(1-1/phi))^(phi/(phi-1));
cstar=(at^(1/rho)*ctstar^(1-1/rho)+(1-at)^(1/rho)*cnstar^(1-1/rho))^(rho/(rho-1));
ctstar=((1-ah)^(1/phi)*chstar^(1-1/phi)+(ah)^(1/phi)*cfstar^(1-1/phi))^(phi/(phi-1));

p = (at*pt^(1-rho)+(1-at)*pn^(1-rho))^(1/(1-rho));
pt = (ah*ph^(1-phi)+(1-ah)*pf^(1-phi))^(1/(1-phi));
pstar = (at*ptstar^(1-rho)+(1-at)*pnstar^(1-rho))^(1/(1-rho));
ptstar = (ah*phstar^(1-phi)+(1-ah)*pfstar^(1-phi))^(1/(1-phi));

ph =1;

yn(-1)= cn;
ynstar(-1)= cnstar;

ph = phstar;
pf = pfstar;

tot = pf/phstar;
rer = pstar/p;
rpn = pn/pt;

% Bond Economy condition

((c/c(-1))^(-rho))*(pt(-1)/pt)=((cstar/cstar(-1))^(-rho))*(ptstar(-1)/ptstar);

% Endowments

yt = (1-rhot)*ytbar + rhot*yt(-1) + sigmat*epsilont +z;
ytstar = (1-rhot)*ytbar + rhot*ytstar(-1) + (sigmat)*(epsilontstar);
yn = (1-rhon)*ytbar + rhon*yn(-1) + sigmat*epsilonn;
ynstar = (1-rhon)*ytbar + rhon*ynstar(-1) + (sigmat)*(epsilonnstar);

% News shocks

z=bb(-1);
bb=bbb(-1);
bbb=rho_z*bbb(-1)+e_z;

end;

%----------------------------------------------------------------
% 4. Initial Values 
%----------------------------------------------------------------

initval;

rpn =1;
c =0.783045; 
ct = 0.335853;
cn = 0.5; 
ch =0.285475;
cf = 0.0503779;
cstar =0.783045;
cnstar =0.5;
ctstar = 0.335853; 
chstar =0.0503779;
cfstar =0.285475;
p =1.66299;
ph =1;
pf =1;
pt = 1;
pn = 1.93268;
pstar =1.66299;
phstar =1; 
pfstar =1;
ptstar = 1;
pnstar = 1.93268;

yt =0.5;
ytstar =0.5;
yn = 0.5;
ynstar = 0.5;

tot =1;
rer =1;
epsilont = 0;
epsilontstar = 0;
epsilonn = 0;
epsilonnstar = 0;
z=0;
bb=0;
bbb=0;

end;

%----------------------------------------------------------------
% 5. Shocks 
%----------------------------------------------------------------

shocks;
var epsilont = 1;
var epsilontstar = 1;
var epsilonn = 1;
var epsilonnstar = 1;
var e_z = 1;
end;

%----------------------------------------------------------------
% 6. Solve the model 
%----------------------------------------------------------------

steady;

stoch_simul(hp_filter = 1600, order = 1);

%----------------------------------------------------------------
% 7. Some Results
%----------------------------------------------------------------

%statistic1 = 100*sqrt(diag(oo_.var(1:6,1:6)))./oo_.mean(1:6);
%dyntable('Relative standard deviations in %',strvcat('VARIABLE','REL. S.D.'),M_.endo_names(1:6,:),statistic1,10,8,4);