// Dynare program for a medium-scale DSGE model
// Model setup is based on Levin, Onatski, Williams and Williams (2005)
// Following them, we normalize some of the shocks
//  so that their amplitude does not
//  depend on other parameters. (e_b, e_i, e_g, ita_p and ita_w).

clc

var c r pi q rk i k u l w y e_b e_i e_g e_a e_l e_u pi_target
 cf rf qf rkf if kf uf lf wf yf ygap ygappro ut wmrs;
varexo nu_b nu_i nu_g nu_a nu_l nu_u ita_w ita_p ita_r nu_pi;

parameters beta, delta, alpha, iy, ky, gy, cy, lambdaw, rk0;
parameters theta sigma kai zetainv psi phi1 xip xiw gammap gammaw lfix;
parameters ri rpi ry rdpi rdy;
parameters sigma_a sigma_pi sigma_b sigma_g sigma_l sigma_i;
parameters  sigma_r sigma_q sigma_p sigma_w sigma_u;
parameters rho_a rho_pi rho_b rho_g rho_l rho_i rho_u;
parameters uadj;

beta = 0.995; //discount rate, Fukunaga (2002)
delta = 0.06;  //depreciation rate, JEM
alpha = 0.37;	//capital share, JEM
cy = 0.6; 	//consumption-output ratio
gy = 0.2;	//government expenditure-output ratio
iy = 1-cy-gy;	//investment-output ratio
ky = iy/delta;	//capital-output ratio
lambdaw = 0.2; //wage markup, Levin et.al. (2005)
rk0=1/beta-(1-delta); //steady state rental rate

theta=0.7;	//habit formation, LOWW(2005)'s prior
sigma=1;	//inverse long-run IES, SW(2003)'s prior
zetainv=4;  //adjustment costs, SW(2003)'s prior
psi=1;	//inverse of the elsticity of utilization cost
phi1=0.075;	//fixed cost share in production, LOWW(2005)
gammaw=0.5;	//wage indexation, LOWW's prior
xiw=0.375;	//Calvo wage no-revise prob, LOWW's prior
kai=2;		//inverse labor supply elastcity, SW(2003)'s prior
gammap=0.5;	//price indexation, LOWW's prior
//gamma=0.04;     //coefficient on MC in inflation dynamics, refer to Altig et.al.
xip=0.375;	//Calvo price no-revise prob, LOWW's prior
lfix=0;		//labor overhead

ri=1;		//policy, lag interest rate
rpi=0.5;	//policy, pi
ry=0.01;	//policy, y
rdpi=0.1;	//policy, diff pi
rdy=0.1;	//policy, diff y

rho_b=0.85;	//persistence of preference, LOWW(2005)
rho_i=0.85;	//persistence of investment, LOWW(2005)
rho_g=0.85;	//persistence of gov expenditure, LOWW(2005)
rho_a=0.85;	//persistence of productivity, LOWW(2005)
rho_l=0.85;	//persistence of labor supply, LOWW(2005)
rho_pi=0.85; 	//persistence of inflation target shock, LOWW(2005)
rho_u=0.85; 	

sigma_a=1;	//SD of productivity shock
sigma_pi=0.1;	//SD of inflation target shock
sigma_b=0.1;	//SD of preference shock
sigma_g=0.4;	//SD of gov. expenditure shock
sigma_l=2;	//SD of labor supply shock
sigma_i=1;	//SD of investment shock
sigma_u=0.5;
sigma_r=0.1;	//SD of monetary policy shock
sigma_q=3;	//SD of equity premium shock
sigma_p=0.2;	//SD of price markup shock
sigma_w=0.2;	//SD of wage markup shock

uadj=1;

model(linear);

//consumption
c = 1/(1+theta+beta*theta^2)*(theta*c(-1)
 +(1+beta*theta^2+beta*theta)*c(+1)-beta*theta*c(+2)
 -(1-theta)/sigma*((1-beta*theta)*(r-pi(+1))))
 +e_b;

