// 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, eeinf,  ew,   eew, elc,
    eli,   A2,     R_e,   mc_c,  mc_i,     eeinfi, el,  eel, gr1, gr2,
    gc_obs,inf_obs,r_obs, gi_obs,infi_obs;
varexo e1,e2,e3,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,b,g,lq,rhoq,f,phiw,phip1_2,phip2_2,phiw_2,
           thetaw;


load fixed_parameter;
//param=load( 'last_good');
load fordynare param;
//load last_good last_good;
//param=last_good;
load numcount;
load param_ind ;
load param_ovr ;
load param_det ;
//*****************************************************************************************
//*****************************************************************************************
//*****************************************************************************************

b=0.6861; theta1=6; theta2=6; 
beta=.99;
b=.7;
phiw=100;
eta=5;
Piss=1;
InfIss=1;
tao=.20;
alpha1=1/3+1/10;
alpha2=1/3-0.05;
delta1=.025;delta2=.025;
sigma=1;
theta1=6;theta2=6;
gamma=0.2;
phik1=10;phip1=2;phip2=2;phik2=10; rhop2=-0.3;rhoq=1; phip1_2=50;phip2_2=100;phiw_2=300;
rhor=0.799;rhopi=4;rhoy=0.1;
rhoac=.1;rhoai=.1;rhoG=.7;
R_ss=Piss/beta;
f=0;
g=3;
lq=.5;
sigac=0.003;
sigai=0.003;
sigm=0.002;
thetaw=6;
if b>1;b=1-e-10;end;

if b<0;b=0;end;

if lq>1;lq=1;end;

if lq<0;lq=0;end;

if sigac<=0;sigac=1e-9;end;

if sigai<=0;sigai=1e-9;end;

if sigm<=0;sigm=1e-9;end;

if phip1<=0;phip1=1e-9;end;

if phip2<=0;phip2=1e-9;end;

if phik1<=0;phik1=1e-9;end;

if phik2<=0;phik2=1e-9;end;

if phip1_2<=0;phip1_2=1e-9;end;

if phip2_2<=0;phip2_2=1e-9;end;

if rhoac<=0;rhoac=1e-9;end;

if rhoai<=0;rhoai=1e-9;end;

if rhor<=0; rhor=1e-9;end;

//if rhop2<-0.35; rhop2=-0.35+1e-9;end;


model;
# rhop=(rhopi)*(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);
eeinf=Inf(+2);
eeinfi=InfI(+2);
ew=W(+1);
eew=W(+2);
elc=Lc(+1);
eli=Li(+1);
R_e=((Q-Q(-1))*K+(Dc+Di))/(Q*K);
mc_c=KsiC/Lam;
mc_i=KsiC/Lam;
el=l(+1);
eel=l(+2);
Ac=gr1+Ac(-1);
A2=gr2+A2(-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)*eW(-1);
  // eta*(Lc+Li)^gamma=elam(-1)*(ew(-1)-phiw*(ew(-1)/W(-1)-1)*(Lc+Li))+beta*eelam(-1)*(eew(-1)/ew(-1)-1)*eew(-1)*(elc+eli);
 //eta*(L)^gamma=elam(-1)*(w-phiw*(w/W(-1)-1)*(L))+beta*eelam(-1)*(eew(-1)/w-1)*eew(-1)*(el);
eta*eL(-1)*(1+gamma)*thetaw=(thetaw-1)*elam(-1)*[eW(-1)*el(-1)+(phiw*(eW(-1)/W(-1)-1)*eW(-1)/W(-1)+phiw_2*(eW(-1)/W(-1)*W(-2)/W(-1)-1)*eW(-1)*W(-2)/W(-1)^2)*el(-1)]-beta*eelam(-1)*[phiw*(ew/ew(-1)-1)*ew/ew(-1)+w*phiw_2*[ew*w(-1)/ew(-1)^2-1]*ew*w(-1)/ew(-1)^2]*el+beta^2*eelam*[phiw_2*(eew*ew(-1)/ew^2-1)*eew*ew(-1)/ew^2*eel];
 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;
  Dc               =Y1-(1-lq+lq*R)*W*Lc-Kai*Ic-phip1/2*(Inf/Piss-1)^2*Y1-phip1_2*(Inf/Inf(-1)-1)^2*Y1;

 //elam(-1)*(1-lq+lq*er(-1))*W         =KsiC*(1-alpha1)*Kc(-1)^alpha1*Lc^(-alpha1);
