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 !!