clear;
clc;
close all;

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

% --- Reproducibility header ---


%% |---------------- 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 
% 
y_fore    = zeros(anumber,nhor,3);
pred_Q    = zeros(anumber,nhor,length(quants),3); 
probs1    = zeros(anumber,nhor,3); 
probs2    = zeros(anumber,nhor,3);  
PITS      = zeros(anumber,nhor,3);     
QScore05  = zeros(anumber,nhor,3);
QScore10  = zeros(anumber,nhor,3);    
QScore20  = zeros(anumber,nhor,3);
QScore50  = zeros(anumber,nhor,3);
QScore80  = zeros(anumber,nhor,3);
QScore90  = zeros(anumber,nhor,3);     
QScore95  = zeros(anumber,nhor,3);
MSPE_mean = zeros(anumber,nhor,3);     
MSPE_med  = zeros(anumber,nhor,3);
CRPS      = zeros(anumber,nhor,3,2);

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

aux = T_thres:T_full;
parfor sample_p = 1:anumber 

   %rng(12345)
    sample = aux(sample_p);

    disp('  ')           
    fprintf(['<strong>YOU ARE RUNNING SAMPLE ' num2str(sample) ' OF ' num2str(T_full) '</strong> \n'] )    
    disp('  ')        
    y_true  = y(sample+1:sample+nhor,:);
    
    %% Generate Forecasts & Metrics
    predmean       = zeros(R,1);

    y_fore_   = zeros(nhor,3);
    CRPS_     = zeros(nhor,3,2);  
    pred_Q_   = zeros(nhor,length(quants),3);
    probs1_   = zeros(nhor,3);
    probs2_   = zeros(nhor,3);
    PITS_     = zeros(nhor,3);
    QScore05_ = zeros(nhor,3);
    QScore10_ = zeros(nhor,3);
    QScore20_ = zeros(nhor,3);
    QScore50_ = zeros(nhor,3);
    QScore80_ = zeros(nhor,3);
    QScore90_ = zeros(nhor,3);
    QScore95_ = zeros(nhor,3);
    MSPE_mean_ = zeros(nhor,3);
    MSPE_med_ = zeros(nhor,3);


    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);

        [draws] = f_sample_sel(X1_fore,X2_fore(:,:),out,10000,0,sampler);
        CRPS_(hor,1,1)  = f_CRPS_nrg(draws',y_true(hor));  
        CRPS_(hor,1,2)  = f_CRPS_pwm(draws',y_true(hor));

        [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';


        for r = 1:R
            predmean(r) = f_predict(X1_fore,X2_fore,beta_mixing_sim(:,:,r),beta_kernel_sim(:,:,r),tau_sim(:,:),[],'mean');
        end
        
        y_fore_(hor,1)     = mean(predmean);
       [pred_Q_(hor,:,1)]  = f_predict2(X1_fore,X2_fore,beta_mixing_sim,beta_kernel_sim,tau_sim,2,quants,'quantile');
        probs1_(hor,1)     = mean(1-f_predict(X1_fore,X2_fore,beta_mixing_sim,beta_kernel_sim,tau_sim,2,'cdf'));
        probs2_(hor,1)     = mean(1-f_predict(X1_fore,X2_fore,beta_mixing_sim,beta_kernel_sim,tau_sim,5,'cdf'));
         PITS_(hor,1)      = mean(f_predict(X1_fore,X2_fore,beta_mixing_sim,beta_kernel_sim,tau_sim,y_true(hor),'cdf'));
        QScore05_(hor,1)   = (y_true(hor) - pred_Q_(hor,1,1))*(double(y_true(hor)<=pred_Q_(hor,1,1)) - 0.05);
        QScore10_(hor,1)   = (y_true(hor) - pred_Q_(hor,2,1))*(double(y_true(hor)<=pred_Q_(hor,2,1)) - 0.10);
        QScore20_(hor,1)   = (y_true(hor) - pred_Q_(hor,3,1))*(double(y_true(hor)<=pred_Q_(hor,3,1)) - 0.20);
        QScore50_(hor,1)   = (y_true(hor) - pred_Q_(hor,4,1))*(double(y_true(hor)<=pred_Q_(hor,4,1)) - 0.5);
        QScore80_(hor,1)   = (y_true(hor) - pred_Q_(hor,5,1))*(double(y_true(hor)<=pred_Q_(hor,5,1)) - 0.80);
        QScore90_(hor,1)   = (y_true(hor) - pred_Q_(hor,6,1))*(double(y_true(hor)<=pred_Q_(hor,6,1)) - 0.90);
        QScore95_(hor,1)   = (y_true(hor) - pred_Q_(hor,7,1))*(double(y_true(hor)<=pred_Q_(hor,7,1)) - 0.95);
        MSPE_mean_(hor,1)  = (y_true(hor) - y_fore_(hor,1))^2;
        MSPE_med_(hor,1)   = (y_true(hor) - pred_Q_(hor,3,1))^2;
        predmean = f_predict(X1_fore,X2_fore,out.beta_mixing,out.beta_kernel,out.tau,[],'mean');      

     
        
        %% |---- QR Mitchell
        [yforeq,beta_mat] = VBTVPQR(y(1+hor:sample),X1_bar,X1_fore,max_iter,0.05:0.05:0.95,0);
        post_mean_quants              = yforeq;   
        pred_Q_(hor,:,2)              = post_mean_quants([1,2,4,10,16,18,19]);
        samples                       = QR1_adjusted(0.05:0.05:0.95,post_mean_quants);
        CRPS_(hor,2,1)                = f_CRPS_nrg(samples,y_true(hor));  
        CRPS_(hor,2,2)                = f_CRPS_pwm(samples,y_true(hor));

        y_fore_(hor,2)                = mean(samples);
        PITS_(hor,2)                  = ksdensity(samples,y_true(hor),'Function','cdf');
        QScore05_(hor,2)              = (y_true(hor) - post_mean_quants(1))*(double(y_true(hor)<=post_mean_quants(1)) - 0.05);
        QScore10_(hor,2)              = (y_true(hor) - post_mean_quants(2))*(double(y_true(hor)<=post_mean_quants(2)) - 0.1);
        QScore20_(hor,2)              = (y_true(hor) - post_mean_quants(3))*(double(y_true(hor)<=post_mean_quants(4)) - 0.20);
        QScore50_(hor,2)              = (y_true(hor) - post_mean_quants(10))*(double(y_true(hor)<=post_mean_quants(10)) - 0.5);
        QScore80_(hor,2)              = (y_true(hor) - post_mean_quants(16))*(double(y_true(hor)<=post_mean_quants(16)) - 0.80);
        QScore90_(hor,2)              = (y_true(hor) - post_mean_quants(18))*(double(y_true(hor)<=post_mean_quants(18)) - 0.9);
        QScore95_(hor,2)              = (y_true(hor) - post_mean_quants(19))*(double(y_true(hor)<=post_mean_quants(19)) - 0.95);
        MSPE_mean_(hor,2)             = (y_true(hor) - post_mean_quants(10))^2;
        MSPE_med_(hor,2)              = (y_true(hor) - post_mean_quants(10))^2;

        
     
        %% |---- 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);
         predsAR           = 1*beta_mat(end,:)' + sqrt(sigma_t(end)).*randn(1,10000);
         CRPS_(hor,3,1)    = f_CRPS_nrg(predsAR,y_true(hor));
         CRPS_(hor,3,2)    = f_CRPS_pwm(predsAR,y_true(hor));

         pred_Q_(hor,:,3)  = norminv(quants,1*beta_mat(end,:)',sqrt(sigma_t(end)));
         probs1_(hor,3)    = 1-normcdf(2,1*beta_mat(end,:)',sqrt(sigma_t(end)));
         probs2_(hor,3)    = 1-normcdf(5,1*beta_mat(end,:)',sqrt(sigma_t(end)));
         PITS_(hor,3)      = normcdf(y_true(hor),1*beta_mat(end,:)',sqrt(sigma_t(end)));%ksdensity(squeeze(pred_Q(sample-T_thres+1,hor,:,3)),y_true(hor),'Function','cdf');
         QScore05_(hor,3)  = (y_true(hor) - pred_Q_(hor,1,3))*(double(y_true(hor)<=pred_Q_(hor,1,3)) - 0.05);
         QScore10_(hor,3)  = (y_true(hor) - pred_Q_(hor,2,3))*(double(y_true(hor)<=pred_Q_(hor,2,3)) - 0.1);
         QScore20_(hor,3)  = (y_true(hor) - pred_Q_(hor,3,3))*(double(y_true(hor)<=pred_Q_(hor,3,3)) - 0.20);
         QScore50_(hor,3)  = (y_true(hor) - pred_Q_(hor,4,3))*(double(y_true(hor)<=pred_Q_(hor,4,3)) - 0.5);
         QScore80_(hor,3)  = (y_true(hor) - pred_Q_(hor,5,3))*(double(y_true(hor)<=pred_Q_(hor,5,3)) - 0.80);
         QScore90_(hor,3)  = (y_true(hor) - pred_Q_(hor,6,3))*(double(y_true(hor)<=pred_Q_(hor,6,3)) - 0.9);
         QScore95_(hor,3)  = (y_true(hor) - pred_Q_(hor,7,3))*(double(y_true(hor)<=pred_Q_(hor,7,3)) - 0.95);
         MSPE_mean_(hor,3) = (y_true(hor) - y_fore_(hor,3))^2;
         MSPE_med_(hor,3)  = (y_true(hor) - pred_Q_(hor,3,3))^2;

         
    

    end

     at = sample-T_thres+1;
     CRPS(sample_p,:,:,:)    = CRPS_;
     pred_Q(sample_p,:,:,:)  = pred_Q_;
     probs1(sample_p,:,:)    = probs1_;
     probs2(sample_p,:,:)    = probs2_;
     PITS(sample_p,:,:)      = PITS_;
     QScore05(sample_p,:,:)  = QScore05_;
     QScore10(sample_p,:,:)  = QScore10_;
     QScore20(sample_p,:,:)  = QScore20_;
     QScore50(sample_p,:,:)  = QScore50_;
     QScore80(sample_p,:,:)  = QScore80_;
     QScore90(sample_p,:,:)  = QScore90_;
     QScore95(sample_p,:,:)  = QScore95_;
     MSPE_mean(sample_p,:,:) = MSPE_mean_;
     MSPE_med(sample_p,:,:)  = MSPE_med_;
end

save(infl_id + "_H_"+ string(H_max) + "_priorM_"+ string(prior_mixing) + "_priorK_"+ string(prior_kernel) + "_" +string(init_method)+"_.mat");

%% Compute RS test statistics

rvec = [0:0.001:1]; 
for h = 1:nhor
    [outtest] = rs_test_adj(PITS(:,h,1), rvec, h>1, 1);
    teststats(h,1,1,1) = outtest.teststat.ks;
    teststats(h,1,2,1) = outtest.teststat.cvm;
    rej_tab(h,:,1,1,1) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,1,2,1) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,2), rvec, h>1, 1);
    teststats(h,2,1,1) = outtest.teststat.ks;
    teststats(h,2,2,1) = outtest.teststat.cvm;
    rej_tab(h,:,2,1,1) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,2,2,1) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,3), rvec, h>1, 1);
    teststats(h,3,1,1) = outtest.teststat.ks;
    teststats(h,3,2,1) = outtest.teststat.cvm;
    rej_tab(h,:,3,1,1) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,3,2,1) = outtest.teststat.cvm > outtest.crit_value.cvm;


    critvals(h,:,1,1) = outtest.crit_value.ks;
    critvals(h,:,2,1) = outtest.crit_value.cvm;


    [outtest] = rs_test_adj(PITS(:,h,1), rvec, h>1, 6);
    teststats(h,1,1,2) = outtest.teststat.ks;
    teststats(h,1,2,2) = outtest.teststat.cvm;
    rej_tab(h,:,1,1,2) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,1,2,2) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,2), rvec, h>1, 6);
    teststats(h,2,1,2) = outtest.teststat.ks;
    teststats(h,2,2,2) = outtest.teststat.cvm;
    rej_tab(h,:,2,1,2) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,2,2,2) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,3), rvec, h>1, 6);
    teststats(h,3,1,2) = outtest.teststat.ks;
    teststats(h,3,2,2) = outtest.teststat.cvm;
    rej_tab(h,:,3,1,2) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,3,2,2) = outtest.teststat.cvm > outtest.crit_value.cvm;


    critvals(h,:,1,2) = outtest.crit_value.ks;
    critvals(h,:,2,2) = outtest.crit_value.cvm;


    [outtest] = rs_test_adj(PITS(:,h,1), rvec, h>1, 7);
    teststats(h,1,1,3) = outtest.teststat.ks;
    teststats(h,1,2,3) = outtest.teststat.cvm;
    rej_tab(h,:,1,1,3) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,1,2,3) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,2), rvec, h>1, 7);
    teststats(h,2,1,3) = outtest.teststat.ks;
    teststats(h,2,2,3) = outtest.teststat.cvm;
    rej_tab(h,:,2,1,3) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,2,2,3) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,3), rvec, h>1, 7);
    teststats(h,3,1,3) = outtest.teststat.ks;
    teststats(h,3,2,3) = outtest.teststat.cvm;
    rej_tab(h,:,3,1,3) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,3,2,3) = outtest.teststat.cvm > outtest.crit_value.cvm;


    critvals(h,:,1,3) = outtest.crit_value.ks;
    critvals(h,:,2,3) = outtest.crit_value.cvm;



    [outtest] = rs_test_adj(PITS(:,h,1), rvec, h>1, 2);
    teststats(h,1,1,4) = outtest.teststat.ks;
    teststats(h,1,2,4) = outtest.teststat.cvm;
    rej_tab(h,:,1,1,4) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,1,2,4) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,2), rvec, h>1, 2);
    teststats(h,2,1,4) = outtest.teststat.ks;
    teststats(h,2,2,4) = outtest.teststat.cvm;
    rej_tab(h,:,2,1,4) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,2,2,4) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,3), rvec, h>1, 2);
    teststats(h,3,1,4) = outtest.teststat.ks;
    teststats(h,3,2,4) = outtest.teststat.cvm;
    rej_tab(h,:,3,1,4) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,3,2,4) = outtest.teststat.cvm > outtest.crit_value.cvm;


    critvals(h,:,1,4) = outtest.crit_value.ks;
    critvals(h,:,2,4) = outtest.crit_value.cvm;



    [outtest] = rs_test_adj(PITS(:,h,1), rvec, h>1, 5);
    teststats(h,1,1,5) = outtest.teststat.ks;
    teststats(h,1,2,5) = outtest.teststat.cvm;
    rej_tab(h,:,1,1,5) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,1,2,5) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,2), rvec, h>1, 5);
    teststats(h,2,1,5) = outtest.teststat.ks;
    teststats(h,2,2,5) = outtest.teststat.cvm;
    rej_tab(h,:,2,1,5) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,2,2,5) = outtest.teststat.cvm > outtest.crit_value.cvm;

    [outtest] = rs_test_adj(PITS(:,h,3), rvec, h>1, 5);
    teststats(h,3,1,5) = outtest.teststat.ks;
    teststats(h,3,2,5) = outtest.teststat.cvm;
    rej_tab(h,:,3,1,5) = outtest.teststat.ks > outtest.crit_value.ks;
    rej_tab(h,:,3,2,5) = outtest.teststat.cvm > outtest.crit_value.cvm;


    critvals(h,:,1,5) = outtest.crit_value.ks;
    critvals(h,:,2,5) = outtest.crit_value.cvm;


