clear;
clc;
close all;

addpath('data')
addpath('functions')
addpath('figures')
addpath('results')


%% |---------------- 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_risk_PCEx.mat');
meas = "PCEx";

% settings
% infl_id     = 'CPILFESL'; 
% lags        = 4;
% start_date  = '01-Apr-1973';
% end_date    = '01-Oct-2024';
% standardize = 2; 
% 
% % selected indicators
% select = {'OECD_6NMEIndustrialProduction',...
%           'Energy',...
%           'Food',...
%           'Metals',...
%           'EBP',...
%           'MARn',...
%           '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);
% 
% save('data_for_server_20250429_risk_CPIx.mat');

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

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

hor          = 1;            % forecast horizon
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
max_iter     = 10000;
tol          = 1e-6;         % estimation tolerance 
R            = 1000;         % number of posterior draws
M            = 1000;         % number of draws for decompositions

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

X1 = x;
X2 = x;

%% Forecasting 
nhor  = length(hor);  
tdates = [];
fdates = [];
yts    = zeros(anumber,nhor);

yyear = 0;   

for sample = T_thres:T_full
    disp('  ')           
    fprintf(['<strong>YOU ARE RUNNING SAMPLE ' num2str(sample) ' OF ' num2str(T_full) '</strong> \n'] )    
    disp('  ')
    if sample+hor<=T_full
        y_true  = y(sample+hor,:);
        tdates = [tdates; string(dates_(sample+hor,:))'];
    else
        y_true = NaN;
        if month(datetime(tdates(end-4))+calmonths(3))==10
            yyear = yyear+1;
            tdates = [tdates; string(datetime(tdates(end-4))+calmonths(3)+calyears(yyear))];
        else
            tdates = [tdates; string(datetime(tdates(end-4))+calmonths(3)+calyears(yyear))];
        end
    end
    yts(sample-T_thres+1,:) = y_true;
    fdates = [fdates; string(dates_(sample))];
    
    %% Generate Forecasts & Metrics
    predmean       = zeros(R,1);

    %% |---- 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,max_iter,tol,prior);
    res(sample-T_thres+1,nhor).DR = out;
    
    %% |---- 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(sample-T_thres+1,hor,:,2)              = post_mean_quants([1,2,4,10,16,18,19]);
    samples                                       = QR1_adjusted(0.05:0.05:0.95,post_mean_quants);

    outQRM.beta                    = beta_mat;
    outQRM.samples                 = samples;
    outQRM.post_mean_quants        = post_mean_quants;
    res(sample-T_thres+1,nhor).QRM = outQRM;

    %% |---- UCSV
    [y_fore(sample-T_thres+1,hor,4),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);
    
    outUC.beta                    = beta_mat;
    outUC.sigmat                  = sigma_t;
    outUC.preds                   = predsAR;
    res(sample-T_thres+1,nhor).UC = outUC;  
   
end

save(["results_realtime_decompositions_h1_"+string(meas)])

%% historical densities example
 
% target quarters
dens_dates = {'01-Jan-2021',...
              '01-Oct-2022'};


bin_edges = [-Inf -20:0.5:20 Inf];
B=length(bin_edges)-1;
bins_center = 0.5 * (bin_edges(1:end-1) + bin_edges(2:end));

groups = [1 1 1 1 2 3 4 4 4 5 6 7 7]; % specified without the constant
clamped = 0;
clampid = 80;
method  = "random";
parfor t = 1:length(dens_dates)
    sample = find(dates_==dens_dates(t));
    fprintf(["vintage: " + string(t) + " of " + string(length(dens_dates)) + "\n"])

    f = @f_DR;
    if clamped ==1
        [shaps(t).DR] = shap_dens_comp3(x(sample,:),x(clampid:sample-1,:),f,res(sample-T_thres+1,1).DR,1,bin_edges,method,1,0,groups,[]);
    else
        [shaps(t).DR] = shap_dens_comp3(x(sample,:),x(1:sample-1,:),f,res(sample-T_thres+1,1).DR,1,bin_edges,method,1,0,groups,[]);
    end

    f = @f_QR3;
    if clamped ==1
        [shaps(t).QRM] = shap_dens_comp3(x(sample,:),x(clampid:sample-1,:),f,squeeze(res(sample-T_thres+1,1).QRM.beta(end,:,:)),1,bin_edges,method,1,0,groups,[]);
    else
        [shaps(t).QRM] = shap_dens_comp3(x(sample,:),x(1:sample-1,:),f,squeeze(res(sample-T_thres+1,1).QRM.beta(end,:,:)),1,bin_edges,method,1,0,groups,[]);
    end

end

%save('dens_decomp_results')
save("density_decomposition_h1_"+string(meas)+"_" +method +"_clamped_"+ string(clamped))

labels = {"Dynamic Persistence";...
          "Domestic Business Cycle";...
          "Global Business Cycle";...
          "Commodity Prices";...
          "Financial Conditions";...
          "Monetary Policy";...
          "Inflation Expectations";...
          "Forecast Density";...
          "Background Data Density";...
          "Δ Density"};



colors = [ ...    
    0.1216, 0.4667, 0.7059;  % Blue
    1.0000, 0.4980, 0.0549;  % Orange
    0.1725, 0.6275, 0.1725;  % Green
    0.8392, 0.1529, 0.1569;  % Red
    0.5804, 0.4039, 0.7412;  % Purple
    0.5490, 0.3373, 0.2941;  % Brown
    0.8902, 0.4667, 0.7608;  % Pink
    0.4980, 0.4980, 0.4980;  % Gray
    0.60 0.70 0.50;  % Olive
    0.00 0.47 0.52]; % Cyan/Teal


%% Plotting density decompositions

for vint = 1:size(dens_dates,2)
    
    %%%%%%%%%%% DR chart

    figure()
    %tcl = tiledlayout(1,2);
    %nexttile(tcl)
    if method == "knn"
        dp = shaps(vint).DR.p_full - shaps(vint).DR.phi0cond;
    else
        dp = shaps(vint).DR.p_full - shaps(vint).DR.phi0;
    end
    bardata = shaps(vint).DR.phi(2:end,:)';
    b=bar(bins_center, bardata', 'stacked');
    for i = 1:7
        b(i).FaceColor = colors(i,:);
        b(i).EdgeColor = 'none';
    end
    %title('Density regression');
    ylabel('Δ Probability');
    hold on; 
    plot(bins_center, shaps(vint).DR.p_full, 'Color',colors(9,:), 'LineWidth', 2);
    if method == "knn"
        plot(bins_center, shaps(vint).DR.phi0cond, 'Color',colors(10,:), 'LineWidth', 2); 
    else
        plot(bins_center, shaps(vint).DR.phi0, 'Color',colors(10,:), 'LineWidth', 2); 
    end
    plot(bins_center, dp, '--','Color','black', 'LineWidth', 2); 
    xlim([-5,10])
    grid on;
    
    % Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
    set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
    
    % Set renderer to 'painters' for vector output (PDF/EPS)
    set(gcf, 'Renderer', 'painters');
    
    % Save as vector PDF (preferred for LaTeX and journals)
    drawnow;
    print(gcf, "shap_dens_decomposition_plot_" + string(year(datetime(dens_dates(vint)))) + "_" + string(month(datetime(dens_dates(vint)))) + "_" + meas + "_" + method +"_DR.pdf", '-dpdf', '-painters');
    
    set(gcf,'PaperPositionMode','auto');
    print(gcf, "shap_dens_decomposition_plot_" + string(year(datetime(dens_dates(vint)))) + "_" + string(month(datetime(dens_dates(vint)))) + "_" + meas + "_" + method +"_DR.eps", '-depsc', '-painters','-loose');
    
    %%%%%%%%%%% QR chart

    bardata = [shaps(vint).QRM.phi(2:end,:)'];
    if method == "knn"
        dp = shaps(vint).QRM.p_full - shaps(vint).QRM.phi0cond;
    else
        dp = shaps(vint).QRM.p_full - shaps(vint).QRM.phi0;
    end
    
    figure()
    %nexttile(tcl)
    b=bar(bins_center, bardata', 'stacked');
    for i = 1:7
        b(i).FaceColor = colors(i,:);
        b(i).EdgeColor = 'none';
    end
    %title('Quantile regression');
    ylabel('Δ Probability');
    hold on; 
    h1=plot(bins_center, shaps(vint).QRM.p_full, 'Color',colors(9,:), 'LineWidth', 2); 
    if method == "knn"
        h2=plot(bins_center, shaps(vint).QRM.phi0cond, 'Color',colors(10,:), 'LineWidth', 2);  
    else
        h2=plot(bins_center, shaps(vint).QRM.phi0, 'Color',colors(10,:), 'LineWidth', 2);  
    end
    h3=plot(bins_center, dp, '--','Color','black', 'LineWidth', 2); 
    xlim([-5,10])
    %legend(labels)
    grid on;
    % Construct a Legend with the data from the sub-plots
    %hL = legend(labels,'NumColumns',4); 
    % Move the legend to the right side of the figure
    %hL.Layout.Tile = 'south';
    
    % Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
    set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
    
    % Set renderer to 'painters' for vector output (PDF/EPS)
    set(gcf, 'Renderer', 'painters');
    
    % Save as vector PDF (preferred for LaTeX and journals)
    drawnow;
    print(gcf, "shap_dens_decomposition_plot_" + string(year(datetime(dens_dates(vint)))) + "_" + string(month(datetime(dens_dates(vint)))) + "_" + meas + "_" + method +"_QR.pdf", '-dpdf', '-painters');

    set(gcf,'PaperPositionMode','auto');
    print(gcf, "shap_dens_decomposition_plot_" + string(year(datetime(dens_dates(vint)))) + "_" + string(month(datetime(dens_dates(vint)))) + "_" + meas + "_" + method +"_QR.eps", '-depsc', '-painters','-loose');
end



%% Loss decomposition

loss_params.lower = 2;
loss_params.upper = 2;
loss_params.alpha = 2;
loss_params.beta  = 2;
loss_params.w     = 0.5;


T_Thres = 144;
numSamples = T_full - T_thres + 1;
RM = repmat(struct('DR', [], 'QR', []), numSamples, 1);
clamped = 1;
clampid = 80;
method = "random";
parfor k = 1:numSamples
    sample = T_thres + k - 1;
    fprintf("vintage: %d of %d\n", sample, T_full);

    f1 = @f_loss_DR;
    f2 = @f_loss_QR;

    x_sample = x(sample,:);
    if clamped == 1
        x_hist = x(clampid:sample-1,:);
    else
        x_hist = x(1:sample-1,:);
    end

    res_DR = res(k, nhor).DR;
    res_QR_beta = squeeze(res(k, nhor).QRM.beta(end,:,:));

    RM_k = struct();
    RM_k.DR = shap_loss_comp3(x_sample, x_hist, f1, res_DR, loss_params, 1, method, 1, 0, groups, 20000, 1);
    RM_k.QR = shap_loss_comp3(x_sample, x_hist, f2, res_QR_beta, loss_params, 1, method, 1, 0, groups, 10000, 1);
    RM(k) = RM_k;
end

% save results
save("loss_decomposition_h1_clamped_" + string(clamped) + "_"+ string(meas)+ "_"+ string(loss_params.alpha) + string(loss_params.beta) + "_" + method + "_" );


%% Loss decomposition charts

for k = 1:numSamples
    BRs(:,k,1) = RM(k).DR.phi(:,3);
    BRs(:,k,2) = RM(k).QR.phi(:,3);
    if method == "knn"
        phi0s(:,k,1)  = RM(k).DR.phi0cond(:,3);
        phi0s(:,k,2)  = RM(k).QR.phi0cond(:,3);
    else
        phi0s(:,k,1)  = RM(k).DR.phi0(:,3);
        phi0s(:,k,2)  = RM(k).QR.phi0(:,3);
    end
    pfulls(:,k,1)  = RM(k).DR.p_full(:,3);
    pfulls(:,k,2)  = RM(k).QR.p_full(:,3);
end

figure()
tcl = tiledlayout(2,1);
nexttile(tcl)
bardata = [BRs(2:end,:,1)];
b=bar(dates_(T_thres:T_full), bardata', 'stacked');
for i = 1:7
b(i).FaceColor = colors(i,:);
end
title('Density regression');
ylabel('Δ Balance of Risk');
hold on;
plot(dates_(T_thres:T_full), pfulls(:,:,1)-phi0s(:,:,1), '--','Color','red', 'LineWidth', 2);
grid on;
nexttile(tcl)
bardata = [BRs(2:end,:,2)];
b=bar(dates_(T_thres:T_full), bardata', 'stacked');
for i = 1:7
b(i).FaceColor = colors(i,:);
end
title('Quantile regression');
ylabel('Δ Balance of Risk');
hold on;
plot(dates_(T_thres:T_full), pfulls(:,:,2)-phi0s(:,:,2), '--','Color','red', 'LineWidth', 2);
grid on;
hL = legend(labels,'NumColumns',4);
% Move the legend to the right side of the figure
hL.Layout.Tile = 'south';
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "shap_decomposition_loss_plot" +string(meas)+ ".pdf", '-dpdf', '-painters');
set(gcf,'PaperPositionMode','auto');
print(gcf, "shap_decomposition_loss_plot" +string(meas)+ ".eps", '-depsc', '-painters','-loose');



ids = find(dates_=='01-Jan-2019');
figure()
bardata = [BRs(2:end,ids(1)-T_thres+1:end,1)];
b=bar(dates_(ids:T_full)+calmonths(12), bardata', 'stacked');
for i = 1:7
    b(i).FaceColor = colors(i,:);
    b(i).EdgeColor = 'none';
end
ylabel('Δ Balance of Risk');
hold on;
plot(dates_(ids:T_full)+calmonths(12), pfulls(:,ids(1)-T_thres+1:end,1)-phi0s(:,ids(1)-T_thres+1:end,1), '--','Color','black', 'LineWidth', 2);
grid on;
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "shap_decomposition_loss_plot_DR_clamped_" + string(clamped) + "_" + string(loss_params.alpha) + string(loss_params.beta) + "_" + method + "_" + string(meas)+ ".pdf", '-dpdf', '-painters');
set(gcf,'PaperPositionMode','auto');
print(gcf, "shap_decomposition_loss_plot_DR_clamped_" + string(clamped) + "_" + string(loss_params.alpha) + string(loss_params.beta) + "_" + method + "_" + string(meas)+ ".eps", '-depsc', '-painters','-loose');

figure()
bardata = [BRs(2:end,ids(1)-T_thres+1:end,2)];
b=bar(dates_(ids:T_full)+calmonths(12), bardata', 'stacked');
for i = 1:7
    b(i).FaceColor = colors(i,:);
    b(i).EdgeColor = 'none';
end
ylabel('Δ Balance of Risk');
hold on;
plot(dates_(ids:T_full)+calmonths(12), pfulls(:,ids(1)-T_thres+1:end,2)-phi0s(:,ids(1)-T_thres+1:end,2), '--','Color','black', 'LineWidth', 2);
grid on;
%hL = legend(labels,'NumColumns',4);
% Move the legend to the right side of the figure
%hL.Layout.Tile = 'south';
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "shap_decomposition_loss_plot_QR_clamped_" + string(clamped) + "_" +string(loss_params.alpha) + string(loss_params.beta) + "_" + method + "_" + string(meas)+ ".pdf", '-dpdf', '-painters');
set(gcf,'PaperPositionMode','auto');
print(gcf, "shap_decomposition_loss_plot_QR_clamped_" + string(clamped) + "_" +string(loss_params.alpha) + string(loss_params.beta) + "_" + method + "_" + string(meas)+ ".eps", '-depsc', '-painters','-loose');


%% Balance of risk charts

% standard version
[beta_mixing_sim,beta_kernel_sim,tau_sim] = f_simul_VB(res(sample-T_thres+1,nhor).DR,500);
beta_kernel_sim = permute(beta_kernel_sim,[2,3,1]);
beta_mixing_sim = permute(beta_mixing_sim,[2,3,1]);
tau_sim = tau_sim';

loss_params.lower = 2;
loss_params.upper = 2;
loss_params.alpha = 2;
loss_params.beta  = 2;
loss_params.w     = 0.5;

aux = [];
for i = 1:size(x,1)
    
    for ii = 1:500
        aux.beta_kernel = beta_kernel_sim(:,:,ii);
        aux.beta_mixing = beta_mixing_sim(:,:,ii);
        aux.tau= tau_sim(:,ii);
        RM_ts(i,:,ii) = f_loss_DR(x(i,:),aux,loss_params);
    end
    RM_Qts(i,:,1) = f_loss_QR(x(i,:),squeeze(res(sample-T_thres+1,nhor).QRM.beta(end,:,:)),loss_params);
    
    if i>4
        RM_UCts(i,:,1) = f_loss_TVPSV(1,[res(sample-T_thres+1,nhor).UC.beta(i-4,:) sqrt(res(sample-T_thres+1,nhor).UC.sigmat(i-4,:))],loss_params);
    end
end



% --------- long version
groupColors = get(groot,'defaultAxesColorOrder');
figure()
plot(dates_,squeeze(mean(RM_ts(:,3,:),3)),'LineWidth',1.25);
hold on
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),05,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),95,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',.15,'LineStyle',"None")
plot(dates_,RM_Qts(:,3,1),'--','Color',groupColors(2,:));
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),25,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),75,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',0.25,'LineStyle',"None")
plot(dates_,RM_Qts(:,3,1),'--','Color',groupColors(2,:),'LineWidth',1.25);
hold off
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 20, 8]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "OS_Figure4_RMchart_RiskAverse_long" + meas  + "_.pdf", '-dpdf', '-painters');
set(gcf,'PaperPositionMode','auto');

