%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var c ct cn labt labn kt kn pn pc eps i y k kdt kdn;
varexo e;


parameters alpha beta sigma omega delta neta rho epsilon mu At An vega;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

alpha = 0.30;
beta  = 0.95;
sigma = 2;
omega = 0.5;
delta = 0.05;
neta  = 0.50;
mu = -0.3;
At = 1;
An = 1.52;
rho = 0.9;
epsilon = 0.74;
vega = 1;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 
c = (omega*ct^(-mu)+(1-omega)*cn^(-mu))^(-1/mu);
pn = ((1-omega)/omega)*(cn/ct)^(-(1+neta));
pc = (omega^(1/(1+neta))+(1-omega)*(1/(1+neta))*pn^(neta/(1+neta)))^((1+neta)/neta);
k = (kt^(-vega)+kn^(-vega))^(-1/vega);
kdt = (kt^(-vega-1))*(kt^-vega+kn^-vega)^((-1/vega)-1);
kdn = (kn^(-vega-1))*(kt^-vega+kn^-vega)^((-1/vega)-1);
(c^(-sigma)/pc) = beta*(c(+1)^(-sigma)/pc(+1))*(At*exp(eps)*alpha*(kt(+1)/labt(+1))^(alpha-1)+1-delta)/kdt(+1);
At*exp(eps)*(1-alpha)*(kt/labt)^alpha = pn*An*exp(eps)*(1-neta)*(kn/labn)^neta;
(At*exp(eps)*alpha*(kt/labt)^(alpha-1))/kdt = pn*(An*exp(eps)*neta*(kn/labn)^(neta-1))/kdn; 
ct(-1)+(k)-(1-delta)*(k(-1))=At*exp(eps)*kt(-1)^alpha*labt(-1)^(1-alpha);
cn = At*exp(eps)*kn^neta*labn^(1-neta);
labn+labt = 1;
eps = rho*eps(-1)+e;
k = i(-1)+(1-delta)*k(-1);
y = At*exp(eps)*kt^alpha*labt^(1-alpha)+pn*At*exp(eps)*kn^neta*labn^(1-neta);
end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
  ct    =     	0.5;
  cn    =     	0.5;
  c     =       (omega*ct^(-mu)+(1-omega)*cn^(-mu))^(-1/mu);
  pn    =       ((1-omega)/omega)*(cn/ct)^(-(1+neta));
  pc    =       (omega^(1/(1+neta))+(1-omega)*(1/(1+neta))*pn^(neta/(1+neta)))^((1+neta)/neta);
  kt    =       1;
  kn    =    	1;
  k      =      (kt^(-vega)+kn^(-vega))^(-1/vega);
  kdt   =       (kt^(-vega-1))*(kt^-vega+kn^-vega)^((-1/vega)-1);
  kdn   =       (kn^(-vega-1))*(kt^-vega+kn^-vega)^((-1/vega)-1);     
  labt  =       0.5;
  labn  =   	0.5;
  eps   =      	0.1;
end;

shocks;
var e;
stderr epsilon^2;
end;

steady;

stoch_simul(hp_filter = 1600, order = 1);