%This is the baseline model for the EPR article

%
%VAR declares the endogenous variables
%
var 

//IS Block
xtil, x, YA,

//PC Block
pitil, pi, 

//Policy Rule
ira,

//Natural Rates
rn_t, lambdaAn, YAn, 

//Efficent Rates
re_t, lambdaAe, YAe, 

//Shocks
delta, gamma_t, u, pistar_t, 

//Definitions
pi4Q, gY4Q, ra_t, rna_t, rngap, xn, rea_t, regap, xe,

//Observables
gYa, pia, ir;

%
%VAREXO declares the shocks
%
varexo edelta, egamma , eu, ei, epistar;

%
%PARAMETERS declares the parameters
%
parameters 

//Tastes and Technology
omega, xi, eta, zeta, 

//Policy
rho, phipi, phix

//Intercepts
pistar, ra, gammaa, 

//Shocks
rhopistar, rhodelta, rhogamma, rhou;

%
%Here is the MODEL
%
model;

//Calibrated parameters
# beta = 0.99;

//Parameter transformations
# gamma = gammaa/400;
# phigammatil = exp(gamma)/(exp(gamma)-beta*eta);
# etagammatil = exp(gamma)/(exp(gamma)-eta);
# phigamma = phigammatil*etagammatil;
# etagamma = eta/exp(gamma);
# r = ra/400;
# pibar=pistar/400;

//Model's equations

//IS Block
xtil = xtil(+1) - phigamma^(-1)*(ir - pi(+1) - rn_t);
xtil = x - (beta*etagamma)*(x(+1) - etagamma*x);
x = YA - YAn;

//PC Block
pitil = beta*pitil(+1) + xi*(omega*x + phigamma*xtil);
pitil = pi - zeta*pi(-1)-(1-zeta)*pibar;

//Policy Rule
ira = rho*ira(-1) + (1-rho)*(rea_t + pistar_t + phipi*(pi4Q - pistar_t) + phix*xe) + ei;

//Natural Rates
lambdaAn = lambdaAn(+1) - gamma_t(+1) + rn_t - r - delta(+1);
lambdaAn = -phigamma*(YAn - etagamma*(YA(-1) - gamma_t) - beta*etagamma*(YAn(+1) + gamma_t(+1) - etagamma*YAn)) + ((beta*eta)/(exp(gamma) - beta*eta))*delta(+1);
YAn = (1/omega)*(lambdaAn-(1/xi)*u);

//Efficient Rates
lambdaAe = lambdaAe(+1) - gamma_t(+1) + re_t - r - delta(+1);
lambdaAe = -phigamma*(YAe - etagamma*(YA(-1) - gamma_t) - beta*etagamma*(YAe(+1) + gamma_t(+1) - etagamma*YAe)) + ((beta*eta)/(exp(gamma) - beta*eta))*delta(+1);
YAe = omega^(-1)*lambdaAe;

//Shocks
delta = rhodelta*delta(-1) + (edelta/400);
gamma_t = rhogamma*gamma_t(-1) + (egamma/400);
u = rhou*u(-1) + (eu/400);
pistar_t = (1-rhopistar)*pistar + rhopistar*pistar_t(-1) + epistar;

//Definitions
pi4Q = 100*(pi + pi(-1) + pi(-2) + pi(-3));
gY4Q = 100*(YA - YA(-4) + gamma_t + gamma_t(-1) + gamma_t(-2) + gamma_t(-3) + 4*gamma);
ra_t = 400*(ir - pi(+1));
rna_t = 400*rn_t;
rngap = ra_t - rna_t;
xn = 100*x;
rea_t = 400*re_t;
regap = ra_t - rea_t;
xe = 100*(YA - YAe);

//Observation Equations
gYa = gammaa + 400*(YA - YA(-1) + gamma_t);
pia = 400*pi;
ira = 400*ir;

end;

%
%Now the ESTIMATION
%

%The observable variables
varobs gYa pia ira;

%The priors on the parameters
estimated_params;

omega , gamma_pdf , 1 , 0.2;
xi , gamma_pdf , 0.1 , 0.05;
eta , beta_pdf , 0.6 , 0.2;
zeta , beta_pdf , 0.6 , 0.2;
rho , beta_pdf , 0.7 , 0.15;
phipi , normal_pdf , 1.5 , 0.25;
phix , normal_pdf , 0.5 , 0.2;
pistar , normal_pdf , 2 , 1;
ra , normal_pdf , 2 , 1;
gammaa , normal_pdf , 3 , .35;
rhopistar , beta_pdf , 0.95 , 0.04;
rhodelta , beta_pdf , 0.5 , 0.2;
rhogamma , beta_pdf , 0.5 , 0.2;
rhou , beta_pdf, 0.5, 0.2;
stderr edelta, inv_gamma1_pdf, 0.5 , 2;
stderr egamma, inv_gamma1_pdf, 0.5 , 2;
stderr eu, inv_gamma1_pdf, 0.5 , 2;
stderr epistar, inv_gamma1_pdf, 0.2 , 1;
stderr ei, inv_gamma1_pdf, 0.5 , 2;

end;

%Set the scale for the graphs (labels every 16 quarters)
options_.XTick = 1:16:92;
options_.XTickLabel = strvcat('84','88','92','96','00','04');

%Run the estimation the first time
estimation(datafile = usedata,
           mh_replic = 100000,
           mh_nblocks = 5,
           mh_drop = 0.2,
           mh_jscale = 0.48,
           nodiagnostic,
           moments_varendo,
           conf_sig = 0.90,
           smoother,
           filtered_vars,
           forecast=8,
           irf=16,
           bayesian_irf);

%Run the estimation when mh_file and/or mode file already exist
% estimation(datafile=usedata,
%            mode_file=model3_apr27_2007_mode,
%            load_mh_file,
%            mh_replic=0,
%            nodiagnostic,
%            moments_varendo,
%            conf_sig = 0.90,
%            smoother,
%            filtered_vars,
%            forecast=8,
%            irf=16,
%            bayesian_irf);

%Save Necessary Structures
%global bayestopt_