LUCAS model with external habits

Good Morning,

I am working on the Lucas’ model (1978) and I would like to include external habits and more especially the ones of the Campbell and Cochrane (2000) in a Dynare code.
I have changed the code however, the residuals of the equation of “log surplus” are not equal to 0 meaning that something is wrong. Nevertheless I don’t see what might be the error.

%% declaration of variable
var p c rk r ep lb m dEc x s z lambda pdr nu;
%% declaration of exogenous variables
varexo e_a ;
% declaration of deep parameters
parameters gamma rho_a  sigma g smax sbar phi delta alpha;

% calibration
gamma	= 2.00;		% risk aversion
rho_a   = 0.8;		% AR(1) term productivity
sigma	= 0.055;	% std productivity shock
g		= 1.02;		% annual growth rate of consumption (here 2%)
smax    = 0.094;    % Maximum surplus consumption ratio Smax
sbar    = 0.057;    %steady state surplus consumption ratio
phi     = 0.87;     % persistence coefficient
delta   = 0.89;     %subjective discount factor

%%%
%%% Dynamic Equations
%%%
model;
    %surplus consumption ratio
    s=(c-x)/c;
    % sensitivity function 
    lambda =  (1/sbar) * sqrt(1- 2*(log(s)-log(sbar))) - 1;
    % consumption shock
    log(c) = g+log(c(-1))+nu;
    % log surplus consumption ratio evolution
    log(s(+1))=(1-phi)*log(sbar)+phi*log(s)+lambda*(nu);
	% risky asset price
	m(+1) * (1+rk(+1)) = 1;
	% riskless asset rate of return
	m(+1) * (1+r) = 1;
	% marginal utility of consumption U'(c)
	lb = (c*s)^-gamma;
	% pricing kernel between t-1 and t // stochastic discount factor // intertemporal marginal rate of substitution
	m = delta*g^-gamma*(lb/lb(-1));
	% risky asset rate of return between t-1 and t
	1+rk 	= (z+p)/p(-1);
	% equity premium
	ep = rk-r;
    % price dividend ratio
    pdr = m*(c(+1)/c)*(1+pdr(+1));
    %productivity shock 
    log(z) = rho_a*log(z(-1))+sigma*e_a;
    % expected consumption growth
    dEc =  log(g*c(+1)/c);
    % exogenous production
	z = c;
    


end; 

%%% steady state of the model
steady_state_model;
	c	= 1;
	z	= 1;
    nu = -g;
    r	= g^gamma/delta-1;
    rk  = r;
    s   = sbar;
    x   =c*(1-s);
	p	= z/(g^gamma/delta-1);
	lb 	= (c*s)^-gamma;
	m 	= delta*g^-gamma;
	ep 	= 0;
	dEc	= log(g);
    lambda = (1/s) *(1) - 1;
    pdr=(delta*g^-gamma)/(1-delta*g^-gamma);
    
    

end;

%%% varcov matrix
%%% = identity matrix as we assume e~N(0,1)
shocks;
var e_a; stderr 1;


end;

Thanks in advance for your help !!

There is a problem with the growth trend. Take the equation with the residual:

log(s(+1))=(1-phi)*log(sbar)+phi*log(s)+lambda*(nu);

The steady state of s cannot be sbar, because lambda*(nu) is not 0.

I thank you very much for your answer.
In fact I change my previous code in order to model a shock of consumption instead of productivity shock. Here my new code, I still have a mistake but don’t know where it is from. May be you will know … ?
Thanks in advance for your help.

Regards,

CB

close all;

%% declaration of variable
var p rk r ep lb m x s lambda dc;
%% declaration of exogenous variables
varexo e_nu ;
% declaration of deep parameters
parameters gamma sigmanu g smax sbar phi delta;

% calibration
gamma	= 2.00;		% risk aversion
g		= 1.02;		% annual growth rate of consumption (here 2%)
phi     = 0.87;     % persistence coefficient
delta   = 0.89;     %subjective discount factor
sbar    = 0.057;    %steady state surplus consumption ratio
smax    = 0.094;    % Maximum surplus consumption ratio Smax
sigmanu = 0.015^2;

%%%
%%% Dynamic Equations
%%%
model;
    %surplus consumption ratio
    s=1-x;
    % sensitivity function 
    lambda = (1/sbar) * sqrt(1- 2*(log(s)-log(sbar))) - 1;
    % consumption shock
    dc=exp(g+sigmanu*e_nu);
    % log surplus consumption ratio evolution
    log(s)=(1-phi)*log(sbar)+phi*log(s(-1))+log(lambda*dc/exp(g));
	% risky asset price
	m(+1) * (1+rk(+1)) = 1;
	% riskless asset rate of return
	m(+1) * (1+r) = 1;
	% marginal utility of consumption U'(c)
	lb = (s)^-gamma;
	% pricing kernel between t-1 and t // stochastic discount factor // intertemporal marginal rate of substitution
	m = delta*((s/s(-1))*g)^-gamma;
	% risky asset rate of return between t-1 and t
	1+rk 	= (1+p)/p(-1);
	% equity premium
	ep = rk-r;
 
   

end; 

%%% steady state of the model
steady_state_model;
    s   = sbar;
    x = 1-s;
    lambda = 1;
	dc = exp(g);
    lb 	= s^-gamma;
    r	= g^gamma/delta-1;
    rk  = r;
	p	= 1/(g^gamma/delta-1);
	m 	= delta*g^-gamma;
	ep 	= 0;

    

end;

%%% varcov matrix
%%% = identity matrix as we assume e~N(0,1)
shocks;
var e_nu; stderr 1;


end;

%%% checking steady state
resid(1);
%%%% computing the steady state
steady;
%%% checking the stability of the model
check;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Running stochastic simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('RUNNING STOCHASTIC SIMULATIONS')
%%% plotting IRF + second moment statistics
stoch_simul(order=1) p rk r ep;
% varcov matrix reporting
disp('varcov matrix for endogenous variables:')
oo_.var	% this is showing the content of the varcov matrix


%%%% sampling 1000 draws of the model
stoch_simul(order=2,noprint,periods=1000,nograph) rk r ep;
% plotting simulated variables
figure
subplot(2,1,1)
plot(1:options_.periods,100*rk,1:options_.periods,100*r)
title('Risky & Riskless asset returns (pp)')
subplot(2,1,2)
plot(1:options_.periods,100*ep)
ylim([min(100*ep) max(100*ep)])
title('Equity Premium (pp)')

%Gamma
gammas =[1:1:50];
EQPrem=zeros(length(gammas),2);
for i = 1:length(gammas) 
gamma=gammas(i);  
stoch_simul(order=2,nograph,noprint) ep; 
mean_ep = oo_.mean(1);
EQPrem(i)= mean_ep; 
end

figure 
plot(gammas,EQPrem);

I don’t get what you are doing.

    lambda = (1/sbar) * sqrt(1- 2*(log(s)-log(sbar))) - 1;

in steady state reduces to

    lambda = (1/sbar) - 1;

which is not satisfied with lambda=1