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

parameters 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 of the model

Rbar=(1/beta)-(1-delta);
YKbar=Rbar/(1-alpha);
KYbar=(1-alpha)/Rbar;
IYbar=(1/YKbar)*delta;
CYbar=1-IYbar;
Hbar= (alpha*1/CYbar)^(1+1/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-alpha)*exp(eps_a+eps_b)^(alpha-1)*h^alpha;   %Production function
log(A)=log(A(-1))+eps_a;    %Technology shock
log(B)=log(B(-1))+eps_b;   %%Labor supply shock    
r = (1-alpha)*y/k(-1)*exp(eps_a+eps_b);    
w = alpha*y/(h*log(B));    %% Labor demand
c^(-1)= exp(eps_a(+1)*eps_b(+1))^(-1)*beta*c(+1)^(-1)*(1+r(+1)-delta);   
h^(1/psi)= c^(-1)*w;    %%Labor supply 
y_h = y/h;    %Productivity
c + k = y + (1-delta)*k(-1)*exp(eps_a+eps_b)^(-1);   %budget constraint
i = y-c;



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;
y_h=YHbar;
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(nocheck);

check;


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