// Estimate a DSGE model with four shocks(TFP shock,Labor_augmented tech shock,IST shock,government expenditure shock) 
// (both shocks feature an anticipated component and an unanticipated component).
//
// This version is the log-linearized equilibrium conditions.
// written by ziguan zhuang

var y c i h u mu_y mu_k x_g q lam k z mu_x mu_a g 
Sz1 Sz2 Sz3 Sx1 Sx2 Sx3 Sa1 Sa2 Sa3 Sg1 Sg2 Sg3;
varexo eps_z0 eps_z1 eps_z2 eps_z3 eps_x0 eps_x1 eps_x2 eps_x3 
eps_a0 eps_a1 eps_a2 eps_a3 eps_g0 eps_g1 eps_g2 eps_g3;

parameters beta theta_c theta_l chi sigma delta0 delta1 delta2 kappa alpha rho_z rho_x rho_a rho_g rho_xg 
y_bar c_bar i_bar h_bar u_bar mu_y_bar mu_k_bar x_g_bar q_bar lam_bar k_bar z_bar mu_x_bar mu_a_bar g_bar;

//calibration
beta = 0.973;
alpha = 0.3;
sigma = 2;
delta0 = 0.025;
u_bar = 1;
mu_a_bar = 0.9957;
mu_y_bar = 1.0045;


//estimation
//theta_c = 0.85;
//theta_l = 0.56;
//chi = 6.1;
//delta2 = 0.11;
//kappa = 5;

//rho_z = 0.89;
//rho_x = 0.14;
//rho_a = 0.52;
//rho_g = 0.98;
//rho_xg = 0.99;

//steady-state
z_bar=1;

mu_x_bar=mu_y_bar*(mu_a_bar^(alpha/(1-alpha)));

mu_k_bar=mu_x_bar*(mu_a_bar^(1/(alpha-1)));

x_g_bar=mu_y_bar^(1/(rho_xg-1));

q_bar=1;

delta1 = (1/beta)*(1-beta*(1-delta0)*mu_a_bar*(mu_y_bar^(-sigma)))*(1/mu_a_bar)*(mu_y_bar^(sigma));

h_bar=(1-alpha)*(mu_y_bar^sigma-beta*theta_c)/((chi/(1-theta_l))*(mu_y_bar-theta_c)*(mu_y_bar^(sigma-1)-beta*theta_l)
*(0.8-(1-(1-delta0)/mu_k_bar)*alpha*beta*mu_y_bar^(1-sigma)*(1-beta*(1-delta0)*mu_a_bar*mu_y_bar^(-sigma))^(-1))
+(1-alpha)*(mu_y_bar^sigma-beta*theta_c));

y_bar=((alpha*beta*mu_y_bar^(-sigma)*(1-beta*(1-delta0)*mu_a_bar*mu_y_bar^(-sigma)))^(-1))^(alpha/(1-alpha))
*h_bar*mu_a_bar^(alpha/((1-alpha)^2))*mu_x_bar^(alpha/(alpha-1));

g_bar=0.2*y_bar/x_g_bar;

k_bar=(alpha*beta*mu_y_bar^(1-sigma)*(1-beta*(1-delta0)*mu_a_bar*mu_y_bar^(-sigma)))^(-1)*y_bar;

i_bar=(1-(1-delta0)/mu_k_bar)*k_bar;

c_bar=y_bar-i_bar-g_bar*x_g_bar;

lam_bar=c_bar^(-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma))
*((1-theta_c/mu_y_bar)^(-sigma)-beta*theta_c*(mu_y_bar-theta_c)^(-sigma));


model;
//equation 1
y = z+alpha*u+alpha*k(-1)+(1-alpha)*h+(alpha/(1-alpha))*mu_a-alpha*mu_x;

//eq 2
c_bar*c + i_bar*i + g_bar*x_g_bar*(g+x_g)=y_bar*y;

//eq 3
k_bar*k=((1-delta0)*k(-1)-(1-delta0)*mu_k-delta1*u)*k_bar/mu_k_bar+i_bar*i
+i_bar*mu_k_bar*kappa*(i(-1)-i-mu_k)/2;

//eq 4
lam_bar*lam=(-sigma)*c_bar^(-sigma)*(1-theta_c/mu_y_bar)^(-sigma-1)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma))
*(c-(theta_c/mu_y_bar)*c(-1)+(theta_c/mu_y_bar)*mu_y)
+chi*(1-sigma)*c_bar^(-sigma)*(1-theta_c/mu_y_bar)^(-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma)-1)*h_bar*(theta_l*h(-1)-h)
-beta*theta_c*(-sigma)*c_bar^(-sigma)*(mu_y_bar-theta_c)^(-sigma-1)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma))*(mu_y_bar*c(+1)-theta_c*c+mu_y_bar*mu_y(+1))
-beta*theta_c*chi*(1-sigma)*c_bar^(-sigma)*(mu_y_bar-theta_c)^(-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma)-1)*h_bar*(theta_l*h-h(+1));

//eq 5
(1-alpha)*(y_bar/h_bar)*lam_bar*(lam+y-h)=chi*(1-sigma)*c_bar^(1-sigma)*(1-theta_c/mu_y_bar)^(-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma)-1)
*(c-(theta_c/mu_y_bar)*c(-1)+(theta_c/mu_y_bar)*mu_y)
+chi*(chi*(1-sigma)-1)*c_bar^(1-sigma)*(1-theta_c/mu_y_bar)^(1-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma)-2)*h_bar*(theta_l*h(-1)-h)
-beta*chi*theta_l*(1-sigma)*c_bar^(1-sigma)*(mu_y_bar-theta_c)^(-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma)-1)*(mu_y_bar*c(+1)-theta_c*c+mu_y_bar*mu_y(+1))
-beta*chi*theta_l*(chi*(1-sigma)-1)*c_bar^(1-sigma)*(mu_y_bar-theta_c)^(1-sigma)*((1-theta_l)*(1-h_bar))^(chi*(1-sigma)-2)*h_bar*(theta_l*h-h(+1));

