close all;  
write_latex_dynamic_model;
write_latex_static_model;

% Declaration
var A R k h c n eta $\eta$ lambda $\lambda$ mu $\mu$ ;
varexo eA $\epsilon^A$ eR $\epsilon^r$;
parameters r alpha $\alpha$ theta $\theta$ delta $\delta$ phi $\phi$ omega $\omega$ sigma $\sigma$ psi $\psi$ rhoA $\rho_A$ rhoR $\rho_r$;

%% Calibration, see Schmitt-Grohe & Uribe p168, Neumeyer & Perri p365
alpha=0.38;
theta=1;
delta=0.036;
phi=0.05;
omega=1.6;
sigma=2; %gamma
psi=0.0306;
rhoA=0.4;
rhoR=0.8;
r=0.04;

%% Model
model;

mu*(c-h^omega/omega)^(-sigma)-lambda-eta*mu*(-psi)*(1+c-h^omega/omega)^(-psi-1)=0;
mu*(c-h^omega/omega)^(-sigma)*(-h^(omega-1))+lambda*(1-alpha)*A*k(-1)^alpha*h^(-alpha)/(1+(R(-1)-1)*theta)-eta*mu*(-psi)*(1+c-h^omega/omega)^(-psi-1)*(-h^(omega-1))=0;
lambda*(1+phi*(k-k(-1)))=lambda(+1)*(alpha*A*k(-1)^(alpha-1)*h^(1-alpha)+1-delta+phi*(k(+1)-k));
lambda=lambda(+1)*R;
((c(+1)-h(+1)^omega/omega)^(1-sigma)-1)/(1-sigma)-eta(+1)*(1+c(+1)-h(+1)^omega/omega)^(-psi)+eta=0;
mu(+1)=mu*(1+c-h^omega/omega)^(-psi);

R=1+r*exp(n);
log(A)=rhoA*log(A(-1))+eA;
n=rhoR*n(-1)+eR;

end;

%% Computation
initval;
k=29.3069;
c=3.3751;
h=2.1857;
R=1.04;
n=0;
A=1;
eta=-16;
mu=1.0243;
lambda=10;
end;

shocks;
var eA=0.0198^2; % percent??
var eR=0.006342^2;
end;

steady;
stoch_simul(order=1,periods=2100,irf=40,nofunctions,nocorr);