var y 
    r 
    p
    g 
    x 
    a 
    e  
    z;
//*************************************************************************************************    
varexo eps_a eps_e eps_z eps_r;
//*************************************************************************************************

parameters  //parameters to be estimate
            omega 
            alpha_x 
            alpha_p  
            rho_p 
            rho_g 
            rho_x 
            rho_a 
            rho_e 
            sig_a 
            sig_e 
            sig_z 
            sig_r  
            //fixed paramters
            z_ss 
            p_ss 
            betta 
            psi; 
//fixed parameters,not to be estimated
z_ss=1.0048;
p_ss=1.0086;
betta=0.99;
psi=0.1;

/*
The following model equations are the stationary model equations in log-linearized form.
*/

model(linear);
//1. Stochastic process for preference shock(15)
a=rho_a*a(-1)+eps_a;
//2. Stochastic process for markup shock(16)
e=rho_e*e(-1)+eps_e;
//3. Stochastic process for technology shock(17)
z=eps_z;
//4.expectational IS curve(23)
x=alpha_x*x(-1)+(1-alpha_x)*x(+1)-(r-p(+1))+(1-omega)*(1-rho_a)*a;
//5.New Keynesian Pillips curve(24)
p=betta*(alpha_p*p(-1)+(1-alpha_p)*p(+1))+psi*x-e;
//6.output gap
x=y-omega*a;
//7.output growth
g=y-y(-1)+z;
//8.Taylor rule
r-r(-1)=rho_p*p+rho_g*g+rho_x*x+eps_r;
end;
//note:equation 5,7,8 are observation functions.

shocks;
var eps_a; stderr sig_a;
var eps_e; stderr sig_e;
var eps_z; stderr sig_z;
var eps_r; stderr sig_r;
end;  

varobs g p r; 

//For MLE
estimated_params;
// full sample
omega,  0.625,0,inf; 
alpha_x, 0.25,0,1; 
alpha_p, 0.25,0,1;
rho_p, 0.3515;
rho_g, 0.2556;
rho_x, 0.0285;
rho_a, (100*0.0287)^2/(1+(100*0.0287)^2);
rho_e, (100*0.0676)^2/(1+(100*0.0676)^2);
stderr eps_a, 0.0187; 
stderr eps_e, 0.0014;
stderr eps_z, 0.0067;
stderr eps_r,  0.0032;
end;
// post-1980
//omega, sqrt(0.625/0.375); 
//alpha_x, sqrt(0.25/0.75); 
//alpha_p, sqrt(0.25/0.75) ; 
//rho_p, 0.6017;
//rho_g, 0.4240;
//rho_x, 0.0770 ;
//rho_a, 0.0205 ;
//rho_e, 0.0882 ;
//sig_a, 0.0150; 
//sig_e, 0.0008;
//sig_z, 0.0000;
//sig_r, 0.0030 ;

estimated_params_init(use_calibration);
end;     


estimation(optim=('MaxIter',2500),datafile=gpr_data,mode_check,prefilter=0,prior_trunc=0,mode_compute=9,nodisplay) y r p g x a e z;

identification(advanced=0,max_dim_cova_group=3);
stoch_simul(irf=20,periods=200) g p r x;
shock_decomposition;