elam(-1)*(1-lq-f*Qc+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;
  Di               =Kai*Y2-(1-lq+lq*R)*W*Li-Kai*Ii-phip2/2*(InfI*Inf/(Piss*InfIss)-1)^2*Kai*Y2-phip2_2*(InfI*Inf/(InfI(-1)*Inf(-1))-1)^2*Y2;

// eeLam(-2)*(1-lq+lq*eer(-2))*W          =eKsii(-1)*(1-alpha2)*eKi(-2)^alpha2*Li^(-alpha2);
eelam(-2)*(1-lq-f*Qi+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            =(Qc+Qi)/K;
 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);
 
//elam(-1)*phip1*(eInf(-1)/Piss-1)*eInf(-1)/Piss*ec(-1)=elam(-1)*(1-theta1)*ec(-1)+KsiC*theta1*eC(-1)-elam(-1)*phip1*(Inf/Piss-1)*eInf(-1)/Piss*Y1+beta*eeLam(-1)*phip1*(eeInf(-1)/Piss-1)*eeInf(-1)*eec/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;
elam(-1)*phip1*(einf(-1)/Piss-1)*einf(-1)/Piss*ec(-1)=elam(-1)*(1-theta1)*ec(-1)+KsiC*theta1*ec(-1)-elam(-1)*phip1*(Inf/Piss-1)*einf(-1)/Piss*Y1-elam(-1)*phip1_2*(Inf/Inf(-1)-1)*Inf/Inf(-1)*Y1+beta*eelam(-1)*phip1*(eeinf(-1)/Piss-1)*eeinf(-1)*eec(-1)/Piss+2*beta*phip1_2*eelam(-1)*(eeinf(-1)/einf(-1)-1)*eeinf(-1)/einf(-1)*eec(-1)-beta^2*phip1_2*(eeinf/einf-1)*eeinf/einf*eec;

 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);
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)-phip2_2*(einfi(-1)*einf(-1)/(InfI(-1)*Inf(-1))-1)*einfi(-1)*einf(-1)/(InfI(-1)*Inf(-1))*Kai*(eii(-1)+eic(-1))+2*beta*phip2_2*(einfi*einf/(einfi(-1)*einf(-1))-1)*einfi*einf/(einfi(-1)*einf(-1))*ekai*(eii+eic)-beta^2*phip2_2*(einfi(+1)*einf(+1)/(einfi*einf)-1)*einfi(+1)*einf(+1)/(einfi*einf)*ekai(+1)*(eii(+1)+eic(+1));
 R/R_ss=(ey(-1)/0.738461)^rhoy1*(einf/Piss)^rhop*(einfi)^rhopi2*(R(-1)/R_ss)^rhor*(eq/Q(-1))^rhoq1*exp(e3);
//R/R_ss=(ey(-1)/0.738461)^rhoy1*(eeinf(-1)/Piss)^rhop*(eeinfi(-1))^rhopi2*(R(-1)/R_ss)^rhor*(eq(-1)/4.30752)^rhoq1*exp(e3);
 Kai=Kai(-1)*InfI;
 L=Lc+Li;
 K=Kc+Ki;
 Iv=Ic+Ii;
 Y=Y1+Kai*Y2;
 gr1=rhoac*gr1(-1)+e1;
