%
% RBC Model with Habit Formation

close all;

% (1) Variables
%----------------------------------------------------------------


% Output (y)
% Consumption (c)
% Capital (k)
% Investment (i)
% Employment (l)
% Interest rate (r)
% Wage (w)
% Technology (z)
% Government expenditure (g)
% Technology shock (e_z)
% Government spending shock (e_g)


var y c k i l r w z g;
varexo e_z e_g;

parameters beta psi gbar lambda delta alpha rho_z rho_g gbar lambda;


% (2) Calibration
%----------------------------------------------------------------


alpha   = 0.36;        % capital share
beta    = 0.99;        % discount factor
delta   = 0.025;       % capital depreciation rate
psi     = 1.75;        % disutility of labor scaling factor
rho_z   = 0.95;        % technology coefficient
rho_g   = 0.90;        % government expenditure coefficient
gbar    = 0.22;        % mean government expenditure
lambda  = 0.00;        % habit formation parameter
                       % lambda>0 (habit formation in consumption)
                       % lambda=0 (baseline model without habit formation)

% (3) Model
%----------------------------------------------------------------


model;
%-----Euler equation----------
 1/(c-lambda*c(-1))-lambda*beta*(1/(c(+1)-lambda*c)) = beta*(1+alpha*exp(z(+1))*k^(alpha-1)*l(+1)^(1-alpha)-delta)*(1/(c(+1)-lambda*c)-lambda*beta*(1/(c(+2)-lambda*c(+1))));

%-----Labor market condition--
 psi/(1-l) = exp(z)*(1-alpha)*k(-1)^alpha*l^(-alpha)*(1/(c-lambda*c(-1))-lambda*beta*(1/(c(+1)-lambda*c)));
 
%-----Accounting identity----- 
 y = c+i+g;

%-----Production function-----
 y = exp(z)*k(-1)^alpha*l^(1-alpha);

%-----Capital accumulation----
 i = k-(1-delta)*k(-1);

%-----Technology process------
 z = rho_z*z(-1)+e_z;

%----Government spending------
 log(g) = (1-rho_g)*log(gbar)+rho_g*log(g(-1))+e_g;

%----Interest rate------------

 r = 1+alpha*exp(z)*k(-1)^(alpha-1)*l^(1-alpha)-delta;

%----Wage---------------------

 w = (1-alpha)*exp(z)*k(-1)^alpha*l^(-alpha);
end;


% (4) Simulation
%----------------------------------------------------------------


initval;
 k   = 10;
 c   = 0.7;
 l   = 0.3;
 z   = 0; 
 e_z = 0;
 e_g = 0;
 g   = 0.2;
end;


shocks;
var e_z; stderr 0.008;
var e_g; stderr 0.0074;
end;

steady;

stoch_simul(hp_filter = 1600, order = 1);


% (5) 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);