//Version V3

var log_c log_n1 log_n2 log_l1 delta  log_y log_z log_A log_w1 log_w2;
varexo eta eta3;

external_function(name=display);

parameters beta  rho rho2 rho3 rho4 epsilon xi
   MeanZ MeanAlpha MeanA MeanGamma popn log_alpha log_gamma thetta thetap ;
beta=.95;
epsilon = 2;
xi=2;

popn=1;
MeanZ=.5;
MeanAlpha=3; 
MeanA=1;
MeanGamma=1;
log_alpha=log(MeanAlpha);
log_gamma=log(MeanGamma);
thetta=2;
thetap=1;

rho = 0.9;
rho2=0.95;
rho3=0.90;
rho4=.95;

display(0.12^(thetta-1));
display(0.12^(2-1));

model;


// MRS between labor(s) and exp(log_c)onsumption
exp(log_z)^(1/xi)*(exp(log_A)^(xi-1)*exp(log_y)/(exp(log_l1(-1))*exp(log_n1)))^(1/xi)=exp(log_c)*exp(log_gamma)*exp(log_n1)^(1/epsilon);
(1-exp(log_z))^(1/xi)*(exp(log_A)^(xi-1)*exp(log_y)/((popn-exp(log_l1(-1)))*exp(log_n2)))^(1/xi)=exp(log_c)*exp(log_gamma)*exp(log_n2)^(1/epsilon);

//Firm MPLs
exp(log_w1)-(exp(log_z)*exp(log_A)^(xi-1)*exp(log_y)/(exp(log_l1(-1))*exp(log_n1)))^(1/xi);
exp(log_w2)-((1-exp(log_z))*exp(log_A)^(xi-1)*exp(log_y)/((popn-exp(log_l1(-1)))*exp(log_n2)))^(1/xi);





//EE
delta^(thetta-1)*exp(log_alpha)*thetta=exp(log_c)*beta*(
        -exp(log_gamma(+1))*epsilon/(1+epsilon)*(exp(log_n1(+1))^((1+epsilon)/epsilon)-exp(log_n2(+1))^((1+epsilon)/epsilon))
        +thetta*exp(log_alpha(+1))*delta(+1)^(thetta-1)/exp(log_c(+1))
        +(exp(log_z(+1))^(1/xi)*exp(log_n1(+1))*(exp(log_A(+1))^(xi-1)*exp(log_y(+1))/(exp(log_l1)*exp(log_n1(+1))))^(1/xi)
            -(1-exp(log_z(+1)))^(1/xi)*exp(log_n2(+1))*(exp(log_A(+1))^(xi-1)*exp(log_y(+1))/((popn-exp(log_l1))*exp(log_n2(+1))))^(1/xi))/exp(log_c(+1))
        );

// Production fcn
exp(log_y)=exp(log_A)*(exp(log_z)^(1/xi)*(exp(log_n1)*exp(log_l1(-1)))^((xi-1)/xi)
        +(1-exp(log_z))^(1/xi)*(exp(log_n2)*(popn-exp(log_l1(-1))))^((xi-1)/xi))^(xi/(xi-1));

// Budget constraint
exp(log_c)+exp(log_alpha)*delta^(thetta)=exp(log_y);

// Law of motion
exp(log_l1)=exp(log_l1(-1))+delta;

// EXOGENOUS PROCESS FOR TECHNOLOGY
   log_z=log(MeanZ)*(1-rho)+rho*log_z(-1)+eta;
   log_A=log(MeanA)*(1-rho3)+rho3*log_A(-1)+eta3;
end;
    
initval;
@#for num in [ "1", "2" ]
  log_n@{num}=log(1);
@#endfor
  delta=0;
  log_l1=log(popn/2);
 
end;

shocks;
    var eta3; stderr .01;
    var eta ; stderr .01;
end;

resid;
steady;

stoch_simul(hp_filter=1600, irf=30, periods=200, ar=1, nocorr,order=1);