// Jing Han's Two Sector Model with Monetary Policy
close all;
var C,     Lc,     Li,    Ic,    Ii,       Kc,   Ki,   Kai,  W,   Qc, 
    Qi,    Dc,     Di,    Lam,   OmeC,     OmeI, KsiC, KsiI, Inf, InfI, 
    R,     Y1,     Y2,    L,     K,        Y,    Ac,   Ai,   Iv,  Q,
    ec,    elam,   ekc,   eki,   eac,      eai,  einf, einfi,ey,  eqc,
    eqi,   eomec,  eomei, eksii, eksic,    eelam,er,   eec,  ekai,eekai,
    eic,   eii,    ey2,   eey2,  eer,      eq,
    gc_obs,inf_obs,r_obs, gi_obs,infi_obs;
varexo epsi1,epsi2,epsi3,err1,err2,err3,err4;

parameters beta, tao, alpha1, alpha2, delta1, delta2, sigma,  theta1, theta2, gamma,
           Piss, eta,  phip1, phip2, phik1, phik2, rhoy, rhopi, rhop2,rhor, rhoac, 
           sigac,rhoai,sigai,sigm,rhoM,rhoG,InfIss,R_ss,rhop,b,g,lq,rhoq;


beta=.99;
b=.7;
eta=5;
Piss=1;
InfIss=1;
tao=.20;
alpha1=1/3;
alpha2=1/3;
delta1=.025;delta2=.025;
sigma=1;
theta1=6;theta2=6;
gamma=0.2;
phik1=10;phip1=30;phip2=30;phik2=10; rhop2=-2;rhoq=0;
rhor=0.6;rhopi=2;rhoy=0.2;
rhoac=.8;rhoai=.8;rhoG=.7;
R_ss=Piss/beta;
g=10;
lq=.9;
sigac=0.01;
sigai=0.01;
sigm=0.01;


model;
# rhop=(rhopi+1)*(1-rhor);
# rhoy1=rhoy*(1-rhor);
# rhopi2=(rhop2)*(1-rhor);
# rhoq1= rhoq*(1-rhor);
 //1. 
ec=c(+1);
 //2.
 elam=lam(+1);
// 3. 
ekc=kc(+1);
// 4. 
eki=ki(+1);
// 5. 
eac=ac(+1);
// 6. 
eai=ai(+1);
// 7. 
einf=inf(+1);
// 8. 
einfi=infi(+1);
// 9.
 ey=y(+1);
//10. 
eqc=qc(+1);
//11. 
eqi=qi(+1);
//12.
 eomec=omec(+1);
//13. 
eomei=omei(+1);
//14. 
eksii=ksii(+1);
//15. 
eksic=ksic(+1);
//16. 
eelam=lam(+2);
//17.
 er=r(+1);
//18. 
eec=c(+2);
ekai=kai(+1);
eekai=kai(+2);
eic=ic(+1);
eii=ii(+1);
//ey1=y1(+1);
ey2=y2(+1);
eey2=y2(+2);
eer=r(+2);
eq=q(+1);
 //(ec(-1)-b*c(-1))^(-sigma)-beta*b*(eec(-1)-b*ec(-1))^(-sigma)       =eLam(-1);
