//Broken mod file

//Replication of Schitt-Grohe and Uribe 2012 Econometrica "Whats News in Business Cycles"


//***KNOWN PROBLEMS: 
//h is not properly assigned in code
//observables are probably not properly equated to variables in model
//May only need ME on Y see page 18
//value of xg different than in paper
//initial values are set ad-hoc 
//priors do not specify an upper or lower bound

//*****Changes from previous code
//Changed timing on xg
//MUW calibrated to 1.15 instead of 0.15
//MUX calibrated
//DELTA1 SET TO 0.4 arbirary because prior definition resulted in a negative number
//DELTA1=1/(BETA*MUA*MUY^(-SIGMA))-1-DELTA0;

 
var G C f_delta f_v gc gg gh giv gpa gtfp gy gov I K la mua muk muw mux mu_y xg xi PI q s stock u w Y z z_I zeta
epz4 epz8 epmua4 epmua8 epg4 epg8 epmux4 epmux8 epmuw4 epmuw8 epzeta4 epzeta8 epz_I4 epz_I8;
external_function(name=diff);

varobs gy gc giv gh gg gtfp gpa;

varexo MEY E_z0 E_z4 E_z8 E_z_I0 E_z_I4 E_z_I8 
E_g0 E_g4 E_g8 E_mua0 E_mua4 E_mua8 E_mux0 E_mux4 E_mux8 E_muw0 E_muw4 E_muw8 
E_zeta0 E_zeta4 E_zeta8;

parameters BETA SIGMA ALPHA_K ALPHA_L DELTA0 DELTA1 DELTA2 GAMA KAPA MUW MUX MUA PSI RHOg RHOmuw RHOmua RHOmux RHOXG RHOz RHOzeta RHOz_I THETA b hrs;
 
//Discount Factor
BETA=0.99;

//Intertemporal elasticity of substitution
SIGMA=1;

//Steady State gross per capita GDP growth rate
MUY=1.0045;

//Steady State gross growth rate of price of investment
MUA=0.9957;

//Capital Share
ALPHA_K=0.3;

//Labor Share
ALPHA_L=0.6;

//Steady State Depreciation Rate
DELTA0=0.025;
DELTA1=0.5;
DELTA2=0;

MUX=MUY*MUA^(ALPHA_K/(1-ALPHA_K));

//Steady State hours
hrs=0.2;

//Steady State Wage Markup per author
MUW=1.15;

PSI=1;


//Many of these equations are in non-linear stationary form and taken from Matlab code that is posted on Econometrica

model;

//Steady State share of covernment consumption in GDP
G=0.2*Y/xg;

//Production technology
Y =z * (u*K(-1)/muk)^ALPHA_K * hrs^(1-ALPHA_K-ALPHA_L); 

// Depreciation rate
f_delta = DELTA0 + DELTA1 * (u-1) + DELTA2/2*(u-1)^2; 
 
// Investment adjustment cost
s = KAPA/2 * (xi - muk)^2;

//Capital-trend gross growth rate
muk = mua^(1/(ALPHA_K-1))* mux;

//GDP-trend gross growth rate
mu_y = mua^(ALPHA_K/(ALPHA_K-1))*mux;

// Evolution of investment 
xi = I * muk / I(-1);

//Evolution of capital 
K =(1-f_delta)*K(-1)/muk + z_I *  I * (1-s); 

//Evolution of S (stock) eq (3) in text
stock = (C-b*C(-1)/mu_y)^GAMA*(stock(-1)/mu_y)^(1-GAMA);

//Law of motion of trend in government spending
xg = xg(-1)^RHOXG / mu_y;

// wage rate
w =  z * (1-ALPHA_K-ALPHA_L) * (u*K(-1)/muk/hrs)^ALPHA_K * hrs^(-ALPHA_L); 
 
//Output=Resource constraint
z * (u*K(-1)/muk)^ALPHA_K * hrs^(1-ALPHA_K-ALPHA_L)= C + I + gov * xg;

//argument of period utility function
f_v = (C-b*C(-1)/mu_y) - PSI*hrs^THETA* stock; 


//FOC

