clear all;

var c k y g h u z y_h i w r lambda;

varexo eps_z eps_g;

parameters mu_g sigma rho_g delta psi alpha rho_z beta;

% Parameter Values 

beta = 1/1.02;
psi = 0.36;
alpha = 0.68;
sigma = 2;
delta = 0.05;
mu_g = log(1.0066);
rho_z = 0.95;
rho_g = 0.01;

% Compute Steady State Values

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

Model;
y = exp(z)*k^(1-alpha)*(exp(g)*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); 
i = y - c;
r = exp(z)*(1-alpha)*k^(-alpha)*(exp(g)*h)^alpha;
w = exp(z)*alpha*k^(1-alpha)*h^(alpha-1)*exp(g)^alpha;
y_h = y/h;
lambda = psi*c^(psi*(1-sigma)-1)*(1-h)^(1-psi)*(1-sigma);
lambda = beta*lambda(+1)*(1+r(+1)-delta);
(1-psi)*c = psi*(1-h)*w;
c+k(+1)=y+(1-delta)*k;
end; 

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


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

check;

stoch_simul(order=1);