% ----- for EPS (charts come out pale and lines might not plot, so here's a robust version)
groupColors = get(groot,'defaultAxesColorOrder');
figure()
hold on
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),05,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),95,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',.15,'LineStyle',"None")
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),25,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),75,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',0.45,'LineStyle',"None")
hMean = plot(dates_,squeeze(mean(RM_ts(:,3,:),3)),'Color',groupColors(1,:),'LineWidth',1.25);
hQ    = plot(dates_,RM_Qts(:,3,1),'--','Color',groupColors(2,:),'LineWidth',1.25);
hold off
% Force both lines to front
uistack(hMean, 'top')
uistack(hQ, 'top')
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 20, 8]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "OS_Figure4_RMchart_RiskAverse_long" + meas  + "_.eps", '-depsc', '-painters','-loose');

% --------- short version
ddate = {'01-Jan-2019'};
ids   = find(dates_ == ddate);
groupColors = get(groot,'defaultAxesColorOrder');
figure()
plot(dates_(ids:end),squeeze(mean(RM_ts(ids:end,3,:),3)),'LineWidth',1.25);
hold on
shade(dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),05,3)),'w',...
             dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),95,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',.15,'LineStyle',"None")
plot(dates_(ids:end),RM_Qts(ids:end,3,1),'--','Color',groupColors(2,:));
shade(dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),25,3)),'w',...
             dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),75,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',0.25,'LineStyle',"None")
