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

UKdata = loadData(197212, 202204);
% 1y-10y yields, 198301 - 202203, 471obs
Y = UKdata{122:end,3:12};
% core CPI inflation
infl = UKdata{122:end,'Inflation'};

% USE Target Rate to help decompose inflation
% The United Kingdom adopted inflation targeting in October 1992 
% Before 199210, no target rate
% 199210-1997 a range between 1% and 4%
% 199801-200312   2.5%
% 200401-now      2%
UKdata.TR = nan(size(UKdata,1),1);
ind1 = find(UKdata.yyyymm == 199210);
ind2 = find(UKdata.yyyymm == 200401);
UKdata.TR(ind1:ind2-1) = 2.5 ;
UKdata.TR(ind2:end) = 2;
TR = UKdata{122:end,'TR'};

% index = find(~isnan(TR), 1, 'first');

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

load theta_start.mat
load theta0.mat
pars_mode = theta2pars(theta0);
pars = theta2pars(theta_start);
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.1;
pars.vc = 0.1;
pars.phi1 = 1.9;
pars.phi2 = -1.2;
pars.vptr = 0.01;

%%%%%======================================================================
%%%%%============================ 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 = [4.61 6.47 2.56]';
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 = [1.2; -0.5];
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_UK2.mat")

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




