theta=[0.999076	5.75	0.76025	0.567383	0.001281	0.0197260625	0.996577	0.00034	0.001028	0.02	0.846875	0.000345]';
%theta= [0.998826	2.250000	1.729	1.629883	0.001281	0.0184760625	0.991577	0.000340	0.001028	0.024966	0.973750	0.000970]';

var c d h ik k l q rf uc vk yy nk dchat dc m da ddem;
varexo ea ed;

parameters	ALFA BETTA DELTA GAM PSI NBAR TAU NU MUA SIGA MUD SIGD; 

set_param_value('BETTA',theta(1));
set_param_value('GAM',theta(2));
set_param_value('PSI',theta(3)); 
set_param_value('TAU',theta(4));
set_param_value('MUA',theta(5));
set_param_value('SIGA',theta(6));
set_param_value('MUD',theta(9));
set_param_value('SIGD',theta(10));

DELTA   = 0.008333333333333; % Depreciation
ALFA    = 0.345000000000000; % Capital share 
NBAR    = 0.22; % 0.180000000000000; % Steady state NUmber of hours worked;
NU      = 0.256884124396632; % 0.204786351767522; % Share of consumption;
			
model;
#a0		= (exp(MUA) + DELTA - 1)^(1/TAU); % Adjustment cost coefficient
#b0		= (1/(1-TAU))*(exp(MUA) + DELTA - 1); % Adjustment cost coefficient
#THETA  = (1-GAM)/(1-(1/PSI));

0 = - da + MUA + SIGA*ea;
0 = - ddem + MUD + SIGD*ed;

0 = - exp(yy) + exp(c)+exp(ik); 
0 = - exp(yy) + exp(k*ALFA)*exp((da+h)*(1-ALFA));

0 = - exp(vk) + exp((uc(+1)+dchat(+1))*(1-GAM)); 
0 = - exp(uc*(1-1/PSI)) + (1-BETTA) + exp(ddem(+1))*BETTA*exp(vk*(1/THETA));

0 = - exp(m) + (exp(ddem) * BETTA * exp(dc*(-1)) * exp(dchat)^(1-1/PSI) * exp((uc+dchat)*(1/PSI-GAM))/exp(vk(-1)*(1-1/THETA)));
0 = - 1 + exp(rf + m(+1)); 

0 = - (1-ALFA)*exp(yy-h) + (1/NU-1)*exp(c-l);
0 = - 1 + exp(h)+exp(l);

0 = - dchat + NU*(c-c(-1)) +  (1-NU)*(l-l(-1)) + NU*da(-1); // used for ease of notation
0 = - dc + c-c(-1) + da(-1);

0 = - exp(k+da(-1)) + (1-DELTA + exp(nk(-1)))*exp(k(-1));
0 = - exp(nk) + b0 + (a0/(1-(1/TAU)))*(exp((ik-k)*(1-(1/TAU))));
0 = - exp(q) + 1/(a0*(exp((ik-k)*(-1/TAU))));

0 = - exp(d) + ALFA*exp(yy-k) + (exp(nk) - DELTA)*exp(q) - exp(ik-k);
0 = - 1 + ((exp(q(+1))+exp(d(+1)))/exp(q)) * exp(m(+1));

end; 

% Initial values 
steady_state_model;

h = log(NBAR);
l = log(1-NBAR);

da = MUA;
ddem = MUD;
dc = MUA;
dchat = MUA*NU;
m = (log(BETTA) + MUD + MUA*NU*(1-1/PSI) + MUA*(-1));

k = log(((1/(BETTA*exp(MUA*NU*(1-1/PSI))*exp(MUA*(-1))*exp(MUD)) -1+DELTA)/ALFA)^(1/(ALFA-1)))+MUA+h;
yy = log(exp(k)^ALFA * (exp(MUA)*exp(h))^(1-ALFA));
ik = log(exp(k)*(exp(MUA)+DELTA-1));
c = log(exp(yy) - exp(ik));
rf = -m;
q = 0;
nk = log((((exp(MUA) + DELTA - 1)^(1/TAU))/(1-1/TAU))*((exp(ik)/exp(k))^(1-1/TAU)) + ((1/(1-TAU))*(exp(MUA) + DELTA - 1))) ; 
d = log(ALFA*exp(yy)/exp(k) + (exp(nk)-DELTA)*exp(q) - (exp(ik)/exp(k))); 
uc = (PSI/(PSI-1))*log((1-BETTA)/(1-BETTA*exp(MUA*NU*(1-1/PSI))*exp(MUD)));
vk = log((exp(uc))^(1-GAM));
      
end; 

shocks;
var ea; stderr 1;
var ed; stderr 1;
end;

steady(nocheck);
check;
stoch_simul(order=3, periods=0,drop=0,irf=0,noprint,nograph,nocorr,nofunctions,nomoments);