// this code is wrriten by luoshuixing for his master paper
// written on april 11th 2014
// w-wage ,l-labor, R-rate, bh-debt of household, pii-inflation
// F-transfer, ch-consumption of household,  T-tax
// Y-GDP, A-productivity, K-private capital, KG-government capital,
// X-mark up, be-debt of firm, ce-consumption of firm, qk-price of capital
// I- investment, ksi-collateral ratio, phi-investment shock coef,
// piitilde-deviation of pii, Xtilde-deviation of X,
//  Rtilde-deviation of R, Ytilde-deviation of Y, GC-consumtion of gov
// GI-investment of gov, GItilde-deviation of GI,
var w l R bh pii F ch  T Y A K KG X be ce I ksi phi piitilde Xtilde
 Rtilde Ytilde GC GI GItilde ksitilde chii Y_obs bh_obs pii_obs GI_obs ;
varexo e_a e_phi e_ksi e_gi;

parameters bet etaa  gam mju alp_g rho_a deltaa rho_phi deltaa_p r_r r_pii rho_ksi rho_y rho_gy omegaa ;
// fixed parameters
bet=0.99;
gam=0.98;
mju=0.5;
deltaa=0.025;
omegaa=0.67;
r_r=0.999;
r_pii=1.01;

// Estimated parameters initialisation
etaa=1.5;
alp_g=0.05;
rho_a=0.5;
rho_phi=0.4;
deltaa_p=0.6;
rho_ksi=0.6;
rho_y=-0.5;
rho_gy=-2.0;


model;
# aa=(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))));
# cece=(mju/1.2-deltaa*aa-(1/bet-1)*0.60*aa);
# chch=((1-mju)/1.2+(1/bet-1)*0.60*aa+1/6-0.16);
# Spii=log(1);
# SX=log(1.2);
# SR=log(1/bet);
# Sksi=log(0.60);
# SY=(mju/(1-mju-alp_g))*log(aa)+(1-mju)/(etaa*(1-mju-alp_g))*log((1-mju)/(1.2*chch)+alp_g/(1-mju-alp_g)*log(aa/18));
# SK=log(aa*exp(Y));
# SI=log(deltaa*aa*exp(Y));
# SKG=log(1/18*aa*exp(Y));
# Sce=log(cece*exp(Y));
# Sch=log(chch*exp(Y));
# Sl=log(((1-mju)/(1.2*chch))^(1/etaa));
# SGI=log(1/18*deltaa*aa*exp(Y));
# SGC=log((1-1/18*deltaa*aa-deltaa*aa-cece-chch)*exp(Y));
# Schii=1/exp(ce)*log((1-gam/bet));
# Sw=(1-mju)*exp(Y)/(1.2*exp(l));
# Sbh=log(0.60*aa*exp(Y));
# Sbe=log(0.60*aa*exp(Y));
# Sphi=log(1);
# Sqk=log(1);
# SA=log(1);
# ST=log(0.16*exp(Y));
# SF=log(1/6*exp(Y));
# Spiitilde=0;
# SXtilde=0;
# SRtild=0;
# SYtilde=0;
# SGItilde=0;
# Sksitilde=0;
// for the household
   // budget constraint
exp(w)*exp(l)+exp(R(-1))*exp(bh(-1))/exp(pii)+exp(F)=exp(ch)+exp(bh)+exp(T);
exp(F)=(exp(X)-1)/exp(X)*exp(Y);
  // F.O.C for household
1/exp(ch)=(1/exp(ch(+1)))*(exp(R)*bet/exp(pii(+1)));
exp(w)/exp(ch)=exp(l)^(etaa-1);
//for the firm
  //production technology
exp(Y)=exp(A)*exp(K(-1))^mju*exp(l)^(1-mju)*exp(KG(-1))^alp_g;
A=rho_a*A(-1)+e_a;
  // budget constraint
exp(Y)/exp(X)+exp(be)=exp(ce)+exp(R(-1))*exp(be(-1))/exp(pii)+exp(w)*exp(l)+exp(I);
exp(be)=exp(ksi)*exp(K(-1));
  // law of motion
exp(K)=(1-deltaa)*exp(K(-1))+exp(I)*exp(phi);
phi=rho_phi*phi(-1)+e_phi;
  //F.O.C for the firm