(c-b*c(-1))^(-sigma)-beta*b*(eec(-1)-b*c)^(-sigma)       =eLam(-1);
 eta*((Lc+Li)^gamma-g/2*(l/l(-1)-1)^2)=eLam(-1)*W;
 Y1               =exp(Ac)*Kc(-1)^alpha1*Lc^(1-alpha1);
 Y2               =exp(Ai)*Ki(-1)^alpha2*Li^(1-alpha2);
 Dc               =Y1-(1-lq+lq*r)*W*Lc-Kai*Ic-phip1/2*(Inf/Piss-1)^2*Y1;
 eLam(-1)*(1-lq+lq*er(-1))*W         =KsiC*(1-alpha1)*Kc(-1)^alpha1*Lc^(-alpha1);
 eLam(-1)*Kai       =OmeC*(1-phik1*(Ic/Kc(-1)-delta1));
 Di               =Kai*Y2-(1-lq+lq*r)*W*Li-Kai*Ii-phip2/2*(InfI*Inf/(Piss*InfIss)-1)^2*Kai*Y2;
 eeLam(-2)*(1-lq+lq*eer(-2))*W          =eKsiI(-1)*(1-alpha2)*eKi(-2)^alpha2*Li^(-alpha2);
 eelam(-2)*eKai(-1)        =eOmeI(-1)*(1-phik2*(Ii/eKi(-2)-delta2));
 Y1               =phip1/2*(Inf/Piss-1)^2*Y1+C;
 Y2               =phip2/2*(InfI*Inf/(InfIss*Piss)-1)^2*Y2+(Ic+Ii);
 Lam/er(-1)         =beta*eLam/eInf;
 Lam*Qc        =beta*eLam*(eQc+Dc/eInf);
 Lam*Qi        =beta*eLam*(eQi+Di/eInf);
 Q/Y             =Qc/C+Qi/Iv;
 Kc               =(1-delta1)*Kc(-1)+Ic-phik1/2*(Ic/Kc(-1)-delta1)^2*Kc(-1);

// OmeC             =beta*eOmeC*(1-delta1-phik1/2*(Ic(+1)/Kc-delta1)^2+phik1*(eIc/Kc-delta1)*Ic(+1)/Kc)+beta*eKsiC*alpha1*exp(Ac(+1))*Kc^(alpha1-1)*Lc(+1)^(1-alpha1);
   OmeC             =beta*eOmeC*(1-delta1-phik1/2*(eIc/Kc-delta1)^2+phik1*(eIc/Kc-delta1)*eIc/Kc)+beta*eKsiC*alpha1*exp(eAc)*Kc^(alpha1-1)*Lc(+1)^(1-alpha1);
 
Lam*phip1*(Inf/Piss-1)*Inf/Piss*ec(-1)=Lam*(1-theta1)*ec(-1)+KsiC*theta1*eC(-1)-Lam*phip1*(Inf/Piss-1)*Inf/Piss*Y1+beta*eLam*phip1*(eInf/Piss-1)*eInf*ec/Piss;
 //Lam*phip1*(eInf(-1)/Piss-1)*eInf(-1)/Piss*Y1=Lam*(1-theta1)*eC(-1)+KsiC*theta1*eC(-1)-eLam(-1)*phip1*(eInf(-1)/Piss-1)*eInf(-1)/Piss*eC(-1)+beta*eLam*phip1*(eInf/Piss-1)*eInf*ec/Piss;
 Ki               =(1-delta2)*Ki(-1)+Ii-phik2/2*(Ii/Ki(-1)-delta2)^2*Ki(-1);
 OmeI             =beta*eOmeI*(1-delta2-phik2/2*(Ii(+1)/eKi(-1)-delta2)^2+phik2*(Ii(+1)/eKi(-1)-delta2)*Ii(+1)/eKi(-1))+beta*eKsiI*alpha2*exp(eAi)*eKi(-1)^(alpha2-1)*Li(+1)^(1-alpha2);
 //eeLam(-2)*phip2*(eInfI(-1)*eInf(-1)/(InfIss*Piss)-1)*eInfI(-1)*eInf(-1)/(InfIss*Piss)*Kai*(eii(-1)+eic(-1))=eLam(-1)*(1-theta2)*Kai*(eii(-1)+eic(-1))+KsiI*theta2*(eii(-1)+eic(-1))+beta*elam*phip2*(eInfI*eInf/(Piss*InfIss)-1)*eeKai(-1)*(eii+eic)*eInfI*eInf/(Piss*InfIss);
   eeLam(-2)*phip2*(eInfI(-1)*eInf(-1)/(InfIss*Piss)-1)*eInfI(-1)*eInf(-1)/(InfIss*Piss)*Kai*(eii(-1)+eic(-1))=eLam(-1)*(1-theta2)*Kai*(eii(-1)+eic(-1))+KsiI*theta2*(eii(-1)+eic(-1))+beta*eelam(-1)*phip2*(eInfI*eInf/(Piss*InfIss)-1)*eeKai(-1)*(eii+eic)*eInfI*eInf/(Piss*InfIss);
 R/R_ss=(eY(-1)/Y(-1))^rhoy1*(eInf/Piss)^rhop*(eInfI)^rhopi2*(R(-1)/R_ss)^rhor*(q/q(-1))^rhoq1*exp(-epsi3);
 Kai=Kai(-1)*InfI;
 L=Lc+Li;
 K=Kc+Ki;
 Iv=Ic+Ii;
 Y=Y1+Kai*Y2;
 Ac=rhoac*Ac(-1)+epsi1;
 Ai=rhoai*Ai(-1)+epsi2;


 gc_obs=log(C/C(-1))+err1;
 inf_obs=log(inf)+err2;
 r_obs=log(R)+log(inf);
 gi_obs=log(Iv/Iv(-1))+err3;
 infi_obs=log(infi)+err4;

