var y n c i Q sqstar k a sigma_a;

varexo e u;

parameters beta b rho_a rho_z z sigma_z sigma_ss delta alpha theta sqmax psi nu eta kappa;

beta=.9615;
b=2.58;
alpha=0.36;
delta=.1;
rho_a=.9225;
rho_z=.706;
sigma_z=.017;
sigma_ss=.021;
theta=1/5.98;
sqmax=.0223;
psi=delta^theta;
nu=-theta/(1-theta)*delta;
eta=1;
kappa=.1;
z=1;

model;

exp(sqstar)=theta/(1-theta)*exp(z)^((1-theta)/theta)*(psi*exp(Q))^(1/theta);
exp(i)=(psi*exp(Q))^(1/theta)*exp(z)^((1-theta)/theta)*exp(k(-1))*(exp(sqstar)/sqmax)^eta;
exp(k)=(1-delta+nu)*exp(k(-1))+psi/(1-theta)*exp(k(-1))*(exp(z)*exp(i)/exp(k(-1)))^(1-theta)*(exp(sqstar)/sqmax)^(eta*theta);
exp(y)=(exp(a)*exp(n))^(1-alpha)*exp(k(-1))^alpha;
exp(y)=exp(i)+exp(c)+exp(k(-1))*eta/(eta+1)*exp(sqstar)*(exp(sqstar)/sqmax)^eta;
b*exp(c)=(1-alpha)*exp(a)^(1-alpha)*(exp(k(-1))/exp(n))^alpha;
exp(Q)=beta*exp(c)/exp(c(+1))*(alpha*exp(a(+1))^(1-alpha)*(exp(n(+1))/exp(k))^(1-alpha)+(1-delta+nu)*exp(Q(+1))+exp(sqstar(+1))*(exp(sqstar(+1))/sqmax)^eta*(1/(eta+1)));
a=rho_a*a(-1)+sigma_a*e;
sigma_a=(1-rho_z)*sigma_ss+rho_z*sigma_a(-1)+kappa*u;

end;

initval;
a=0;
sqstar=-13; %log(((1-theta)^theta*theta^(1-theta)*sqmax^eta/psi*(delta-nu))^(1+eta-theta));
Q=-1.49; %log(1/psi*(exp(sqstar)*(1-theta)/theta)^theta);
n=-1.1; %n =log(1/3);
k = 2.5; %log(((1-beta*(1-delta+nu))*exp(Q)/(beta*alpha))^(1/(alpha-1))*exp(n));
y=.2; %y = log(exp(n)^(1-alpha)*exp(k)^alpha);
i = -1.3; %log(exp(k)*(delta-nu)*(1-theta)*exp(Q));
c = -.06; %log(exp(y) - exp(i));
sigma_a=sigma_ss;

end;

shocks;
var e = 1;
var u = 1;
end;

stoch_simul(order=3);

ex_ = zeros(40,2);
ex_(1,2) = 1;
irf2 = simult_(oo_.steady_state,oo_.dr,ex_,2);

figure;
plot(irf2(3,:)-oo_.steady_state(3));
title('Consumption');