// as in aus_soe_pibar_test1.mod
// but here I increase mh_replic to 10,000 and use the previously calculated
// model from pibar_test1.mod






% psi here is psif in JP's paper, the deviation from the law of one price. 
% (I'm using the same notation as in Dennis et al).
% Likewise, JP assume that the shock to pif is iid. Here, as in JP, I allow 
% an AR(1) process.

% uq(t) here, as in Dennis is the risk premium shock (AR(1)), whereas in 
% JP it's phi(t)

%dennis_soe.mod
% 30/5/2011

% Model is from Dennis, Leitemo and Soderstrom (mprs and robustness in soe)

% in this first model, the CB only cares about inflation and output.
% The simple Taylor rule only responds to inflation (annual) and output
%----------------------------------------------------------------------
%List of variables:
var c ug q e pstar p r rstar a uq y pf pid mu ua s psi upif ystar 
pistar pi pif pd dq ds de dr pibar rtilda Y PIE R YSTAR PISTAR RSTAR
DQ DS PIF;


%c r pi ug q e pstar p rstar pistar dq a uq  y pid mu ua s pf pd
%pif psi upif ds  ystar pibar dr de;


%---------------------------------------------------------------------
%list of exogenous variables:
varexo epsilon_g epsilon_q epsilon_pid epsilon_a epsilon_pif epsilon_pistar
epsilon_ystar epsilon_rstar epsilon_r;% 9 SHOCKS

%------------------------------------------------------------------
% Parameter and values:
parameters alpha beta chi sigma phi eta h thetad thetaf deltad deltaf rhoa rhog rhopif
rhoq  b111 b112 b113 b121 b122 b123 b131 b132 b133 b211 b212 b213 b221 b222 b223 b231
b232 b233 psi1 psi2 lambdaa nu nus rhor;%sigmaa sigmag sigmapd sigmapif sigmaq sigma_r

lambdaa = 0.1; % weight on output in the loss function
nu = 0.8;% weight on interest rate volatility in the loss function
nus = 0;% weight on nom. exchange rate stabilisation


alpha = 0.185;
beta = 0.99;
chi = 0.01;
sigma = 1.309;
phi = 1.1157;
eta = 0.5824;
h = 0;
thetad = 0.7935;
thetaf = 0.5511;
deltad = 0.0499;
deltaf = 0.0693;
rhoa = 0.6936;
rhog = 0.9257;
rhopif = 0.9352;
rhoq = 0.9384;
rhor = 0.5;

%sigmaa = 0.3665;
%sigmag = 0.1610;
%sigmapid = 0.769;
%sigmapif = 1.5769;
%sigmaq = 0.35;

b111 = 0.1;
b112 = 0.1;
b113 = 0.1;
b121 = 0.1;
b122 = 0.1;
b123 = 0.1;
b131 = 0.1;
b132 = 0.1;
b133 = 0.1;

b211 = 0.1;
b212 = 0.1;
b213 = 0.1;
b221 = 0.1;
b222 = 0.1;
b223 = 0.1;
b231 = 0.1;
b232 = 0.1;
b233 = 0.1;


%sigma_pistar = 0.3498;
%sigma_ystar = 0.4795;
%sigma_rstar = 0.11541;

psi1 = 1.5;
psi2 = 0.5;
%-----------------------------------------------------------------------
model(linear);
%1
(1+h)*c-h*c(-1) = c(+1)-(1-h)/sigma*(r-pi(+1)-ug+ug(+1));

%2
ug = rhog*ug(-1)+epsilon_g;

%3
q = e+pstar -p;

%4
(r-pi(+1))-(rstar-pistar(+1))=dq(+1)-chi*a-uq;

%5
uq = rhoq*uq(-1)+epsilon_q;

%6
a = a(-1)/beta+y-c-alpha*(e+pstar-pf);

%7
pid = beta/(1+beta*deltad)*pid(+1)+deltad/(1+beta*deltad)*pid(-1)+
(1-thetad)*(1-beta*thetad)/(thetad*(1+beta*deltad))*mu+epsilon_pid;

%8
mu = phi*y-(1+phi)*ua+alpha*s+sigma/(1-h)*(c-h*c(-1));

%9
ua = rhoa*ua(-1)+epsilon_a;

%10
s = pf-pd;

%11
pif = beta/(1+beta*deltaf)*pif(+1)+deltaf/(1+beta*deltaf)*pif(-1)+
(1-thetaf)*(1-beta*thetaf)/(thetaf*(1+beta*deltaf))*psi+upif;

%12
psi = e+pstar-pf;

%13
upif = rhopif*upif(-1)+epsilon_pif; 

%14
pi = pid+alpha*ds;

%15
y = (1-alpha)*c+alpha*eta*(2-alpha)*s+ alpha*eta*psi+alpha*ystar;

%16

pistar = b111*pistar(-1)+b112*ystar(-1)+b113*rstar(-1)+
    b211*pistar(-2)+b212*ystar(-2)+b213*rstar(-2)+epsilon_pistar;

%17

ystar = b121*pistar(-1)+b122*ystar(-1)+b123*rstar(-1)+b221*pistar(-2)+
b222*ystar(-2)+b223*rstar(-2)+epsilon_ystar;

%18
rstar = b131*pistar(-1)+b132*ystar(-1)+b133*rstar(-1)+b231*pistar(-2)+
b232*ystar(-2)+b233*rstar(-2)+epsilon_rstar;

%19 See below
%r = rhor*r(-1)+(1-rhor)*(psi1*pibar+psi2*y+epsilon_r);


% up to now 19 equations (same as in paper, exc. that they add a redundant one)
% now need to add auxiliary equations

%20
pistar = pstar-pstar(-1);
%21
pi = p-p(-1);

%22
pif = pf-pf(-1);

%23
pid = pd-pd(-1);

%24
dq = q-q(-1);

%25
ds = s-s(-1);

%26
de = e-e(-1);

%27
dr = r-r(-1);

%28
pibar = pi+pi(-1)+pi(-2)+pi(-3);

rtilda = 4*r;

% mpr
%R = (1-rhor)*(psi1*PIE+psi2/4*y)+rhor*R(-1)+epsilon_r;
R = (1-rhor)*(psi1*pibar+psi2/4*y)+rhor*R(-1)+epsilon_r;


% OBSERVED VARIABLES
y = Y;
4*pi = PIE;
4*r = R;
dq = DQ;
ds = DS;
4*pif = PIF;
ystar = YSTAR;
4*rstar = RSTAR;
4*pistar = PISTAR;
 
end;

%--------------------------------------------------------------------------------------

shocks;
var epsilon_g;
stderr 0.161;%sigmag;

var epsilon_q;
stderr 0.35;%sigmaq;

var epsilon_pid;
stderr 0.769; %sigmapid;

var epsilon_a;
stderr 0.3665;%sigmaa;

var epsilon_pif;
stderr 1.5769;%sigmapif;

var epsilon_pistar;
stderr 0.3498;%sigma_pistar;

var epsilon_ystar;
stderr 0.4795;%sigma_ystar;

var epsilon_rstar;
stderr 0.11541;%sigma_rstar;

var epsilon_r;
stderr 0.15;%sigma_r;

end;



estimated_params;
//alpha calibrated for each economy
//beta calibrated
//chi calibrated
sigma, gamma_pdf, 1.2, 0.4;
phi, gamma_pdf, 1.5, 0.75;
thetad, beta_pdf, 0.5, 0.1;
thetaf, beta_pdf, 0.5, 0.1;
eta, gamma_pdf, 1.5, 0.75;
h, beta_pdf, 0.5, 0.25;
deltad, beta_pdf, 0.5, 0.25;
deltaf, beta_pdf, 0.5, 0.25;
rhoa, beta_pdf, 0.8, 0.1;
rhog, beta_pdf, 0.8, 0.1;
rhopif, beta_pdf, 0.8, 0.1;
rhoq, beta_pdf, 0.8, 0.1; % risk premium shock
rhor, beta_pdf, 0.5, 0.1; 

b111, normal_pdf, 0.59, 0.2;
b122, normal_pdf, 0.9, 0.1;
b133, normal_pdf, 0.9, 0.1;

b112, normal_pdf, 0, 0.3;
b113, normal_pdf, 0, 0.3;
b121, normal_pdf, 0, 0.3;
b123, normal_pdf, 0, 0.3;
b131, normal_pdf, 0, 0.3;
b132, normal_pdf, 0, 0.3;

b211, normal_pdf, 0, 0.25;
b212, normal_pdf, 0, 0.15;
b213, normal_pdf, 0, 0.15;
b221, normal_pdf, 0, 0.15;
b222, normal_pdf, 0, 0.25;
b223, normal_pdf, 0, 0.15;
b231, normal_pdf, 0, 0.15;
b232, normal_pdf, 0, 0.15;
b233, normal_pdf, 0, 0.25;



stderr epsilon_a,  INV_GAMMA_PDF, 0.5, inf;% second number is s.e.
stderr epsilon_g, INV_GAMMA_PDF, 0.5, 2;
stderr epsilon_pid, INV_GAMMA_PDF, 0.5, 2;
stderr epsilon_pif, INV_GAMMA_PDF, 0.5, inf;
stderr epsilon_q, INV_GAMMA_PDF, 0.5, inf;

stderr epsilon_pistar, INV_GAMMA_PDF, 0.5, inf;
stderr epsilon_ystar, INV_GAMMA_PDF, 0.5, inf;
stderr epsilon_rstar, INV_GAMMA_PDF, 0.5, inf;

stderr epsilon_r, INV_GAMMA_PDF, 0.5, inf;

//psi1, normal_pdf, 2.0, 0.5; %Joe uses normal_pdf but this gives odd values
psi1, gamma_pdf, 1.5, 0.3; % as in JP
psi2, gamma_pdf, 0.5, 0.13;







end;



%---------------------------------------------------------------------------------------

%stoch_simul(nograph);
%---------------------------------------------------------------------------------------
// US data 1983q1 to 2002q4
% setting lik_init = 2 (instead of 1) if some variables are I(1)
varobs Y PIE R RSTAR DQ DS PIF  YSTAR PISTAR;
%estimation(datafile= australiandata1990q1to2007q2b, first_obs=1, mode_compute=4,lik_init=2,
%mh_nblocks = 5, mh_replic=2000, mh_drop=0.3, mh_jscale=0.4, conf_sig=0.9, nograph, mode_check);

estimation(datafile= aus_dynaredata, first_obs=1, mode_compute=6,
mh_nblocks = 2, mh_replic= 30000, mh_drop=0.3, mh_jscale=0.4, conf_sig=0.9, lik_init = 2, mode_check);