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;

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);