//asset price
q = -(r-pi(+1))+1/(1-delta+rk0)*((1-delta)*q(+1)+rk0*rk(+1));

//investment
i = 1/(1+beta)*(i(-1)+beta*i(+1)+1/zetainv*q)
 -e_i;

//capital
k = (1-delta)*k(-1)+delta*i-rk0*(ut+e_u);

//capital utilization
ut = psi*(rk-q) - (1+psi)*e_u;

uadj*u = ut;

//labour supply
(1+lfix)*l = -w+rk+ut+k(-1);

//market clearing condition
y = cy*c+e_g+delta*ky*i;

//production function
y = (1+phi1)*(e_a+alpha*(ut+k(-1))+(1-alpha)*(1+lfix)*l);

//wage
w = 1/(1+beta)*(beta*w(+1)+w(-1)+beta*pi(+1)
 -(1+beta*gammaw)*pi+gammaw*pi(-1)
 -lambdaw*(1-beta*xiw)*(1-xiw)/(lambdaw+(1+lambdaw)*kai)/xiw
 *(w-kai*l-e_l+beta*theta/(1-beta*theta)
 *(1-rho_b)*e_b*(1+theta+beta*theta^2)*sigma
 /(1-rho_b*(1+beta*theta)+rho_b^2*beta*theta)/(1-theta)
 -sigma/(1-theta)/(1-beta*theta)*((1+beta*theta^2)*c
 -theta*c(-1)-beta*theta*c(+1))))
 +ita_w;

//price
pi = 1/(1+beta*gammap)*(beta*pi(+1)+gammap*pi(-1)
 +(1-beta*xip)*(1-xip)/xip*(alpha*rk+(1-alpha)*w-e_a))
 +ita_p;

//policy rule, r is a nominal interest rate
//responding to lagged pi and y
r = ri*r(-1)+(1-ri)*pi_target+rpi*(pi(-1)-pi_target)
 +ry*ygap(-1)+rdpi*(pi-pi(-1))+rdy*(ygap-ygap(-1))+ita_r;
//r = ri*r(-1)+rpi*(pi(-1))+ry*y(-1)+rdpi*(pi-pi(-1))+rdy*(y-y(-1))+ita_r;

e_b = rho_b*e_b(-1)+ nu_b;
e_i = rho_i*e_i(-1)+ nu_i;
e_g = rho_g*e_g(-1)+ nu_g;
e_a = rho_a*e_a(-1)+ nu_a;
e_l = rho_l*e_l(-1)+ nu_l;
e_u = rho_u*e_u(-1)+ nu_u;
pi_target = rho_pi*pi_target(-1)+nu_pi;

//No friction economy
cf = 1/(1+theta+beta*theta^2)*(theta*cf(-1)
 +(1+beta*theta^2+beta*theta)*cf(+1)-beta*theta*cf(+2)
 -(1-theta)/sigma*((1-beta*theta)*rf))
 +e_b;

qf = -rf+1/(1-delta+rk0)*((1-delta)*qf(+1)+rk0*rkf(+1));

if = 1/(1+beta)*(if(-1)+beta*if(+1)+1/zetainv*qf)
 -e_i;

kf = (1-delta)*kf(-1)+delta*if-rk0*(uf+e_u);

uf = psi*(rkf-qf) - (1+psi)*e_u;

(1+lfix)*lf = -wf+rkf+uf+kf(-1);

yf = cy*cf+e_g+delta*ky*if;

yf = (1+phi1)*(e_a+alpha*(uf+kf(-1))+(1-alpha)*(1+lfix)*lf);

0 = wf-kai*lf-e_l+beta*theta/(1-beta*theta)
 *(1-rho_b)*e_b*(1+theta+beta*theta^2)*sigma
 /(1-rho_b*(1+beta*theta)+rho_b^2*beta*theta)/(1-theta)
 -sigma/(1-theta)/(1-beta*theta)*((1+beta*theta^2)*cf
 -theta*cf(-1)-beta*theta*cf(+1));

0 = alpha*rkf+(1-alpha)*wf-e_a;

