%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dmitri Bogonos
% dmitri.bogonos@rwi-essen.de
% 02.02.2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------
var y theta Lu k c gammaStern rB Ls Ws phi j Lambda i Chi yk;
varexo epsChi epsTheta;
parameters alpha beta delta nu chiN chiS eta rho varphi kappa iota aChi bChi aTheta bTheta;

%----------------------------------------------------------------
% 2. Parameter choice
%----------------------------------------------------------------
beta = 0.99;
varphi = 0.9995;
alpha = 0.33;
eta = 3.4;
nu = 1.67;
iota = (1-alpha)/(nu-1);
kappa = (nu+alpha-2)/(nu-1);
delta = 0.015;
bChi = 0.9;
rho = 0.4;
aChi = log(0.2)*0.1;
aTheta = 0;
bTheta = 0.9;

load parameterfile;
set_param_value('chiN', chiN);
set_param_value('chiS', chiS);

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------
model;
y = theta*(Lu^(1-alpha))*(k^alpha);
1 = beta*c/(c(+1)*gammaStern)*(1 - delta + alpha*y(+1)/(nu*k));
1 = beta*c/(c(+1)*gammaStern)*(1 + rB);
chiN*(Lu^(eta+1)) = (1-alpha)*y/(nu*c);
chiS*(Ls^(eta)) = Ws/c;
phi = Chi*(Ls^(rho-1));
gammaStern^iota = Chi*Ls + varphi;
j - varphi*Lambda*j(+1)*(gammaStern^kappa) = (nu - 1)*y/nu;
Ws = Lambda*phi*j(+1)*(gammaStern^kappa);
k = (1 - delta)*k(-1)/gammaStern + i;
Lambda = beta/gammaStern;
log(Chi) = aChi + bChi*log(Chi(-1)) + epsChi;
log(theta) = aTheta + bTheta*log(theta(-1)) + epsTheta;
y = c + i;
yk = (gammaStern/beta - 1 + delta)*(nu/alpha);
end;


%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------
initval;
Ls = 0.2;
Lu = 0.8;
Chi = 0.2;
theta = 1;
gammaStern = (Chi*Ls + varphi)^(1/iota);
Lambda = beta/gammaStern;
rB = gammaStern/beta - 1;
yk = (gammaStern/beta - 1 + delta)*(nu/alpha);
k = (yk*(Lu^(alpha-1))/theta)^(1/(alpha-1));
y = theta*(Lu^(1-alpha))*(k^alpha);
i = k*(gammaStern - 1 + delta)/gammaStern;
j = ((nu-1)/nu)*y/(1-varphi*Lambda*(gammaStern^kappa));
phi = Chi*(Ls^(rho-1));
Ws = Lambda*phi*j*(gammaStern^kappa);
c = y - i;
end;

resid;
model_diagnostics;
steady;

shocks;
var epsTheta = 0.0049;
end;

stoch_simul(order=1,irf=200);