plot(dates_(ids:end),RM_Qts(ids:end,3,1),'--','Color',groupColors(2,:),'LineWidth',1.25);
hold off
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 10, 8]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "OS_Figure4_RMchart_RiskAverse_short" + meas  + "_.pdf", '-dpdf', '-painters');

% ----- for EPS (charts come out pale and lines might not plot, so here's a robust version)
ddate = {'01-Jan-2019'};
ids   = find(dates_ == ddate);
groupColors = get(groot,'defaultAxesColorOrder');
figure()
hold on
shade(dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),05,3)),'w',...
             dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),95,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',.15,'LineStyle',"None")
shade(dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),25,3)),'w',...
             dates_(ids:end),...
             squeeze(prctile(RM_ts(ids:end,3,:),75,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',0.45,'LineStyle',"None")
hMean = plot(dates_(ids:end),squeeze(mean(RM_ts(ids:end,3,:),3)),'Color',groupColors(1,:),'LineWidth',1.25);
hQ    = plot(dates_(ids:end),RM_Qts(ids:end,3,1),'--','Color',groupColors(2,:),'LineWidth',1.25);
hold off
% Force both lines to front
uistack(hMean, 'top')
uistack(hQ, 'top')
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 20, 8]);
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 10, 8]);

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

set(gcf,'PaperPositionMode','auto');
print(gcf, "OS_Figure4_RMchart_RiskAverse_short" + meas  + "_.eps", '-depsc', '-painters','-loose');