1/exp(ce)=gam*exp(R)/(exp(ce(+1))*exp(pii(+1)))+exp(chii);
gam*mju*exp(Y(+1))/(exp(ce(+1))*exp(K)*exp(X(+1)))+gam*exp(chii(+1))*exp(ksi(+1))+gam*(1-deltaa)/(exp(phi(+1))*exp(ce(+1)))-1/(exp(phi)*exp(ce));
(1-mju)*exp(Y)/(exp(X)*exp(l))=exp(w);
// retailer pricing
piitilde=bet*piitilde(+1)/(1+bet*deltaa_p)+deltaa_p*piitilde(-1)/(1+bet*deltaa_p)-((1-omegaa)*(1-omegaa*bet)*Xtilde)/(omegaa*(1+bet*deltaa_p));
//goverment and center bank
r_r*piitilde(-1)+(1-r_r)*r_pii*piitilde(-1);
ksitilde=rho_ksi*ksitilde(-1)+rho_y*Ytilde(-1)+e_ksi;
exp(T)=exp(GC)+exp(GI);
exp(KG)=(1-deltaa)*exp(KG(-1))+exp(GI);
GItilde=rho_gy*Ytilde(-1)+e_gi;
exp(Y)=exp(ch)+exp(ce)+exp(GC)+exp(GI)+exp(I);
exp(be)=exp(bh);
Ytilde=Y-SY;
Xtilde=X-log(1.2);
piitilde=pii-log(1);
Rtilde=R+log(bet);
GItilde=GI-SGI;
ksitilde=ksi-log(0.6);
// Measurment equations
 Y_obs=exp(Ytilde);
 bh_obs=exp(bh-Sbh);
 pii_obs=exp(pii);
 GI_obs=exp(GItilde);
end;

steady_state_model;
pii=log(1);
X=log(1.2);
R=log(1/bet);
ksi=log(0.60);
phi=log(1);
qk=log(1);
A=log(1);

piitilde=0;
Xtilde=0;
Rtild=0;
Ytilde=0;
GItilde=0;
ksitilde=0;
Y=(mju/(1-mju-alp_g))*log(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))+(1-mju)/(etaa*(1-mju-alp_g))*log((1-mju)/(1.2*((1-mju)/1.2+(1/bet-1)*0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*0.60*(1-gam/bet))))+1/6-0.16))+alp_g/(1-mju-alp_g)*log((gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))/18));
K=log((gam*mju/(1.2*(1-gam*(1-deltaa)-gam*0.6*(1-gam/bet))))*exp(Y));
I=log(deltaa*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*0.6*(1-gam/bet))))*exp(Y));
KG=log(1/18*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))*exp(Y));
ce=log((mju/1.2-deltaa*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))-(1/bet-1)*0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet)))))*exp(Y));
ch=log(((1-mju)/1.2+(1/bet-1)*0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))+1/6-0.16)*exp(Y));
l=log(((1-mju)/(1.2*((1-mju)/1.2+(1/bet-1)*0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))+1/6-0.16)))^(1/etaa));
GI=log(1/18*deltaa*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))*exp(Y));
GC=log((1-1/18*deltaa*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))-deltaa*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))-((1-mju)/1.2+(1/bet-1)*0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))+1/6-0.16)-((1-mju)/1.2+(1/bet-1)*0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))+1/6-0.16))*exp(Y));
chii=1/exp(ce)*log((1-gam/bet));
w=(1-mju)*exp(Y)/(1.2*exp(l));
bh=log(0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))*exp(Y));
be=log(0.60*(gam*mju/(1.2*(1-gam*(1-deltaa)-gam*exp(ksi)*(1-gam/bet))))*exp(Y));
T=log(0.16*exp(Y));
F=log(1/6*exp(Y));
end;

// shocks
shocks;
var e_a;
stderr 0.02;
var e_phi ;
stderr 0.06;
var e_ksi;
stderr 0.015;
var e_gi;
stderr 0.05;
end;

varobs Y_obs bh_obs pii_obs GI_obs;
//  ESTIMATION
// Definition: 
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF

estimated_params;
stderr e_a,0.5,INV_GAMMA_PDF,0.1,0.02;
stderr e_phi,0.5,INV_GAMMA_PDF,0.4,0.06;
stderr e_ksi,0.5,INV_GAMMA_PDF,0.1,0.02;
stderr e_gi,0.3,INV_GAMMA_PDF,0.2,0.02;

etaa,normal_pdf,1.5,0.1;
alp_g,BETA_PDF,0.05,0.1;
rho_a,BETA_PDF,0.5,0.2;
rho_phi,BETA_PDF,0.3,0.2;
deltaa_p,BETA_PDF,0.6,0.2;
rho_ksi,BETA_PDF,0.6,0.2;
rho_y,NORMAL_PDF,-0.5,0.2;
rho_gy,NORMAL_PDF,-2.0,2;
end;

options_.noconstant = 1;
estimation(datafile=paperobs,mode_compute=4,loglinear, presample=0,prefilter=0,mh_replic=50000,mh_nblocks=4,mh_jscale=0.20,mh_drop=0.2,mode_check,moments_varendo);

stoch_simul(irf=40,print, conditional_variance_decomposition=[1, 4, 8, 16, 40]) w l bh pii ch   Y  K  ce  I GI;