ygap = y-yf;

ygappro = y-(1+phi1)*e_a;

wmrs=w-kai*l-e_l++beta*theta/(1-beta*theta)
 *(1-rho_b)*e_b*(1+theta+beta*theta^2)*sigma
 /(1-rho_b*(1+beta*theta)+rho_b^2*beta*theta)/(1-theta)
  -sigma/(1-theta)/(1-beta*theta)*((1+beta*theta^2)*c
 -theta*c(-1)-beta*theta*c(+1));

end;

initval;
c=0;
r=0;
pi=0;
q=0;
rk=0;
i=0;
k=0;
u=0;
l=0;
w=0;
y=0;

e_b=0;
e_i=0;
e_g=0;
e_a=0;
e_l=0;
e_u=0;
pi_target=0;

cf=0;
rf=0;
qf=0;
rkf=0;
if=0;
kf=0;
uf=0;
lf=0;
wf=0;
yf=0;

end;

shocks;

var nu_b; stderr sigma_b;
var nu_i; stderr sigma_i;
var nu_g; stderr sigma_g;
var nu_a; stderr sigma_a;
var nu_l; stderr sigma_l;
var nu_pi; stderr sigma_pi;
var nu_u; stderr sigma_u;

var ita_w; stderr sigma_w;
var ita_p; stderr sigma_p;
var ita_r; stderr sigma_r;
end;

//steady;

//stoch_simul c,i,k,u,l,y,w,pi,r;

estimated_params; //even S.D. affects a posterior mode
theta,  beta_pdf,0.7,0.15;	//LOWW(2005)'s prior
sigma, normal_pdf, 1, 0.375;	//SW(2003)'s prior
kai, normal_pdf, 2, 0.25;	//SW(2003)'s prior
zetainv, normal_pdf, 4, 1.5; //SW(2003)'s prior, defined non-inverse in SW and OW
psi, normal_pdf, 1, 1;	
phi1, gamma_pdf, 0.075, 0.0125;	//LOWW
xip, beta_pdf, 0.375, 0.1;	//LOWW(2005)
xiw, beta_pdf, 0.375, 0.1;	//LOWW 
gammap, beta_pdf, 0.5, 0.25;	//LOWW
gammaw, beta_pdf, 0.5, 0.25;	//LOWW

ri, normal_pdf, 1, 0.15;     
rpi, normal_pdf, 0.5, 0.2;     
ry, normal_pdf, 0.01, 0.1;    
rdpi, normal_pdf, 0.1, 0.1;    
rdy, normal_pdf, 0.1, 0.5;     

stderr nu_a, inv_gamma_pdf,1,inf;
stderr nu_pi, inv_gamma_pdf,0.1,inf; 
stderr nu_b, inv_gamma_pdf,0.1,inf;
stderr nu_g, inv_gamma_pdf,0.4,inf;
stderr nu_l, inv_gamma_pdf,2,inf;
stderr nu_i, inv_gamma_pdf,1,inf;
stderr nu_u, inv_gamma_pdf,0.5,inf;
stderr ita_r, inv_gamma_pdf,0.1,inf;
stderr ita_p, inv_gamma_pdf,0.2,inf;
stderr ita_w, inv_gamma_pdf,0.2,inf;

rho_a, beta_pdf, 0.85, 0.1;		
rho_pi, beta_pdf, 0.85, 0.1;	
rho_b, beta_pdf, 0.85, 0.1;	
rho_g, beta_pdf, 0.85, 0.1;	
rho_l, beta_pdf, 0.85, 0.1;	
rho_i, beta_pdf, 0.85, 0.1;
rho_u, beta_pdf, 0.85, 0.1;	

end;

//estimated_params_bounds;
//phi1,0,;
//psi,0,;
//end;

varobs y c i w l u pi r;

//estimation(datafile=fsdat, nobs=60, linear, mode_compute=0, mode_file=mode, mh_replic=0, mh_nblocks=1);
estimation(datafile=fsdat, nobs=60, mh_replic=1000, mh_nblocks=1);