var c k y g h u z y_h i w r linv ly lcons lw lh lk ly_h;
varexo eps_z eps_g;

parameters mu_g sigma rho_g delta psi alpha rho_z beta Ibar Ybar Cbar Wbar Hbar Kbar YHbar;

% Parameter Values 

beta = 0.95;
psi = 0.33;
alpha = 0.68;
sigma =1.5;
delta = 0.05;
mu_g = log(1.0066);
rho_z = 0.6;
rho_g = 0.1;


% Compute Steady State Values

Rbar=(1/(beta*exp(mu_g)^(psi*(1-sigma)-1))-(1-delta));
YKbar=Rbar/((1-alpha)*exp(mu_g));
IYbar=(1/YKbar)*(1-(1-delta)*exp(mu_g)^(-1));
CYbar=1-IYbar;
Hbar=(((1-psi)*CYbar+psi*alpha)/(psi*alpha))^(-1);
Kbar=(((exp(mu_g)^(alpha-1)*(Hbar^alpha))/YKbar))^(1/alpha);
Ybar=Kbar^(1-alpha)*Hbar^alpha*exp(mu_g)^(alpha-1);
Cbar=CYbar*Ybar;
Ibar=IYbar*Ybar;
Wbar=alpha*(Ybar/Hbar);
YHbar = Ybar/Hbar;



Model;
y = exp(z)*k(-1)^(1-alpha)*exp(g)^(alpha-1)*h^alpha; 
z = rho_z*z(-1)+eps_z;     
g = (1-rho_g)*mu_g+rho_g*g(-1)+eps_g;   
u = (c^psi*(1-h)^(1-psi))^(1-sigma)/(1-sigma);  
r = (1-alpha)*y/k(-1)*exp(g);
w = alpha*y/h;
c^(psi*(1-sigma)-1)*(1-h)^(1-psi)*(1-sigma) = exp(g(+1))^(psi*(1-sigma)-1)*beta*(c(+1)^(psi*(1-sigma)-1)*(1-h(+1))^(1-psi)*(1-sigma)*(1+r(+1)-delta));
(1-psi)*c = psi*(1-h)*w;
y_h = y/h;
c + k = y + (1-delta)*k(-1)*exp(g)^(-1);
i+ c = y;

linv = log(i) - log(Ibar);
ly =log(y) - log(Ybar);
lcons =log(c) - log(Cbar);
lw =log(w) - log(Wbar);
lh = log(h) - log(Hbar);
lk = log(k) - log(Kbar);
ly_h = log(y_h) - log(YHbar);
end; 

initval;
z = 0;
g = mu_g;
c = Cbar;
k = Kbar;
y = Ybar;
h = Hbar;
i = Ibar;
u = (c^psi*(1-h)^(1-psi))^(1-sigma)/(1-sigma);
y_h=y/h;
w=Wbar; 
r = Rbar;
linv = 0;
ly = 0;
lcons = 0;
lw = 0;
lh = 0;
lk = 0;
ly_h = 0;
end;

shocks;
var eps_g; stderr 1; 
var eps_z; stderr 1;
end;
resid(1);
steady;

check;


stoch_simul(order=1,periods=10000,drop=1800,nomoments,nofunctions) y h  w y_h c i k;