// 1) Definition of variables 
 
var y c k invest l y_l z thheta varthheta log_y log_c log_invest log_l log_k log_y_l;
varexo e ethheta evarthheta; 
parameters alppha betta delta psii rho sigma rhothheta sigmathheta rhovarthheta sigmavarthheta;
 
// 2) Calibration
 
alppha = 0.36657;
betta = 0.95563;
delta = 0.03852;
psii = 3.01412;
rho = 0.961879;
sigma = 0.017812;
rhothheta = 0.9;
sigmathheta = 0.01;
rhovarthheta = 0.9;
sigmavarthheta = 0.01;

 
// 3) Model
// In this case Thheta and Varthheta follow an AR(1) process in the natual logs. This means that rhoThheta and rhoVarthheta and sigmaThheta and sigmaVarthheta come from the AR(1) regression of ln Thheta and ln Varthheta.
 
model;
(1/c) = betta*((thheta(+1)/thheta)*(1/c(+1))*(1+alppha*(k^(alppha-1))*((exp(z(+1))*l(+1))^(1-alppha))-delta));
psii*varthheta*c/(1-l) = (1-alppha)*(k(-1)^alppha)*(exp(z)^(1-alppha))*(l^(-alppha));
c+invest = y;
y = (k(-1)^alppha)*(exp(z)*l)^(1-alppha);
invest = k-(1-delta)*k(-1);
y_l = y/l;
z = rho*z(-1)+e;
log(thheta) = rhothheta*log(thheta(-1))+ethheta;
log(varthheta) = rhovarthheta*log(varthheta(-1))+evarthheta;
log_y=log(y);
log_invest=log(invest);
log_c=log(c);
log_l=log(l);
log_k=log(k);
log_y_l=log(y_l);
end;

steady_state_model;
l = ((1 - alppha)*(1/betta - (1 - delta)))/((psii*(1/betta - (1 - delta) - alppha*delta) + (1 - alppha)*(1/betta - (1 - delta))));
k = ((alppha/(1/betta - (1 - delta)))^(1/(1 - alppha)))*l;
invest = delta*k;
y = (k^alppha)*(l^(1-alppha));
c = y - invest;
y_l = y/l;
z = 0;
thheta = 1;
varthheta = 1;
log_y=log(y);
log_invest=log(invest);
log_c=log(c);
log_l=log(l);
log_k=log(k);
log_y_l=log(y_l);
end;
 
// 4) Computation

shocks;
var e = sigma^2;
var ethheta = sigmathheta^2;
var evarthheta = sigmavarthheta^2;
end;

steady;
check;

stoch_simul(order=1,hp_filter=100);