close all
clc; clf;

// Endogenous variables
var chat lambdahat yhat ahat pihat qhat xhat ehat uhat mhat muhat ghat zhat rhat;

// Exogenous variables
varexo epsilon_a epsilon_z epsilon_u epsilon_e ;

parameters z_ss r_ss gamma alpha psi rho_a rho_u rho_e delta_r phi beta rho_mm rho_mpi rho_mx rho_r rho_pi rho_x ;

z_ss = 1.07;          
r_ss = 0.04;          
gamma = 0.5;
rho_a = 0.75;
alpha = 0.5;         
psi =   0.1;
phi =   10;
rho_u = 0.75;
rho_e = 0.5;
beta =  0.99;
delta_r = 15;
rho_r = 0.75;
rho_pi = 0.4;
rho_x = 0.2;

//rho_mm = 0;
//rho_mpi = 0;
//rho_mx = 0;

model(linear);

chat = yhat;

(z_ss-beta*gamma)*(z_ss-gamma)*lambdahat = gamma*z_ss*yhat(-1) - (z_ss^2+beta*gamma^2)*yhat + beta*gamma*z_ss*yhat(+1) + (z_ss-beta*gamma*rho_a)*(z_ss-gamma)*ahat - gamma*z_ss*zhat ;

lambdahat = rhat + lambdahat(+1) - pihat(+1) ;

0 = gamma*z_ss*qhat(-1) - (z_ss^2 + beta*gamma^2)*qhat + beta*gamma*z_ss*qhat(+1) + beta*gamma*(z_ss-gamma)*(1-rho_a)*ahat - gamma*z_ss*zhat ;

xhat = yhat - qhat ;

(1+beta*alpha)*pihat = alpha*pihat(-1) + beta*pihat(+1) - psi*lambdahat + psi*ahat + ehat ;

//rhat = rho_r*rhat(-1) + rho_pi*pihat(-1) + rho_x*xhat(-1) + epsilon_r ;

delta_r*(r_ss-1)*lambdahat - delta_r*(r_ss-1)*ahat - uhat + phi*zhat = phi*mhat(-1) - (1+(1+beta)*phi)*mhat + beta*phi*mhat(+1) - delta_r*rhat ; 

muhat = zhat + mhat - mhat(-1) + pihat ; 

ghat = yhat - yhat(-1) + zhat; 

ahat = rho_a*ahat(-1) + epsilon_a ; 

uhat = rho_u*uhat(-1) + epsilon_u ; 

ehat = rho_e*ehat(-1) + epsilon_e ; 

zhat = epsilon_z ; 

muhat = muhat(-1) + -1.25*xhat(-1);

end;

check;


//----------------------------------------------------------------
// 5. Shocks 
//----------------------------------------------------------------

shocks;
var epsilon_a; stderr 0.0125; 
var epsilon_z; stderr 0.0125; 
var epsilon_u; stderr 0.0125; 
var epsilon_e; stderr 0.0031; 
end;

steady;

check;

stoch_simul(order=1,irf=20); // order =1 because it is a linear model

//----------------------------------------------------------------
// 6. Estimating of parameters
//----------------------------------------------------------------

estimated_params;

beta,   normal_pdf,     1.1890,     0.0089;     // beta follows a normal distribution with mean = 1.189 and teh standard error = 0.0089
alpha,  beta_pdf,       0.3329,     0.0006;     // price indexation
psi,    gamma_pdf,      0.0315,     0.0015;     // Phillips curve slope
gamma,  beta_pdf,       0.6741,     0.0059;     // habit formation *
delta_r,gamma_pdf,      6,          0.9366;     // annual money demand semi-elasticity * 5.5958
phi,    gamma_pdf,      18,         0.9394;     // money demand adjustment cost .786
rho_r,  normal_pdf,       0.8726,     0.0089;     // interest rate smoothing
rho_pi, normal_pdf,      1,          0.00002;    // Policy response to inflation
rho_x,  normal_pdf,      0.2129,     0.00003;    // Policy response to output gap
rho_a , normal_pdf,       0.98,       0.0081;     // Preference shock persistence
rho_u , normal_pdf,       0.98,       0.0082;     // Money demand shock persistence
rho_e , normal_pdf,       0.3330,     0.0006;     // Cost push shock persistence
r_ss,      normal_pdf,     0.492,      0.0016;   
z_ss,      normal_pdf,     0.6568,     0.0054;

stderr  epsilon_a,     gamma_pdf,     0.0049,      0.0004,   0;      // preference shock volatility
stderr  epsilon_z,     gamma_pdf,     0.0290,      0.0094,   0;      // productivity shock volatility
stderr  epsilon_u,     gamma_pdf,     0.0697,      0.005,    0;      // money demand shock volatility
stderr  epsilon_e,     gamma_pdf,     0.004,      0.0002,   0;      // cost push shock volatility
//stderr  epsilon_r,     gamma_pdf,     0.0035,      0.0004;           // monetary policy shock volatility

end;

//varobs ghat pihat muhat;
estimated_params_init(use_calibration);

//----------------------------------------------------------------
// 7. Estimating NKM 
//----------------------------------------------------------------

estimation(optim=('MaxIter',200),datafile='Dynare_data.xlsx',prior_trunc= 0,mode_compute=6,presample=4,lik_init=2,mode_check,prefilter=0,
mh_replic=2000,mh_nblocks=2,mh_jscale=0.8,conf_sig=.95,mh_drop=0.45,bayesian_irf,irf=20);