%% clearing the session
clc
clear 
close all

addpath('./mtimesx/')
addpath('./functions/')

%% load data
load dataHeadline


%% model settings
p=6; q=6; qp=[q,p];
model='MAIARSV';
setenv('model',model);
fblks=[20];

%% prepare data
mfirstlag=ones(n,1);
Y_nstd=(rawdata)';
Y=standardize(demean(Y_nstd,2),2);

%% prior building
t_end_train=length(dateslab);
prior=build_prior_std_MAIAR_SV(Y',fblks,qp,t_end_train,mfirstlag);
prior.mB0=make_B0_restricted(prior.mb0,fblks);
prior.vb0=n*prior.vb0;


%% setting GS properties
store_bvec=[1,1,1,1,1,0];
c_metro=2.5;   
l_c=2e4;
logfilepath='./log/MAIARSV_infla.txt';


%% running GS
Y_est=Y;
tic
    [gammas,as,b0_tildes,acc,gs,sigmas,Qss]=...
    MAIARSV(Y_est,qp,fblks,prior,c_metro,l_c,logfilepath,store_bvec);
toc

disp('Acceptance rates for B0:')
disp(acc')

%% posterior summaries


quants=[0.05,0.16,0.5,0.84,0.95];
nrp=size(as,1); n_r=size(b0_tildes,1); m=size(gs,1); nq=size(gammas,1);

msigmas=median(sigmas,3);
mgammas=median(gammas,2);
mas=median(as,2);
mb0=median(b0_tildes,2);
mgs=median(gs,2);

qsigmas=quantile(sigmas,quants,3);
qgammas=quantile(gammas,quants,2);
qas=quantile(as,quants,2);
qb0=quantile(b0_tildes,quants,2);
qgs=quantile(gs,quants,2);

nbins=60;
c_hist_as=cell(2,nrp);
c_hist_gammas=cell(2,nq);
c_hist_b0s=cell(2,n_r);
c_hist_gs=cell(2,m);

for j=1:nq
    [c_hist_gammas{2,j},c_hist_gammas{1,j}] = histnorm(gammas(j,:), nbins);
end

for j=1:nrp
    [c_hist_as{2,j},c_hist_as{1,j}] = histnorm(as(j,:), nbins);
end
for j=1:n_r
    [c_hist_b0s{2,j},c_hist_b0s{1,j}] = histnorm(b0_tildes(j,:), nbins);
end

for j=1:m
    [c_hist_gs{2,j},c_hist_gs{1,j}] = histnorm(gs(j,:), nbins);
end


tab_B0_summary=[table( make_B0_restricted(prior.mb0,fblks)','RowNames',short_desc,'VariableNames',{'factors_prior'} ),...
table( make_B0_restricted(mb0,fblks)','RowNames',short_desc,'VariableNames',{'factors_posterior'} )];

disp(tab_B0_summary)

%% build factors

factors=makeFactors(b0_tildes,fblks,Y);
mfactor=median(factors,3);
qfactors=squeeze(quantile(factors,quants,3));