//***********************************************
//****Optimality conditions w.r.t. consumption 
//NOTE:(substituted into lhs of optimality condition w.r.t. leisure)

// Verision #1 FOC wrt C - in the paper is currently used
//la=((zeta*f_v^(-SIGMA))-PI*GAMA*stock(+1)/(C-b*C(-1)/mu_y)-BETA*b*(zeta(+1)*f_v(+1)^(-SIGMA)-PI(+1)*GAMA*(stock(+1)/(C(+1)-b*C/mu_y(+1)))));

//Version #2 FOC wrt C - in code NOT USED HERE
//la=((zeta*f_v^(-SIGMA))-PI*GAMA*stock(+1)/(C-b*C(-1)/mu_y)-BETA*b*mu_y(+1)^(-SIGMA)*(zeta(+1)*f_v(+1)^(-SIGMA)-PI(+1)*GAMA*(stock(+1)/mu_y(+1)/(C(+1)-b*C/mu_y(+1)))^(1-GAMA)));
//***********************************************


//************************************************
//*****Optimality conditions w.r.t. leisure
//NOTE: difference between Version 1 and 2 is the time subscript on stock

//Verision #1 FOC wrt leisure - in the paper 
((zeta*f_v^(-SIGMA))-PI*GAMA*stock(+1)/(C-b*C(-1)/mu_y)-BETA*b*(zeta(+1)*f_v(+1)^(-SIGMA)-PI(+1)*GAMA*(stock(+1)/(C(+1)-b*C/mu_y(+1))))) * w / muw = THETA * PSI * zeta * f_v^(-SIGMA) * hrs^(THETA-1) * stock;
 
//Verision #1 FOC wrt leisure - in the matlab code NOT USED HERE
//((zeta*f_v^(-SIGMA))-PI*GAMA*stock(+1)/(C-b*C(-1)/mu_y)-BETA*b*(zeta(+1)*f_v(+1)^(-SIGMA)-PI(+1)*GAMA*(stock(+1)/(C(+1)-b*C/mu_y(+1))))) * w / muw = THETA * PSI * zeta * f_v^(-SIGMA) * hrs^(THETA-1) * stock(+1);
//************************************************


//*************************************************
//Evolution of the shadow price of S
//NOTE: This equation might not be correctly specified
PI = PSI * zeta * f_v^(-SIGMA) * hrs^THETA + BETA * (1-GAMA) * PI(+1) * mu_y(+1)^(1-SIGMA) * ((C(+1)-b*C/mu_y(+1)) /stock(+1))^GAMA * mu_y(+1)^(GAMA-1);
//*************************************************


//Euler equation for capital
q*la = BETA * mua(+1) * mu_y(+1)^(-SIGMA) * la(+1) * (z(+1) * u(+1) * ALPHA_K * (u(+1)*K/muk(+1))^(ALPHA_K-1) * hrs(+1)^(1-ALPHA_K-ALPHA_L) + q(+1) * (1-f_delta(+1)));
 
//w.r.t. u   
q * (DELTA1+DELTA2*(u-1)) = z*ALPHA_K*(u(+1)*K/muk(+1))^(ALPHA_K-1) * hrs(+1)^(1-ALPHA_K-ALPHA_L);


//**************************************************** 
//w.r.t. investment
//Simple omission (mua(+1)*mu_y(+1)^(-SIGMA) ) on my part might not be right

//Version #1 :My verision
la = q*la * z_I * (1-s-xi*KAPA*(xi-muk)) + BETA * q(+1) * la(+1) * z_I(+1) *  xi(+1)^2 * KAPA*(xi(+1)-muk(+1));

//Verision #2 Authors Verision
//la = q*la * z_I * (1-s-xi*KAPA*(xi-muk)) + BETA * mua(+1)*mu_y(+1)^(-SIGMA) * q(+1) * la(+1) * z_I(+1) *  xi(+1)^2 * KAPA*(xi(+1)-muk(+1));
//****************************************************


//evolution of exogenous shocks

z= RHOz * z(-1)+ E_z0 + epz4(-4) + epz8(-8);
epz4=E_z4;
epz8=E_z8; 

