clear all
periods = 400;
var C, Yc, Yh, H, Nc, Nd, Kc, Kd, L, W, R, Ph, Pl, Xl, Z; 
varexo  ep;

parameters alfa, gama, teta, mh, mc, bet, dk, dh, roh, ac, ah, b;

%Productivity parameters:
alfa = 0.4; gama = 0.15; teta = 0.894; 
% utility and other consumer parameters
mh = 0.044; mc = 0.32; bet = 0.95; dk = 0.05; dh = 0.014;
% shock parameters
roh = 0.95;
%others
ac = ((1-alfa)^(alfa-1)/alfa^alfa);
ah = ((((1-gama)/gama)^(teta*(gama-1))*((1-teta)/(gama*teta))^teta)/(1-gama));
b = 1-mc-mh; 

model;
%Household side
mc/C = L;  
bet*((mh/H(+1))+L(+1)*ph(+1)*(1-dh)) = L*ph; 
b/(1-Nc-Nd) = L*W; 
bet*(L(+1)*(R(+1) + (1 - dk))) = L;  
C(-1) + Kc + Kd - (1-dk)*(Kd(-1) + Kc(-1)) +ph(-1)*(H -(1-dh)*H(1))= R*(Kc(-1)+Kd(-1)) + W*(Nc(-1)+Nd(-1)) + Pl(-1)*Xl(-1); 
%Firm side
(Z^(alfa-1))*R^alfa*(W^(1-alfa))*ac = 1; 
(Z^(teta*(1-gama)))*(R^(teta*gama))*(W^((1-gama)*teta))*(Pl^(1-teta))*ah = Ph; 
(Z^(alfa-1))*(1-alfa)*R^alfa*(W^(-alfa))*ac*Yc = Nc;  
(Z^(teta*(1-gama)))*teta*(R^(teta*gama))*(W^((1-gama)*teta-1))*(Pl^(1-teta))*ah*Yh = Nd; 
(Z^(alfa-1))*alfa*R^(alfa-1)*(W^(1-alfa))*ac*Yc = Kc; 
(Z^(teta*(1-gama)))*gama*teta*(R^(teta*gama-1))*(W^((1-gama)*teta))*(Pl^(1-teta))*ah*Yh = Kd; 
(Z^(teta*(1-gama)))*(1-teta)*(R^(teta*gama))*(W*((1-gama)*teta))*(Pl^(-teta))*ah*Yh = Xl; 
Yc = (Kc^alfa)*(Z*Nc)^(1-alfa);
Yh = (Kd^(gama*teta))*((Z*Nd)^((1-gama)*teta))*(Xl^(1-teta)); %I had an error here in the mod file here
%Technology
Z=roh*Z(-1)*exp(ep);


end;

initval;
%ss values calculated in Mathematica
C = 0.542; 
Yc= 0.674; 
Nc = 0.271; 
Nd = 0.0031; 
Kc = 2.62; 
Kd = 0.0079; 
Pl = 0.00190682;
Yh = 0.0065; 
Ph = 0.932; 
H = 0.461; 
L = 0.58; 
Xl =1; 
W = 1.4859; 
R = 0.10263;
Z = 1; ep = 0;

end;

shocks;
var ep = 0.001;
end;
stoch_simul;