var n k c w z;

varexo e;

parameters alpha beta gamma delta eta rho;

alpha = 0.33;
beta = 0.95;
gamma = 1; %Analytical steady state is valid only if gamma = 1
delta = 0.025;
eta = 2;
rho = 0.9;

%Steady state
n_ss = (((1-alpha)*((1/beta-1+delta)/alpha)^(alpha/(alpha-1))*((1/beta-1+delta)/alpha-delta))/(((1/beta-1+delta)/alpha)^(1/(alpha-1))))^(eta+1);
k_ss = ((1/beta-1+delta)/alpha)^(1/(alpha-1))*n_ss;
c_ss = (1/beta-1+delta)/alpha - delta;

model;

c^(-gamma) = c(+1)^(-gamma)*beta*(alpha*exp(z(+1))*(k/n(+1))^(alpha-1) + 1 - delta);
%n^eta = w*c^(-gamma);
n^eta = (w>1.36654)*w*c^(-gamma) + (w<=1.36654)*1.36654*c^(-gamma);
w = max(exp(z)*(1-alpha)*(k(-1)/n)^alpha,1.36654);
c + k = exp(z)*k(-1)^alpha*n^(1-alpha) + (1-delta)*k(-1);
z = rho*z(-1) - e;

end;

initval;
z = 0;
n = (((1-alpha)*((1/beta-1+delta)/alpha)^(alpha/(alpha-1))*((1/beta-1+delta)/alpha-delta))/(((1/beta-1+delta)/alpha)^(1/(alpha-1))))^(eta+1);
k = ((1/beta-1+delta)/alpha)^(1/(alpha-1))*n;
c = (1/beta-1+delta)/alpha - delta;
w = (1-alpha)*(k/n)^alpha;
end;

steady;
check;


%EP
shocks;
var e; stderr 0.01;
end;
extended_path(periods=200, order=0);