mua/MUA = RHOmua * mua(-1)/MUA +  E_mua0 + epmua4(-4) + epmua8(-8);
epmua4=E_mua4;
epmua8=E_mua8;
 
gov/G = RHOg * gov(-1)/G +  E_g0 + epg4(-4) + epg8(-8);
epg4=E_g4;
epg8=E_g8;
 
mux/MUX  = RHOmux * mux(-1)/MUX +  E_mux0 + epmux4(-4) + epmux8(-8);
epmux4=E_mux4;
epmux8=E_mux8;

 
muw/MUW = RHOmuw * muw(-1)/MUW +  E_muw0 + epmuw4(-4) + epmuw8(-8);
epmuw4=E_muw4;
epmuw8=E_muw8;
 
zeta = RHOzeta * zeta(-1) +  E_zeta0 + epzeta4(-4) + epzeta8(-8);
epzeta4=E_zeta4;
epzeta8=E_zeta8;
 
z_I = RHOz_I*z_I(-1) + E_z_I0 + epz_I4(-4) + epz_I8(-8);
epz_I4=E_z_I4;
epz_I8=E_z_I8;


//Observables gy gc giv gh gtfp gpa are equal to diff(log(x)) in the data

gy=Y-Y(-1)+mu_y+MEY;
gc=C-C(-1)+mu_y;
giv=I-I(-1)+mu_y;
gg=gov-gov(-1)+xg^(RHOXG-1);
gh=hrs-hrs(-1);
gtfp =z-z(-1)+(1-ALPHA_K)+mux;
gpa=mua-mua(-1);
end;



estimated_params;
THETA,gamma_pdf,3.77,3,1;
GAMA, beta_pdf,0.73,0.25,0.00001,0.99999;
KAPA,gamma_pdf,3.92,2,0.0001;
b,beta_pdf,0.909,0.02,0,0.99999;
RHOXG,beta_pdf,0.73,0.25,0.00001,0.99999;


//Measurement Error
stderr MEY,gamma_pdf,1.04,0.5;


//Stationary exogenous shocks
RHOz,beta_pdf,0.91234, 0.25,0.00001,0.99999;
stderr E_z0,gamma_pdf,1.04,0.5;
stderr E_z4,gamma_pdf,0.43,0.5;
stderr E_z8,gamma_pdf,0.43,0.5;

RHOz_I,beta_pdf, 0.73, 0.25,0.00001,0.99999;
stderr E_z_I0,gamma_pdf,1.04,0.5;
stderr E_z_I4,gamma_pdf,0.43,0.5;
stderr E_z_I8,gamma_pdf,0.43,0.5;

RHOg,beta_pdf, 0.73, 0.25,0.00001,0.99999;
stderr E_g0,gamma_pdf,1.04,0.5;
stderr E_g4,gamma_pdf,0.43,0.5;
stderr E_g8,gamma_pdf,0.43,0.5;


//Non-stationary exogenous shocks
RHOmua,beta_pdf, 0.5, 0.25,0.00001,0.99999;
stderr E_mua0,gamma_pdf,1.04,0.5;
stderr E_mua4,gamma_pdf,0.43,0.5;
stderr E_mua8,gamma_pdf,0.43,0.5;

RHOmux,beta_pdf, 0.73, 0.25,0.00001,0.99999;
stderr E_mux0,gamma_pdf,1.04,0.5;
stderr E_mux4,gamma_pdf,0.43,0.5;
stderr E_mux8,gamma_pdf,0.43,0.5;

RHOmuw,beta_pdf, 0.73, 0.25,0.00001,0.99999;
stderr E_muw0,gamma_pdf,1.04,0.5;
stderr E_muw4,gamma_pdf,0.43,0.5;
stderr E_muw8,gamma_pdf,0.43,0.5;

RHOzeta,beta_pdf, 0.5, 0.25,0.00001,0.99999;
stderr E_zeta0,gamma_pdf,1.04,0.5;
stderr E_zeta4,gamma_pdf,0.43,0.5;
stderr E_zeta8,gamma_pdf,0.43,0.5;
end; 

estimation(nograph,datafile=Economet,mh_replic=2000,mh_nblocks=1,mh_drop=0.45,mh_jscale=0.5);