end;


initval;
C=      1.14510493126063;
eai  =1;
eac  =1;
ec= C;
eec= C;
Lc=     0.29072403227433;
Li=     0.13103380441791;
l= Lc+Li;
Ic=     0.15291596637723;
Ii=     0.01604523654428;
Kc=     6.11663865508928;
Ki=     0.64180946177124;
eki=6.11663865508928;
ekc=0.64180946177124;
Kai=    2.00006501098956;
ekai= Kai;
eekai=Kai;
W=      1.80528508772536;
Qc=     51.95913625575858;
Qi=     6.85956157878905;
Q =26.7612;
eq= Q;
Dc=     0.52483976015918;
Di=     0.06928850079585;
Lam=    0.93449580794579;
elam=   Lam;
eelam=  Lam;
OmeC=   1.86905236854460;
OmeI=   1.86905236854460;
KsiC=   0.77874650662149;
eksic= 0.77874650662149;
KsiI=   1.55754364045383;
eksii=  1.55754364045383;
Inf=   Piss;
einf = Piss;
InfI=  InfIss;
einfi=1;
einf=1;
eomec=omec;
eomei=omei;
R=     1.01010101010101;
er=    R;
eer=   R;
Y1=    1.14510493126063;
Y2=    0.33793339017803;
Y=1.5;
ey=1.5;
ey2= Y2;
eey2=Y2;
Iv=Ic+Ii;
gc_obs=0;
inf_obs=0;
r_obs=0;
gi_obs=0;
infi_obs=0;
err1=0;
err2=0;err3=0;err4=0;epsi1=0;epsi2=0;epsi3=0;Ac=0;Ai=0;
end;



endval;
Inf=Piss;
R=R_ss;
InfI=1;
end;

%simul(periods=100);


maxit=200;
steady;
%check;
shocks;
var epsi1;
periods 1;
values 1;
var epsi2;
periods 1;
values 0;
var epsi3;
periods 1;
values 0;
end;

//simul(periods=20);
y_e1 = sigac*y_(:,:);
//pi_e1 = sigac*y_(13,3:22);
//r_e1 = sigac*y_(15,3:22);

shocks;
var epsi1;
periods 1;
values 0;
var epsi2;
periods 1;
values 1;
var epsi3;
periods 1;
values 0;
end;

//simul(periods=21);

//pi_e2 = sigai*y_(13,4:23);
//r_e2 = sigai*y_(15,4:23);

shocks;
var epsi1;
periods 1;
values 0;
var epsi2;
periods 1;
values 0;
var epsi3;
periods 1 2;
values 1;
end;

