var c k y B h u A y_h i w r linv ly lcons lw lh lk ly_h lb;
varexo eps_a eps_b;

parameters gamma delta psi alpha beta Ibar Ybar Cbar Wbar Hbar Kbar YHbar;

% Parameter Values 

beta = 0.95;        
psi = 0.33;          
alpha = 0.68;         
delta = 0.05; 




% Compute Steady State Values

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



Model;
y = k(-1)^(1-alpha)*exp((A*B)/(A(-1)*B(-1)))^(alpha-1)*h^alpha; 
A = gamma+A(-1)+eps_a;  
B = B(-1)+eps_b;       
u = log(c)-(h/B)^(1+1/psi)/(1+1/psi);
r = (1-alpha)*y/k(-1)*exp((A*B)/(A(-1)*B(-1)));
w = alpha*y/(h*B);  
lb =  c^(-1);
lb= exp(A(+1)*B(+1)/(A*B))^(-1)*beta*(lb(+1)*(1+r(+1)-delta));   
h^(1/psi)= c^(-1)*w;  
y_h = y/h;
c + k = y + (1-delta)*k(-1)*exp((A*B)/(A(-1)*B(-1)))^(-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;

c = Cbar;
k = Kbar;
y = Ybar;
h = Hbar;
i = Ibar;
u = log(c)-((h/B)^(1+1/psi))/(1+1/psi);
y_h=y/h;
w=Wbar; 
r = Rbar;
A=0;
B=0;
linv = 0;
ly = 0;
lcons = 0;
lw = 0;
lh = 0;
lk = 0;
ly_h = 0;
end;

shocks;
var eps_a; stderr 1; 
var eps_b; stderr 1;
end;
resid(1);
steady;

check;


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