end

models = ["DR", "QR", "UC-SV"];
j = 1;
for i = [1,2,3,4,8,12]
    tt = squeeze(teststats(i,[1,2,3],:,[1,2,4,5]));
    gg = squeeze(critvals(i,:,:,[1,2,4,5]));
    [formatted(:,:,j)] = significanceStarsMatrix(tt(:,:), repmat(permute(flipud(gg(:,:)),[3,2,1]),4,1),models);
    j = j+1;
end

% Latex table
[formatted(:,1:end-1,1) formatted(:,2:end,2)];
[formatted(:,1:end-1,3) formatted(:,2:end,4)];
[formatted(:,1:end-1,5) formatted(:,2:end,6)];


%% Compute Quantile Scores & CRPS

j = 1;
for i = [1,2,3,4,8,12]
    QS(:,:,j) = [mean(squeeze(QScore05(:,i,:)))' mean(squeeze(QScore10(:,i,:)))' mean(squeeze(QScore20(:,i,:)))' mean(squeeze(QScore50(:,i,:)))' mean(squeeze(QScore80(:,i,:)))' mean(squeeze(QScore90(:,i,:)))' mean(squeeze(QScore95(:,i,:)))' mean(squeeze(CRPS(:,i,:,1)))'];
    j = j+1;
