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

% USdata = loadData(197110, 202204);
% USdata = USdata(2:end,:);
load("./data/USdata_LW2.mat")
% 3m,6m,1y-10y forward rates
Y = USdata{:,3:14};
% core PCE inflation
infl = USdata{:,'Inflation'};


% USE quarterly PTR as a proxy for long-run inflation
frbus = readtable("./data/HISTDATA.TXT");
PTR = frbus(:,{'OBS','PTR'});
PTR.date = datenum(PTR.OBS,'YYQQ');
PTR.year = year(PTR.date);
PTR.month = month(PTR.date);
PTR = table2array(PTR(:,{'year','month','PTR'}));
PTR = kron(PTR, ones(3,1));
[~,m,~] = unique(100*PTR(:,1)+PTR(:,2),'first');
m2 = m+1; m3 = m+2; 
PTR(m2,2) = PTR(m2,2)+1; PTR(m3,2) = PTR(m3,2)+2; 
PTR(m,3) = nan; PTR(m2,3) = nan; 
PTR = [PTR(:,1)*100 + PTR(:,2) PTR(:,3)];
PTR = PTR(PTR(:,1)>197109,:);
% PTR = PTR(PTR(:,1)>197110,:);
ptr = [PTR(:,2); nan(4,1)];
ptr = ptr(2:end,:);

dt = 1;
rlb = 0;
N =3;
nobs = size(Y,1);
mats = [3, 6, 12:12:120]';


load pars.mat;
load theta0 
theta_start = pars2theta(pars1);

pars_mode = theta2pars(theta0);
pars = theta2pars(theta_start);
% pars = pars2;
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.005;
pars.vc = 0.04;
pars.phi1 = 1.8;
pars.phi2 = -0.9;
pars.vptr = 0.002;

%%%%%======================================================================
%%%%%============================ 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.78 4.47 3.85]';
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 = diag([0.01; 0.01; 0.1; 0.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.01^2;
prior.kvar = 0.02^2;

% Priors for inflation decomposition
prior.alpha_e1 = 4;
prior.delta_e1 = 0.007*(prior.alpha_e1-2);
prior.alpha_e2 = 4;
prior.delta_e2 = 0.02*(prior.alpha_e2-2);
prior.alpha_e3 = 4;
prior.delta_e3 = 0.002*(prior.alpha_e3-2);
prior.Phi2 = [1.4; -0.45];
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,RMSE] = MCMCsampler(Y, infl, ptr, pars, prior, mats, dt, rlb, nobs, M);


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

save("./results/temp100000.mat")




