clc
clear
close all
addpath(genpath('./functions'));

% 197210 - 202203, 594 obs
GMdata = loadData(197210, 202203);
infl = GMdata{:,'Inflation'};

% Targeting Rate for Germany 
% because it was one of the first countries (along with Switzerland)
% to implement many of the features of an inflation-targeting regime, even 
% though Germany is not an inflation targeter per se.
% Before 1999, we adopt the target rate in Table 3 of Bernanke and Mihov (NBER 1996)
% ECB set the target rate of inflation around 2% since 1999
% 1975-1976   4.5%
% 1977        4%
% 1978        3.5%
% 1979        3%
% 1980        4%
% 1981        3.75%
% 1982-1983   3.5%
% 1984        3%
% 1985        2.5
% 1986-now    2%
TR = nan(size(GMdata,1),1);
ind1 = find(GMdata.yyyymm == 197501);
ind2 = find(GMdata.yyyymm == 197701);
ind3 = find(GMdata.yyyymm == 197801);
ind4 = find(GMdata.yyyymm == 197901);
ind5 = find(GMdata.yyyymm == 198001);
ind6 = find(GMdata.yyyymm == 198101);
ind7 = find(GMdata.yyyymm == 198201);
ind8 = find(GMdata.yyyymm == 198401);
ind9 = find(GMdata.yyyymm == 198501);
ind10 = find(GMdata.yyyymm == 198601);
TR(ind1:ind2-1) = 4.5;
TR(ind2:ind3-1) = 4;
TR(ind3:ind4-1) = 3.5;
TR(ind4:ind5-1) = 3;
TR(ind5:ind6-1) = 4;
TR(ind6:ind7-1) = 3.75;
TR(ind7:ind8-1) = 3.5;
TR(ind8:ind9-1) = 3;
TR(ind9:ind10-1) = 2.5;
TR(ind10:end) = 2;
GMdata.TR = TR;

% Effective Lower Bound
% For GM, we adopt the ECB deposit rate as the effective lower bound (ELB)
% Lemke and Vladu (2016) estimates rlb=0.1% before 201406
% Acoording to deposit rate, 
% ECB rate = -0.1 201406-201408
% ECB rate = -0.2 201409-201511
% ECB rate = -0.3 201512-201602
% ECB rate = -0.4 201603-201908
% ECB rate = -0.5 201909-end
% we choose the minimal of ECB rate and actual Forward rate as the lower bound 
rlb_ind1 = find(GMdata.yyyymm == 201406,1);
rlb_ind2 = find(GMdata.yyyymm == 201409,1);
rlb_ind3 = find(GMdata.yyyymm == 201512,1);
rlb_ind4 = find(GMdata.yyyymm == 201603,1);
rlb_ind5 = find(GMdata.yyyymm == 201909,1);
rb1 = min(-0.1, min(GMdata{rlb_ind1:rlb_ind2-1,3:14},[],'all'));
rb2 = min(-0.2, min(GMdata{rlb_ind2:rlb_ind3-1,3:14},[],'all'));
rb3 = min(-0.3, min(GMdata{rlb_ind3:rlb_ind4-1,3:14},[],'all'));
rb4 = min(-0.4, min(GMdata{rlb_ind4:rlb_ind5-1,3:14},[],'all'));
rb5 = min(-0.5, min(GMdata{rlb_ind5:end,3:14},[],'all'));
% Set lower bound
GMdata.rlb = zeros(size(GMdata,1),1);
GMdata.rlb(rlb_ind1:rlb_ind2-1) = rb1;
GMdata.rlb(rlb_ind2:rlb_ind3-1) = rb2;
GMdata.rlb(rlb_ind3:rlb_ind4-1) = rb3;
GMdata.rlb(rlb_ind4:rlb_ind5-1) = rb4;
GMdata.rlb(rlb_ind5:end) = rb5;
rlb = GMdata{:,'rlb'};

dt = 1;
N =3;
Y = GMdata{:,3:14};
nobs = size(Y,1);
mats = [3, 6, 12:12:120]';

load theta3.mat
load theta0.mat
pars_mode = theta2pars(theta0);
pars = theta2pars(theta3);
loads0 = SRTSM_loadings(pars.kinfQ, diag(pars.lamQ), pars.omega, mats, dt) ;
pars.aJ = loads0.an;
pars.bJ = loads0.bn;
pars.sigmasJ = loads0.SigmaQ;

pars.vt = 0.3;
pars.vc = 0.02;
pars.phi1 = 1.4;
pars.phi2 = -0.9;
pars.vptr = 0.1;


%%%%%======================================================================
%%%%%============================ PRIORS  =================================
%%%%%======================================================================

%%%% sigma_tau2 ~ IG(alpha_eta/2, delta_eta/2)
prior.alpha_eta = 100;  % tightness. Del Negro et al use 100
prior.mean_sigtau2  = pars_mode.sigtau2;   
prior.delta_eta = prior.mean_sigtau2 * (prior.alpha_eta - 2) ;

%%%% Sigma_tilde ~ IW(kappa, Psi) prior mean = Psi/(kappa-N-1) = Psi
Psi_data = diag(diag(pars_mode.Sigma_tilde));
prior.kappa = N+2 ;    % low precision
prior.Psi = Psi_data ;

%%%% vec(Phi) ~ N(0; Omega)
phi_sigma = [0.8875 3.1485 2.2169]';
prior.Phi = pars_mode.Phi;
omega = 0.2*eye(N);
for i = 1:N
    for j = setdiff(1:N,i)
        omega(i,j) = 0.2*phi_sigma(i)/phi_sigma(j);
    end
end
omega = reshape(omega.^2,[],1);
prior.Omega = diag(omega) ;
prior.Omega_inv = diag(1./omega) ;

%%%% ap ~ N(ap_mean, ap_cov)
prior.ap_mean = [pars_mode.a; pars_mode.p];
prior.ap_cov = diag([0.1; 0.1; 1; 1].^2) ;
prior.ap_cov_inv = prior.ap_cov\eye(4);

%%%% sige2 = meas. error variance ~ IG(alpha_e/2, delta_e/2)
prior.alpha_e = 4;
prior.mean_sige2 = pars_mode.sige2;
prior.delta_e = prior.mean_sige2 * (prior.alpha_e - 2) ;

% kinfQ ~ N(delta_k,var_k)
prior.kmean = pars_mode.kinfQ;
prior.kvar = 0.02^2;

% Priors for inflation decomposition
prior.alpha_e1 = 4;
prior.delta_e1 = 1*(prior.alpha_e1-2);
prior.alpha_e2 = 4;
prior.delta_e2 = 0.01*(prior.alpha_e2-2);
prior.alpha_e3 = 4;
prior.delta_e3 = 100*(prior.alpha_e3-2);
prior.Phi2 = [0.6; 0.2];
prior.Omega2 = [4 0; 0 4];
prior.Omega2_inv = prior.Omega2\eye(2);


%%%%%======================================================================
%%%%%==================== MCMC estimating real rate trend  ================
%%%%%======================================================================

M = 100000;
[kinfQ, lamQ, gamma, Xbar, sigtau2, sige2, Phi, Omega, Sigma_tilde, tau, X, Alpha, Pit, Pic, Vt, Vc, phi1, phi2, Vptr] = MCMCsampler(Y, infl, TR, pars, prior, mats, dt, rlb, nobs, M);
save("./results/mcmc_GM2.mat")

mean(Alpha(M/2:end,:),1)















