var  u, y, ygap, yn, pi, i, rn, n, a, r, ry, iy, piy, m, mr, mry,dy,robs,pinfobs,v; 
varexo e_a e_u e_v; 
parameters cgamma,cpie,cr,ctrend,constepinf,conster,RHO_V,ALFA, BETA, THETA, SIGMA, KAPPA, LAMBDA, PHI, PSI_YA, EPSILON, ETA, RHO, FI_PI, FI_Y, RHO_V, RHO_A, RHO_U, AUX; 

RHO_V = 0.50;
RHO_U = 0.5;
RHO_A = 0.79;
ALFA    =  1/3;           
BETA    =  1.03^(-0.25);   
THETA   =  2/3;           
SIGMA   =  1.58; 
PHI     =  1;
EPSILON =  6;          
ETA     =  4;           
RHO     =  1/BETA-1;
FI_PI   =  1.07;
FI_Y    =  0.10;
LAMBDA  =  ((1-THETA)*((1-BETA*THETA)/THETA)*(1-ALFA)/(1-ALFA+ALFA*EPSILON));
KAPPA   =  LAMBDA*(SIGMA + (PHI+ALFA)/(1-ALFA));
PSI_YA  =  (1+PHI)/(PHI+ALFA+SIGMA*(1-ALFA)); 
AUX     =  (1-ALFA)*(log(EPSILON/(EPSILON-1)-log(1-ALFA)))/(SIGMA*(1-ALFA)+PHI+ALFA);
cgamma=1.004;
cpie=1.005;
cr=cpie/(BETA*cgamma^(-SIGMA));
ctrend=(cgamma-1)*100;
conster=(cr-1)*100;
constepinf=(cpie-1)*100;

model(linear);                              
yn     = PSI_YA*a + AUX; 
ygap   = ygap(+1) - 1/SIGMA*(i-pi(+1)-rn);      
y      = ygap + yn;
r      = i - pi(+1);
pi     = BETA*pi(+1) + KAPPA*ygap+u;             
rn     = RHO + SIGMA*PSI_YA*(a(+1)-a);        
y      = a  + (1-ALFA)*n;                       
m      = pi + y - ETA*i;                         
a      = RHO_A*a(-1)+e_a;                                              
u      = RHO_U*u(-1)+e_u; 
ry     = 4*r;                       
iy     = 4*i;                       
piy    = 4*pi;                      
mr     = m - pi;                   
mry    = 4*mr;                      
i = RHO+FI_PI*pi+FI_Y*ygap+v;
v = RHO_V*v(-1)+e_v;
dy=y-y(-1)+ctrend;
pinfobs = 1*(pi) + constepinf;
robs =    1*(i) + conster;
end;

initval;
ygap  = 0;
y     = 0;
yn    = 0;
m     = 0;
n     = 0;
pi    = 0;
i     = RHO;
rn    = RHO;
a     = 0;
v     = 0;
u     = 0;
e_a   = 0;
e_u   = 0;
e_v   = 0;
r     = RHO;
ry    = 4*r;
iy    = 4*i;
piy   = 4*pi;
end;

shocks;
var e_a; stderr 1.29;
var e_u;stderr 0.25;
var e_v; stderr 0.5;        
end;

steady;
check;

estimated_params;
constepinf,0.7,0.1,2.0,GAMMA_PDF,0.625,0.1;
ctrend,0.3982,0.1,0.8,NORMAL_PDF,0.4,0.10;
stderr e_a,0.4618,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr e_u,0.1455,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr e_v,0.2397,0.01,3,INV_GAMMA_PDF,0.1,2;
RHO_A,.9676 ,.01,.9999,BETA_PDF,0.5,0.20;
RHO_U,.8692,.01,.9999,BETA_PDF,0.5,0.20;
RHO_V,.3,.01,.9999,BETA_PDF,0.5,0.20;
FI_PI,1.7985,1.0,3,NORMAL_PDF,1.5,0.25;
FI_Y,0.0893,0.001,0.5,NORMAL_PDF,0.125,0.05;
THETA,0.7813,0.5,0.95,BETA_PDF,0.5,0.10;
ALFA, 0.24,0.01,1.0,NORMAL_PDF,0.3,0.05;
SIGMA,1.2312,0.25,3,NORMAL_PDF,1.50,0.375;
PHI,2.8401,0.25,10,NORMAL_PDF,2,0.75;
EPSILON, 6.3325,2,15,NORMAL_PDF,4,1.5;
end;

varobs dy robs pinfobs;

estimation(datafile=usdata,mh_nblocks=5,mh_replic=10000,mode_compute=6);

stoch_simul(irf=12, order =1) dy pinfobs robs;