clear;
clc;
close all;

addpath('data')
addpath('functions')
addpath('figures')
addpath('results')
addpath('RossiSekhposyan')
addpath(genpath('VulnerableGrowth'))

%% |---------------- Load Data

% loads data for convenience (here PCE inflation excluding energy and food)
% to run a different aggregate, load/compase the corresponding data
load('data_for_server_20250429_PCEx.mat');


% ======> to generate data, uncomment below block
% 
% settings
% infl_id     = 'PCEPILFE'; 
% lags        = 4;
% start_date  = '01-Apr-1973';
% end_date    = '01-Oct-2024';
% standardize = 2; 
% 
% % selected indicators
% select = {'OECD_6NMEIndustrialProduction',...
%           'Energy',...
%           'Food',...
%           'Metals',...
%           'EBP',...
%           'DFF',...
%           'SPF_linked_1YR',...
%           'PTR'};
% 
% % transformations
% trans = [4,4,4,4,1,1,1,1,1];
% 
% % construct data
% [y,x,n,T,dates_] = load_data(infl_id,lags,start_date,end_date,standardize,select,trans);


%% |---------------- Estimation Settings


nhor         = 12;
H_max        = 5;            % Number of mixing components 
init_method  = 'cluster';   % cluster initialization method random|cluster|standard
iter         = 20;           % Print every "iter^th" iteration
nthin        = 100;
prior_mixing = 2;           % prior for mixing coefficients: 0: Normal Prior | 1: Tipping's ML prior | 2: Horseshoe Prior
prior_kernel = 0;           % prior for mixing coefficients: 0: Normal Prior | 1: Tipping's ML prior | 2: Horseshoe Prior
basis_exp    = 2;           % 0: no basis function expansion | 1: basis function expansion | 2: panel
max_iter     = 10000;
tol          = 1e-6;        % 1e-10;
R            = 1000;
quants       = [0.05,0.1,0.2,0.5,0.8,0.9,0.95]; 

% Start with  first_per% of observations for first period forecast
first_per = 0.5;
T_full    = T-nhor;
T_thres   = round(first_per*T_full); 
anumber   = T_full-T_thres+1;
varnum    = size(x,2);

X1 = x;
X2 = x;

% prior settings
prior.a_tau = 1;
prior.b_tau = 0.1;

prior.Ab_mix  = 1;
prior.tau_mix = 0.1;     
prior.nu_mix  = 0.1;  

sampler       = 2; 

%% Forecasting 

HHs = zeros(anumber,nhor);
quantQR = [0.05,0.1,0.2,0.5,0.8,0.9,0.95];

aux = T_thres:T_full;


sample = T;
for hor = 1:nhor
   
    %% |---- Density Regression
    
     X1_bar         = X1(1:sample-hor,:);
     X2_bar         = X2(1:sample-hor,:);
     X1_fore        = X1(sample,:);
     X2_fore        = X2(sample,:);

    [out] = VB_BDR_pipe(y(1+hor:sample),X1_bar,X2_bar(:,:),H_max,init_method,prior_mixing,prior_kernel,10000,tol,prior);

   
    [beta_mixing_sim,beta_kernel_sim,tau_sim] = f_simul_VB(out,R);
    beta_kernel_sim = permute(beta_kernel_sim,[2,3,1]);
    beta_mixing_sim = permute(beta_mixing_sim,[2,3,1]);
    tau_sim = tau_sim';

   [pred_Q_(hor).quantile(:,:,1)] = f_predict2(X1,X2,beta_mixing_sim,beta_kernel_sim,tau_sim,2,quants,'quantile');
  
    %% |---- QR Mitchell
    [yforeq,beta_mat] = VBTVPQR(y(1+hor:sample),X1_bar,X1_bar,max_iter,quants,0);
    pred_Q_(hor).quantile(:,:,2)              = X1*squeeze(beta_mat(end,:,:));   
      
    %% |---- UCSV
    [y_fore_(hor,3),beta_mat,~,sigma_t,~] =  VB_DVS(y(1+hor:sample),ones(size(X1_bar,1),1),1,1,max_iter,2,.8,.1,.1);
     
    for ii = 1:T
        if ii >hor %length(beta_mat)
            beta = beta_mat(ii-hor,:)';
            sigma = sqrt(sigma_t(ii-hor));
            
            pred_Q_(hor).quantile(ii,:,3)  = norminv(quants,1*beta',sqrt(sigma));
        end
    end
 
end

%save("quantile_chart.mat");

chart_data = pred_Q_(4).quantile(:,[2,4,6],:);
groupColors = get(groot,'defaultAxesColorOrder');


figure()
hold on
for g = 1:3
    plot(dates_,chart_data(:,g,1),'-','Color',groupColors(g,:))
    % Shaded band between the two
    fill([dates_' fliplr(dates_')], [chart_data(:,g,1)' fliplr(chart_data(:,g,2)')], ...
        groupColors(g,:), 'FaceAlpha', 0.15, 'EdgeColor', 'none');
    ax = gca;
    ax.ColorOrder = get(groot, 'defaultAxesColorOrder');
    ax.ColorOrderIndex = 1;
    plot(dates_,chart_data(:,g,2),'--','Color',groupColors(g,:))
    %plot(dates_,chart_data(:,g,3),'--','Color','black')
end
hold off
legend(["10th percentile - DR","","10th percentile - QR","50th percentile - DR","","50th percentile - QR","90th percentile - DR","","90th percentile - QR"],"NumColumns",3,"Location","southoutside","FontSize",8)
% Set renderer to 'painters' for vector output (PDF/EPS)
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 20, 6]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');

% Save as vector PDF (preferred for LaTeX and journals)
print(gcf, 'OS_Figure3_quantile_chart.pdf', '-dpdf', '-painters');
set(gcf,'PaperPositionMode','auto');
print(gcf, "OS_Figure3_quantile_chart.eps", '-depsc', '-painters','-loose');