gr2=rhoai*gr2(-1)+e2;
//Ai=Ac+A2;
Ai=Ac(-1)+A2;
//Ai=Ac;
 //Ai=rhoai*Ai(-1)+e2;


 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;



%simul(periods=100);

//***********************************Just for Estimation***************************************
randn('seed',1234567);
Y_e1=5e-11*normrnd(0,1,[20,1]);Inf_e1=5e-11*normrnd(0,1,[20,1]);R_e1=5e-11*normrnd(0,1,[20,1]);
L_e1=5e-11*normrnd(0,1,[20,1]);C_e1=5e-11*normrnd(0,1,[20,1]); Iv_e1=5e-11*normrnd(0,1,[20,1]);
Kai_e1=5e-11*normrnd(0,1,[20,1]);R_e_e1=5e-11*normrnd(0,1,[20,1]);

Y_e2=5e-11*normrnd(0,1,[20,1]);Inf_e2=5e-11*normrnd(0,1,[20,1]);R_e2=5e-11*normrnd(0,1,[20,1]);
L_e2=5e-11*normrnd(0,1,[20,1]);C_e2=5e-11*normrnd(0,1,[20,1]); Iv_e2=5e-11*normrnd(0,1,[20,1]);
Kai_e2=5e-11*normrnd(0,1,[20,1]);R_e_e2=5e-11*normrnd(0,1,[20,1]);

Y_e3=5e-11*normrnd(0,1,[20,1]);Inf_e3=5e-11*normrnd(0,1,[20,1]);R_e3=5e-11*normrnd(0,1,[20,1]);
L_e3=5e-11*normrnd(0,1,[20,1]);C_e3=5e-11*normrnd(0,1,[20,1]); Iv_e3=5e-11*normrnd(0,1,[20,1]);
Kai_e3=5e-11*normrnd(0,1,[20,1]);R_e_e3=5e-11*normrnd(0,1,[20,1]);

maxit=2000;
steady;
check;


shocks;
var e1;stderr sigac;
var e2;stderr sigai;
var e3;stderr sigm;

end;


stoch_simul(hp_filter=1600,periods=20100,qz_criterium=1,order=1,irf=20,replic=50,noprint)Y C Y2 Inf Kai R  L  R_e Iv;
save ForMatlab Y_e1 Inf_e1 R_e1 L_e1 C_e1 Iv_e1 Kai_e1 R_e_e1... %NT
               Y_e2 Inf_e2 R_e2 L_e2 C_e2 Iv_e2 Kai_e2 R_e_e2... %ET 
               Y_e3 Inf_e3 R_e3 L_e3 C_e3 Iv_e3 Kai_e3 R_e_e3...  %MP;

//save ForMatlab Y_e1 Inf_e1 R_e1 L_e1 C_e1 Iv_e1 Kai_e1... %NT
  //            Y_e2 Inf_e2 R_e2 L_e2 C_e2 Iv_e2 Kai_e2... %ET 
    //          Y_e3 Inf_e3 R_e3 L_e3 C_e3 Iv_e3 Kai_e3...  %MP;


lastparam =  [b theta1 theta2 phik1 phip1 phip2 phik2 rhop2 rhor rhopi rhoy rhoac rhoai sigac sigai sigm phip1_2 phip2_2 lq phiw rhoq  phiw_2 thetaw]';
save lastparam lastparam;
//save est_parameter b theta1 theta2 phik1 phip1 phip2 phik2 rhop2 rhor rhopi rhoy rhoac rhoai rhoG sigac sigai sigm phip1_2 phip2_2 lq phiw;

if sum(abs(Y_e1)>=1e-10)>5 | sum(abs(R_e2)>=1e-10)>5
last_good =       [b theta1 theta2 phik1 phip1 phip2 phik2 rhop2 rhor rhopi rhoy rhoac rhoai sigac sigai sigm phip1_2 phip2_2 lq phiw rhoq phiw_2 thetaw]';
save last_good last_good;
end;