We have an RBC model with external habits and nonseparability in preferences. Our production function contains Z, which grows at rate g at steady state(Z_t=e^(gt)*e^z_t). With the following dynare code, we dont get the right IRF-s. What might be the problem? We suppose that we might need to detrend the above-mentioned equation, but are not sure how to do it. Thanks in advance.
var Y C I K W R z N Z H F
y c i k w r n q h f;
varexo eps_z;
parameters beta alpha delta gamma a g nu zeta qsi fi phi rho sigma teta Zss Yss Css Iss Kss Wss Rss Hss Fss Nss;
beta = 0.99;
alpha = 0.33;
delta = 0.015;
gamma = 7;
a = 0.01;
nu = 1.5;
zeta = 0.5;
qsi = 0.9;
fi = 0.97;
phi = 0.9;
g = 0.005;
rho = 0.552;
sigma = 0.0065;
//Steady-state values
Nss=1/3;
Zss=exp(g);
Rss=(1/beta)-(1-delta);
Kss=((Rss/alpha)^(1/(alpha-1)))*Nss*Zss;
Wss=(1-alpha)*((Kss/Nss)^alpha)*Zss^(1-alpha);
Yss= (Kss^alpha)*(Zss*Nss)^(1-alpha);
Iss=delta*Kss;
Css = Yss - Iss;
Fss = phi*(1-Nss);
Hss = fi*Css;
//teta = (nu*(Css-Hss)(1-Nss-Fss)^(nu-1))/(Wss(a+(1-Nss-Fss)^nu));
teta = Wss*a/(nu*(Css-Hss)(1-Nss-Fss)^(nu-1)-Wss(1-Nss-Fss)^nu);
model;
//teta=1;
//Euler equation
((C-H)^(-gamma))((a+teta(1-N-F)^nu))^(1-gamma)=beta*((C(+1)-H(+1))^(-gamma))(((a+teta(1-N(+1)-F(+1))^(nu)))^(1-gamma))*(R(+1)+(1-delta));
//Intratemporal optimality
//((C-H)nu(1-N-F)^(nu-1))=teta*W*(a+(1-N-F)^nu);
teta = W*a/(nu*(C-H)(1-N-F)^(nu-1)-W(1-N-F)^nu);
//Wage equation
W=(1-alpha)*K(-1)^alpha*N^(-alpha)*Z^(1-alpha);
//Return to capital
R=alpha*K(-1)^(alpha-1)*(N*Z)^(1-alpha);
//Production function
Y=K(-1)^alpha*(Z*N)^(1-alpha);
//Technologicsal innovation
Z = exp(g)*exp(z);
//Productivity
z=rho*z(-1)+eps_z;
//Resource constraint
Y=C+I;
//Capital accumulation
K=(1-delta)*K(-1)+I;
//Habit formation in consumption
H=(1-zeta)*fi*C(-1)+zeta*H(-1);
//Habit formation in leisure
F=(1-qsi)phi(1-N(-1))+qsi*F(-1);
// Define the variables in logs
y=log(Y);
c=log(C);
i=log(I);
k=log(K);
w=log(W);
r=log(R);
n=log(N);
q=log(Z);
h=log(H);
f=log(F);
end;
// Provide steady state values
initval;
N=Nss;
K=Kss;
I=Iss;
Y=Yss;
W=Wss;
R=Rss;
C=Css;
H=Hss;
F=Fss;
Z=Zss;
y=log(Yss);
c=log(Css);
i=log(Iss);
k=log(Kss);
w=log(Wss);
r=log(Rss);
n=log(Nss);
q=log(Zss);
h=log(Hss);
f=log(Fss);
z=0;
end;
//The starting values for the steady state
resid;
//Compute steady state given the starting values
steady;
check;
shocks;
var eps_z; stderr 0.0065;
end;
stoch_simul(order=1, irf =25) y c i n z;