% ---------- Risk neutral risk measure, i.e. with alpha=beta=0
loss_params.lower = 2;
loss_params.upper = 2;
loss_params.alpha = 0;
loss_params.beta  = 0;
loss_params.w     = 0.5;


for i = 1:size(x,1)
    
    for ii = 1:500
        aux.beta_kernel = beta_kernel_sim(:,:,ii);
        aux.beta_mixing = beta_mixing_sim(:,:,ii);
        aux.tau= tau_sim(:,ii);
        RM_ts(i,:,ii) = f_loss_DR(x(i,:),aux,loss_params);
    end
    RM_Qts(i,:,1) = f_loss_QR(x(i,:),squeeze(res(sample-T_thres+1,nhor).QRM.beta(end,:,:)),loss_params);

end

groupColors = get(groot,'defaultAxesColorOrder');
figure()
plot(dates_,squeeze(mean(RM_ts(:,3,:),3)),'LineWidth',1.25);
hold on
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),05,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),95,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',.15,'LineStyle',"None")
plot(dates_,RM_Qts(:,3,1),'--','Color',groupColors(2,:));
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),25,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),75,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',0.25,'LineStyle',"None")
plot(dates_,RM_Qts(:,3,1),'--','Color',groupColors(2,:),'LineWidth',1.25);
hold off
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "RMchart_RiskNeutral_" + meas  + "_.pdf", '-dpdf', '-painters');

% ----- for EPS (charts come out pale and lines might not plot, so here's a robust version)

groupColors = get(groot,'defaultAxesColorOrder');
figure()
hold on
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),05,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),95,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',.15,'LineStyle',"None")
shade(dates_,...
             squeeze(prctile(RM_ts(:,3,:),25,3)),'w',...
             dates_,...
             squeeze(prctile(RM_ts(:,3,:),75,3)),'w',...
             'FillType',[2 1],'FillColor',groupColors(1,:),'FillAlpha',0.45,'LineStyle',"None")
hMean = plot(dates_,squeeze(mean(RM_ts(:,3,:),3)),'Color',groupColors(1,:),'LineWidth',1.25);
hQ    = plot(dates_,RM_Qts(:,3,1),'--','Color',groupColors(2,:),'LineWidth',1.25);
hold off
% Force both lines to front
uistack(hMean, 'top')
uistack(hQ, 'top')
box on
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 20, 8]);
set(gcf,'PaperPositionMode','auto');
print(gcf, "RMchart_RiskNeutral_" + meas  + "_.eps", '-depsc', '-painters','-loose');





