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 = 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,0.99,0.02;
alpha, normal_pdf,0.65,0.02; // price indexation
psi, normal_pdf, 0.7, 0.02; // Phillips curve slope
gamma,beta_pdf, 0.5, 0.02; // habit formation
delta_r,normal_pdf, 12, 0.05; // annual money demand semi-elasticity
phi,normal_pdf, 7.5, 0.05; // money demand adjustment cost
rho_r, normal_pdf, 1.1, 0.02 ; // interest rate smoothing
rho_pi, normal_pdf, 0.7, 0.02 ; // policy response to inflation
rho_x, normal_pdf, 0.15, 0.02; // policy response to output gap
rho_a ,normal_pdf,0.80,0.05; // preference shock persistence
rho_u ,normal_pdf,0.95,0.02; // money demand shock persistence
rho_e ,normal_pdf,0.85,0.02; // cost push shock persistence
r_ss,  normal_pdf, 0.85, 0.05 ;
z_ss,  normal_pdf, 1.1, 0.05 ;
stderr epsilon_a,gamma_pdf, 0.45,0.02; // preference shock volatility
stderr epsilon_z,gamma_pdf,1.8,0.02; // productivity shock volatility
stderr epsilon_u,gamma_pdf, 0.3,0.002; // money demand shock volatility
stderr epsilon_e,gamma_pdf, 0.8,0.05; // cost push shock volatility
//stderr epsilon_r,normal_pdf, 0.015,0.04; // monetary policy shock volatility


end;

varobs ghat pihat muhat;
estimated_params_init(use_calibration);

end;

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

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