var ci cs ki ks bi bs q r mu w rb a kappa_t BETAi BETAs;
varexo nui nuk;
parameters EPSILON  KAPPA_bar BETTA OMEGA N SIGMA RHO_a RHO_kappa SIGMANUi SIGMANUk ETA ZETAi ZETAs Z1; 


ETA       = 0.022;
EPSILON   = 0.39; 
ZETAi     = 0.99; 
ZETAs     = 0.985;
BETTA     = 0.97;
OMEGA     = 0.1;
KAPPA_bar = 0.86; 
N         = 0.5;
SIGMA     = 2; 
RHO_a     = 0.9;
RHO_kappa = 0.9; 
SIGMANUi  = 0.005;
SIGMANUk  = 0.011;
Z1 =  1;

model;

% Discount factors are endogenous
BETAi = ZETAi*(1+exp(ci))^(-ETA);
BETAs = ZETAs*(1+exp(cs))^(-ETA);

% budget constraint 
exp(ci)+(exp(q)*exp(ki)) = exp(w)+(exp(q)+exp(r))*exp(ki(-1))+ exp(bi)-exp(rb(-1))*exp(bi(-1)); 


%borrowing constraint  
exp(bi) = exp(kappa_t)*exp(q)*exp(ki); 

% budget constraint savers 
exp(cs) + exp(q)*exp(ks) = exp(w)+ exp(q)*exp(ks(-1))+ exp(bs) - exp(rb(-1))*exp(bs(-1)) + Z1*(exp(ks(-1)))^OMEGA; 


% FOC  wrt to ki 
exp(ci)^(-SIGMA)= BETAi *exp(ci(+1))^(-SIGMA)*(exp(q(+1))+exp(r(+1)))/exp(q) + exp(kappa_t)*exp(mu); 


% FOC wrt to bi 
exp(ci)^(-SIGMA) = BETAi * exp(ci(+1))^(-SIGMA)*exp(rb) + exp(mu); 



%FOC wrt to ks (For savers)
exp(cs)^(-SIGMA)= BETAs *exp(cs(+1))^(-SIGMA)*(exp(q(+1))+ OMEGA*Z1*(exp(ks))^(OMEGA-1))/exp(q);  


%FOC  wrt to bs
%exp(cs)^(-SIGMA) = BETAs * exp(cs(+1))^(-SIGMA)*exp(rb); 


% FOC wrt to l in the production function 
exp(w)= exp(a)*(EPSILON)*(exp(ki(-1)))^(1-EPSILON);  


% FOC  wrt to k in the production function equation 15
exp(r)= exp(a)*(1-EPSILON)*(exp(ki(-1)))^(-EPSILON); 


% Capital accumulation  
N*exp(ki)+(1-N)*exp(ks) = 1; 


%bonds market 
N*bi + (1-N)*bs = 0; 


% resource constraint
 N*exp(ci) + (1-N)*exp(cs) = exp(a)*(N*exp(ki(-1)))^(1-EPSILON) + (1-N)*Z1*exp(ks(-1))^OMEGA; 


% AR(1)  equation 21
a = RHO_a*a(-1) - nui;

% AR(1)  equation 22
kappa_t = (1-RHO_kappa)*KAPPA_bar + RHO_kappa*kappa_t(-1) + nuk;
end;

% I put the steady state  values calculated by hand
initval;
BETAi = BETTA;
BETAs = BETTA;
ci = log(0.83);
cs = log(1.48);
ki = log(1.6);
ks = log(0.4);
r = log(0.51);
rb = log(1.03);
w = log(0.29);
a = log(1);
kappa_t = log(1);
bi = 8.85;
bs = -8.85;
%the value of bs_s is negative, therefore no log()
mu = log(0.15);
q = log(7.37);
nui = 0;
nuk = 0;
end;

steady;
resid (1);
check;

shocks;
var nui =SIGMANUi^2;
var nuk = SIGMANUk^2;
end;
stoch_simul(order =1,irf=40, hp_filter=1600);