end

QSstr = [repmat(cellstr(models)',1,1,6) cellfun(@(x) sprintf('%.2f', x), num2cell(round(QS([1,2,3],:,:),2)), 'UniformOutput', false)];

% Size info
[nRows, nCols, nLayer] = size(QSstr);

% New array: double the number of columns
B = cell(nRows, 2 * nCols, nLayer);

% Fill original columns into odd positions
B(:, 1:2:end, :) = QSstr;

% Fill symbol into even positions
B(:, 2:2:end,:) = repmat({'&'}, nRows, nCols, nLayer);

QSstr = B;
QSstr(:,end,:) = repmat({'"\\"'},nRows,nLayer);

% Latex table
[QSstr(:,1:end-1,1) QSstr(:,2:end,2) ];
[QSstr(:,1:end-1,3) QSstr(:,2:end,4) ];
[QSstr(:,1:end-1,5) QSstr(:,2:end,6) ];


%% Compute Quantile Scores T-tests

[H,P,CI,STATS] = ttest2(QScore05(:,:,1),QScore05(:,:,3));
t_stat05(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore10(:,:,1),QScore10(:,:,3));
t_stat10(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore20(:,:,1),QScore20(:,:,3));
t_stat20(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore50(:,:,1),QScore50(:,:,3));
t_stat50(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore80(:,:,1),QScore80(:,:,3));
t_stat80(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore90(:,:,1),QScore90(:,:,3));
t_stat90(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore95(:,:,1),QScore95(:,:,3));
t_stat95(1,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(CRPS(:,:,1,1),CRPS(:,:,3,1));
t_CRPS(1,:) = squeeze(STATS.tstat);

[H,P,CI,STATS] = ttest2(QScore05(:,:,2),QScore05(:,:,3));
t_stat05(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore10(:,:,2),QScore10(:,:,3));
t_stat10(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore20(:,:,2),QScore20(:,:,3));
t_stat20(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore50(:,:,2),QScore50(:,:,3));
t_stat50(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore80(:,:,2),QScore80(:,:,3));
t_stat80(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore90(:,:,2),QScore90(:,:,3));
t_stat90(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(QScore95(:,:,2),QScore95(:,:,3));
t_stat95(2,:) = squeeze(STATS.tstat);
[H,P,CI,STATS] = ttest2(CRPS(:,:,2,1),CRPS(:,:,3,1));
t_CRPS(2,:) = squeeze(STATS.tstat);

tstats = [permute(t_stat05,[1,3,2]) permute(t_stat10,[1,3,2]) permute(t_stat20,[1,3,2]) permute(t_stat50,[1,3,2]) permute(t_stat80,[1,3,2]) permute(t_stat90,[1,3,2]) permute(t_stat95,[1,3,2]) permute(t_CRPS,[1,3,2])];

for h=1:nhor
    for i = 1:size(tstats,1)
        for j = 1:size(tstats,2)
            formattedStr = sprintf('%.2f%s', tstats(i,j,h));
            if abs(tstats(i,j,h)) >= 1.96
                formattedStr = ['\textbf{' formattedStr '}'];
            end
            formattedTstat{i,j,h} = formattedStr;
        end
    end
end


% Size info
[nRows, nCols, nLayer] = size(formattedTstat);

% New array: double the number of columns
B = cell(nRows, 2 * nCols, nLayer);

% Fill original columns into odd positions
B(:, 1:2:end, :) = formattedTstat;

% Fill symbol into even positions
B(:, 2:2:end,:) = repmat({'&'}, nRows, nCols, nLayer);

formattedTstat = B;
formattedTstat(:,end,:) = repmat({'"\\"'},nRows,nLayer);
formattedTstat = [repmat(cellstr(models(1:2))',1,1,12) repmat({'&'},2,1,12) formattedTstat ];

% Latex table
[formattedTstat(:,1:end-1,1) formattedTstat(:,2:end,2) ];
[formattedTstat(:,1:end-1,3) formattedTstat(:,2:end,4) ];
[formattedTstat(:,1:end-1,8) formattedTstat(:,2:end,12)];


%% Compute RMSE and DM test statistics

RMSE = array2table(sqrt(squeeze(mean(MSPE_mean(:,[1:4 8 12],[1,2,3])))),'RowNames',["h=1"; "h=2"; "h=3"; "h=4"; "h=8"; "h=12"],'VariableNames',["DR"; "QR_M"; "UC-SV"]);
absRMSE = sqrt(squeeze(mean(MSPE_mean(:,:,[1,2,3]))));
relRMSE = sqrt(squeeze(mean(MSPE_mean(:,:,[1,2]))))./sqrt(squeeze(mean(MSPE_mean(:,:,3))))';

DMtests = zeros(nhor,2);
aux = MSPE_mean(:,:,[1,2,3]);
for m = 1:2
    for h = 1:nhor
        [~,DMtests(h,m)] = dmtest(squeeze(aux(:,h,m)),squeeze(MSPE_mean(:,h,3)));
    end
end

stars = {'*', '**', '***'};
for i = 1:size(relRMSE,1)
    for j = 1:size(relRMSE,2)

        if DMtests(i,j) < 0.01
            star = '$^{***}$';
        elseif DMtests(i,j) < 0.05
            star = '$^{**\phantom{*}}$';
        elseif DMtests(i,j) < 0.1
            star = '$^{*\phantom{**}}$';
        else
            star = '$^{\phantom{***}}$';
        end
        % Bold if significant at 5%
        formattedStr = sprintf('%.2f%s', relRMSE(i,j), star);
        formattedRMSE{i,j} = formattedStr;
    end
end

hors = [ "h=1 $^{\phantom{***}}$"; "h=2 $^{\phantom{***}}$" ; "h=3 $^{\phantom{***}}$" ;...
         "h=4 $^{\phantom{***}}$" ; "h=5 $^{\phantom{***}}$" ; "h=6 $^{\phantom{***}}$" ;...
         "h=7 $^{\phantom{***}}$"; "h=8 $^{\phantom{***}}$"; "h=9 $^{\phantom{***}}$" ;...
         "h=10 $^{\phantom{***}}$" ; "h=11 $^{\phantom{***}}$" ; "h=12 $^{\phantom{***}}$"];


models = ["DR"    "QR"];

formattedRMSE = [["" hors' ]; [models' formattedRMSE']];

% Size info
[nRows, nCols, nLayer] = size(cellstr(formattedRMSE));

% New array: double the number of columns
B = cell(nRows, 2 * nCols, nLayer);

% Fill original columns into odd positions
B(:, 1:2:end, :) = cellstr(formattedRMSE);

% Fill symbol into even positions
B(:, 2:2:end,:) = repmat({'&'}, nRows, nCols, nLayer);

% Latex table
formattedRMSE = B;
formattedRMSE(:,end,:) = repmat({'"\\"'},nRows,nLayer);

% Latex table
formattedRMSE;