//simul(periods=21);
y_e1 = sigm*y_(:,:);
//c_e2 =sigm*y_(1,4:23);
//pi_e2 = sigm*y_(13,4:23);
//r_e2 = sigm*y_(15,4:23);

shocks;
var epsi1;stderr sigac;
var epsi2;stderr sigai;
var epsi3;stderr sigm;
var err1;stderr .01;
var err2;stderr .01;
var err3;stderr .01;
var err4;stderr .01;
end;

stoch_simul(hp_filter=1600,periods=20100,qz_criterium=1,order=1,irf=20,replic=50)C  Inf L Kai  Q  R W  Y Y2;
%forecast(periods=20);





%ML estimation;
%estimated_params;
%theta1,6,0,1000;
%theta2,6,0,1000;
%alpha2,.16,.05,.9;
%phip1,30,20,40;
%phip2,2,.5,5;
%phik1,1.5,.01,3;
%phik2,1.5,.01,3;
%rhopi,1.1,0.2,3;
%rhor,0.6,-0.1,0.95;
%rhoy, .2,-.2,.7;

%rhop2,0,-.6,.6;
%rhoac,.6,0.1,0.99;
%rhoai,.6,0.1,0.99;
%sigac,0.01,.0001,.1;
%sigai,0.01,0.0001,.1;
%sigm,.01,0.0001,0.1;
%end;
%alpha1,.45,.1,.9;

//varobs gc_obs,  inf_obs, r_obs, gi_obs,infi_obs;
//estimation(datafile=TSJH_data_b,mh_replic=50,mh_nblocks=3,mh_jscale=.8);
x=[3 59 61 36 42 56 53 58 47];


xname=['Consumption  ';
       'Output       ';
       'Investment   ';
       'Inflation    ';
       'Pri. of Inv  ';
       'Interest Rate';
       'Asset Price  ';
       'Real Wage    ';
       'Average Hour ';];
%Deterministic Simulation
simy_e1=[C_epsi1./ys_(x(1)), Y_epsi1./ys_(x(2)), Y2_epsi1./ys_(x(3)), Inf_epsi1./ys_(x(4)), Kai_epsi1./ys_(x(5)), R_epsi1./ys_(x(6)), Q_epsi1./ys_(x(7)), W_epsi1./ys_(x(8)), L_epsi1./ys_(x(9))];
simy_e2=[C_epsi2./ys_(x(1)), Y_epsi2./ys_(x(2)), Y2_epsi2./ys_(x(3)), Inf_epsi2./ys_(x(4)), Kai_epsi2./ys_(x(5)), R_epsi2./ys_(x(6)), Q_epsi2./ys_(x(7)), W_epsi2./ys_(x(8)), L_epsi2./ys_(x(9))];
simy_e3=[C_epsi3./ys_(x(1)), Y_epsi3./ys_(x(2)), Y2_epsi3./ys_(x(3)), Inf_epsi3./ys_(x(4)), Kai_epsi3./ys_(x(5)), R_epsi3./ys_(x(6)), Q_epsi3./ys_(x(7)), W_epsi3./ys_(x(8)), L_epsi3./ys_(x(9))];

save model_sim simy_e1 simy_e2 simy_e3;

shocks;

var epsi2;
periods 2 ;
values 1 ;

end;





%forecast(periods=20);

figure;
qq=0:size(simy_e1,1);
for iter=1:length(x);
subplot(3,3,iter);
pick=iter;
plot(simy_e1(:,pick));
title(xname(iter,:));
end;
suptitle('Neutral Technology Shock');

figure;
qq=0:size(simy_e2,2);
for iter=1:length(x);
subplot(3,3,iter);
pick=iter;
plot(simy_e2(:,pick));
title(xname(iter,:));
end;
suptitle('Embodied Technology Shock');

figure;
qq=0:size(simy_e3,2);
for iter=1:length(x);
subplot(3,3,iter);
pick=iter;
plot(simy_e3(:,pick));
title(xname(iter,:));
end;
suptitle('Monetary Shock');