//works with dynare 4.4.3 
//TR(3); xi_pi estimated (or not)

var pi y i a e v z mc phi m
    pi_obs y_obs i_obs;
varexo  e_a e_e e_u e_v;			

parameters  sigma eta h beta alpha omega epsilon rho_a rho_e rho_v; 
parameters  phi_y phi_pi phi_delta1 phi_delta2 phi_deltay rho_i pi_ss;
   omega   =   0.5;
   sigma   =   1.5;         
   eta     =   1;  
   h       =   0.8;       
   beta    =   0.99;      
   alpha   =   2/3;           
   epsilon =   5;

    rho_i   =   0.8;
    phi_pi  =   1.5;      
    phi_y   =   0.5/4;  
 phi_delta1 =   0.3;
 phi_delta2 =   0.3;
 phi_deltay =   0.0625;  
    pi_ss   =   1.005;  //non zero inflation, check
    rho_a   =   0.8;       
    rho_e   =   0.8;       
    rho_v   =   0.8;       


//MODEL DECLARATION

model(linear);
#mju = epsilon/(epsilon-1);
#R = pi_ss/beta;
#psi_cap    = alpha*pi_ss^(epsilon-1) + omega*(1-alpha*pi_ss^epsilon*(pi_ss^(-1)-beta));// 1.1750
#xi_f	= alpha*beta*pi_ss^epsilon/psi_cap;                                            //0.5730
#xi_b	= omega/psi_cap;                                                               //0.4255

#kappa = (1-omega)*(1-alpha*pi_ss^(epsilon-1))*(1-alpha*beta*pi_ss^epsilon)/psi_cap;     //0.0449
#xi_pi = alpha*beta*pi_ss^(epsilon-1)*(pi_ss-1)*(1-omega)*(1-alpha*pi_ss^(epsilon-1))/psi_cap; //4.6081e-04

  y = (1/(1+h))*y(+1)+(h/(1+h))*y(-1)-((1-h)/sigma*(1+h))*(i - pi(+1)+(e(+1)-e));
 pi = xi_f*pi(+1) + xi_b*pi(-1) + kappa*(mc+v) + xi_pi * phi;

 mc = (sigma+eta)*y +eta*z - (eta+1)*a;
phi = ((epsilon-1)*pi(+1)+m(+1)-(1-sigma)*y);
  m = (1-alpha*beta*pi_ss^(epsilon-1))*((1-sigma)*y)+(alpha*beta*pi_ss^(epsilon-1))*((epsilon-1)*pi(+1)+m(+1));
       z = epsilon*alpha*(pi_ss^epsilon)*(1-pi_ss^(-1))/(1-alpha*pi_ss^(epsilon-1))*pi+alpha*pi_ss^epsilon*z(-1);

  i = rho_i *i(-1)+ phi_pi * pi + phi_y*y  +phi_delta1*(pi-pi(-1))
      + phi_delta2*(pi(+1)-pi)+phi_deltay*(y - y(-1))+ e_u; 

      a = rho_a * a(-1) + e_a;
      e = rho_e * e(-1) + e_e;
      v = rho_v * v(-1) + e_v;

//observation equations
pi_obs = pi + log(pi_ss);
 y_obs = y-y(-1);
 i_obs = i + log(R);

end;

/*
check;
steady;

shocks;
var  e_a;  stderr  0.1;            
var  e_e;  stderr  0.1;
var  e_u;  stderr  0.1;
var  e_v;  stderr  0.1;
end;

stoch_simul;
*/


varobs pi_obs y_obs i_obs; 

//specify priors
estimated_params;
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF

stderr e_a,     INV_GAMMA_PDF,0.1, 2;      //technology
stderr e_e,     INV_GAMMA_PDF,0.1, 2;      //demand shock stderr
stderr e_u,     INV_GAMMA_PDF,0.1, 2;      //cost push- rescaled stderr
stderr e_v,     INV_GAMMA_PDF,0.1, 2;      //MP shock

sigma,   NORMAL_PDF, 1, 0.375;
eta,    NORMAL_PDF, 2, 1.5;
h,      BETA_PDF,   0.5, 0.15;
alpha,  BETA_PDF,   0.6, 0.15;
epsilon,NORMAL_PDF, 5, 1.5;

rho_a, BETA_PDF,   0.8,0.1;
rho_e, BETA_PDF,   0.8,0.1;
rho_v, BETA_PDF,   0.8,0.1;

rho_i,      BETA_PDF,   0.8, 0.1;
phi_pi,     GAMMA_PDF,  1.5, 0.1, 1;
phi_y,      GAMMA_PDF,  0.5/4, 0.2;
phi_delta1, GAMMA_PDF,  0.8, 0.2;
%phi_delta2, GAMMA_PDF,  0.3, 0.1;
%phi_deltay, GAMMA_PDF,  0.0625, 0.05;  
end;

//specify initial conditions for maximising the posterior
/*
estimated_params_init;
phi_pi, 1.5;
phi_y,  0.5/4;
phi_delta1, 0.3;
rho_i, 0.8;
eta, 1;
end;
*/

//estimation(datafile=EA,mh_replic=2000, mh_nblocks =2, mh_drop=0.45,mh_jscale=0.55,mode_compute=6) ;
estimation(datafile= EA, mh_replic=5000, mh_nblocks =2,mh_drop=0.45,
mh_jscale=0.3,mode_file= EA_str_mode, mode_compute=0, nobs=[86:90], forecast =4) pi y i ;
//stoch_simul (irf=0, nofunctions)  pi  y  i ;

//load_mh_file,