//eq 6
q_bar*(q+lam)=alpha*beta*mu_y_bar^(-sigma)*k_bar^(alpha-1)*h_bar^(1-alpha)*mu_x_bar^(1-alpha)
*(lam(+1)-sigma*mu_y(+1)+alpha*u(+1)+z(+1)+(alpha-1)*k+(1-alpha)*h(+1)+(1-alpha)*mu_x(+1))
+beta*q_bar*mu_a_bar*(lam(+1)*(1-delta0)-sigma*(1-delta0)*mu_y(+1)+(1-delta0)*q(+1)+(1-delta0)*mu_a(+1)-delta1*u(+1));

//eq 7
q_bar*(delta1*q+delta2*u)=alpha*k_bar^(alpha-1)*h_bar^(1-alpha)*mu_x_bar^(1-alpha)*mu_a_bar^(-1)
*(z+(alpha-1)*u+(alpha-1)*k(-1)+(1-alpha)*h+(1-alpha)*mu_x-mu_a);

//eq 8
lam=q_bar*(q+lam)-kappa*q_bar*mu_k_bar*(mu_k_bar*i+mu_k_bar*mu_k-i(-1))
+kappa*beta*q_bar*mu_a_bar*mu_y_bar^(-sigma)*mu_k_bar^3*(i(+1)+mu_k(+1)-i);

//eq 9
mu_y=(alpha/(alpha-1))*mu_a+mu_x;

//eq 10
mu_k=(1/(alpha-1))*mu_a+mu_x;

//eq 11
x_g=rho_xg*x_g(-1)-mu_y;

//eq 12
z=rho_z*z(-1)+eps_z0+Sz1+Sz2+Sz3;
Sz1 = eps_z1;
Sz2 = eps_z2;
Sz3 = eps_z3;
//Sz4 = eps_z4;

//eq 13
mu_x=rho_x*mu_x(-1)+eps_x0+Sx1+Sx2+Sx3;
Sx1 = eps_x1;
Sx2 = eps_x2;
Sx3 = eps_x3;
//Sx4 = eps_x4;

//eq 14
mu_a=rho_a*mu_a(-1)+eps_a0+Sa1+Sa2+Sa3;
Sa1 = eps_a1;
Sa2 = eps_a2;
Sa3 = eps_a3;
//Sa4 = eps_a4;

//eq 15
g=rho_g*g(-1)+eps_g0+Sg1+Sg2+Sg3;
Sg1 = eps_g1;
Sg2 = eps_g2;
Sg3 = eps_g3;
//Sg4 = eps_g4;

end;

varobs y;

initval;
y = 0;
c = 0;
i = 0;
h = 0;
u = 0;
mu_y = 0;
mu_k = 0;
x_g = 0;
q = 0;
lam = 0;
k = 0;
z = 0;
mu_x = 0;
mu_a = 0;
g = 0;
 
Sz1 =0;
Sz2 =0;
Sz3 =0;
//Sz4 =0;
Sx1 =0;
Sx2 =0;
Sx3 =0;
//Sx4 =0;
Sa1 =0;
Sa2 =0;
Sa3 =0;
//Sa4 =0;
Sg1 =0;
Sg2 =0;
Sg3 =0;
//Sg4 =0;
end;

//steady;
//check;

estimated_params;
theta_c, beta_pdf, 0.5, 0.1; 
theta_l, beta_pdf, 0.5, 0.1;
chi, gamma_pdf, 4, 1;
delta2, uniform_pdf, , ,0.01,10; 
kappa, gamma_pdf, 4, 1; 
rho_z, beta_pdf, 0.7, 0.2; 
rho_g, beta_pdf, 0.7, 0.2; 
rho_xg, beta_pdf, 0.7, 0.2;
rho_x, beta_pdf, 0, 0.1; 
rho_a, beta_pdf, 0.5, 0.1; 

stderr eps_z0, uniform_pdf, , , 0, 8.6603;
stderr eps_x0, uniform_pdf, , , 0, 8.6603;
stderr eps_a0, uniform_pdf, , , 0, 8.6603;
stderr eps_g0, uniform_pdf, , , 0, 8.6603;

stderr eps_z1, uniform_pdf, , , 0, 5;
stderr eps_z2, uniform_pdf, , , 0, 5;
stderr eps_z3, uniform_pdf, , , 0, 5;

stderr eps_x1, uniform_pdf, , , 0, 5;
stderr eps_x2, uniform_pdf, , , 0, 5;
stderr eps_x3, uniform_pdf, , , 0, 5;

stderr eps_a1, uniform_pdf, , , 0, 5;
stderr eps_a2, uniform_pdf, , , 0, 5;
stderr eps_a3, uniform_pdf, , , 0, 5;

stderr eps_g1, uniform_pdf, , , 0, 5;
stderr eps_g2, uniform_pdf, , , 0, 5;
stderr eps_g3, uniform_pdf, , , 0, 5;

end;

estimation(datafile=simuldata_news,nobs=200,first_obs=500,
mh_replic=2000,mh_nblocks=2,mh_drop=0.45,mh_jscale=0.8) y c i h k u q;