% The RED: no habit, no capital 

%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

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

var Y Ll Lb Q B R Cl Cb Lambdab Lambdal A Phi Theta;
varexo ea etheta;

predetermined_variables Ll Lb B R;

parameters L Thetass betal betab gamma sigmab sigmal sigmaw rhoa rhotheta;

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

L=1;
betal=0.99;
betab=0.98;
gamma=0.05;
sigmab=4;
rhoa=0.9;
rhotheta=0.9;
Thetass=1;
sigmal=2;
sigmaw=2;

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

model;
# Lbss=L/2;
# Llss=L-Lbss;
# Qss=betab*gamma*Lbss^(gamma-1)/(1-(Thetass*betal+(1-Thetass)*betab));
# Rss=1/betal;
# Bss=betal*Thetass*Qss*Lbss;
# Clss=(Rss-1)*Bss;
# Lambdalss=Clss^(-sigmal);
# b=Qss*(1-betal)*Lambdalss/(betal*Llss^(-sigmaw));

Ll+Lb=L;
Lambdal=Cl^(-sigmal);
Lambdal*Q=betal*(Lambdal(+1)*Q(+1)+b*((Ll(+1))^(-sigmaw)));
Lambdal=betal*Lambdal(+1)*R(+1);
Cl+Q*(Ll(+1)-Ll)+B(+1)=R*B;
Lambdab=Cb^(-sigmab);
Lambdab*Q=betab*Lambdab(+1)*(Q(+1)+gamma*A(+1)*(Lb(+1))^(gamma-1))+Phi*Theta*Q(+1);
Lambdab=(betab*Lambdab(+1)+Phi)*R(+1);
Cb+Cl=Y;
R(+1)*B(+1)=Theta*Q(+1)*Lb(+1);
Y=A*(Lb^gamma);
log(A)=rhoa*log(A(-1))+ea;
log(Theta)=rhotheta*log(Theta(-1))+etheta;
end;

%----------------------------------------------------------------
% 4. Steasy State
%----------------------------------------------------------------

steady_state_model;
A=1;
Theta=1;
R=1/betal;
Lb=L/2;
Ll=L-Lb;
Y=Lb^gamma;
Q=betab*gamma*Lb^(gamma-1)/(1-(Theta*betal+(1-Theta)*betab));
B=betal*Theta*Q*Lb;
Cl=(R-1)*B;
Cb=Y-Cl;
Lambdal=Cl^(-sigmal);
Lambdab=Cb^(-sigmab);
Phi=(betal-betab)*Lambdab;
%b=Q*(1-betal)*Lambdal/(betal*Ll^(-sigmaw));
end;

%----------------------------------------------------------------
% 5. Simulation
%----------------------------------------------------------------
%model_diagnostics;
steady;

%resid;

check;

%model_info;



shocks;
var ea=1;
var etheta=1;
end;

stoch_simul(loglinear,irf=100,order=1,periods=200);
datatomfile('simudatalevel',[]);
dynatype(typefilename);
rplot Y;