clear;
clc;
close all;

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

load('BH_shocks');

supply = table2timetable(T(:,1:2));
supply = retime(supply,"Quarterly",'mean');


demand = table2timetable(T2(:,1:3));
demand = retime(demand,"Quarterly");

% this file contains decompositions computed for energy prices instead of
% aggregated commodity prices by applying a different grouping in the
% density decomposition functions
load('loss_decomposition_for_BH_comparison.mat')

for k = 1:numSamples
    BRs(:,k,1) = RM(k).DR.phi(:,3);
    BRs(:,k,2) = RM(k).QR.phi(:,3);
end

tab = rmmissing(synchronize(synchronize(synchronize(array2timetable(BRs(:,:,1)',"RowTimes",dates_(T_thres:end)),supply),demand)));
corr1 = corr(table2array(tab(:,["Var4","Var3","oilSupplyShocks","economicActivityShocks","oilConsumptionDemandShocks"])));

figure()
y = zscore(table2array(tab(:,["Var4"])));
x = -1.*zscore(table2array(tab(:,["oilSupplyShocks"])));
t = datenum(tab.Properties.RowTimes);

% Fit a line y = a*x + b
p = polyfit(x, y, 1);          % 1 = linear
yFit = polyval(p, x);

% Plot the regression line on top of scatter
hold on;
plot(x, yFit, 'r-', 'LineWidth', 1.5);  % red line, thicker for visibility

scatter(x,y, 60, t, 'filled'); 
colormap(viridis); 
cb=colorbar; 
% Relabel colorbar ticks with datetime strings
cb.Ticks = linspace(min(t), max(t), 5); % choose 5 ticks
cb.TickLabels = datestr(cb.Ticks, 'yyyy-mm-dd');
% Convert tick positions back to datetime
tickDates = datetime(cb.Ticks, 'ConvertFrom', 'datenum');
% Add zero reference lines
xline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
yline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
% Format as "YYYYQn"
tickLabels = strcat(string(year(tickDates)), "Q", string(quarter(tickDates)));
cb.TickLabels = tickLabels;

% Annotate points far from the center
dist = sqrt(x.^2 + y.^2);
%dist = max(abs([x, y]), [], 2);  % max of abs(x) and abs(y)
threshold = prctile(dist, 95); % top 5% farthest points
idxFar = find(dist > threshold);% Get axis limits

[A,B] = sort(dist,'descend');
idxFar = B(1:4);
xl = xlim(gca);
yl = ylim(gca);

% Axis range for offset
dx = 0.02 * diff(xl);  % 2% of x-axis range
dy = 0.02 * diff(yl);  % 2% of y-axis range
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
for i = idxFar
    lbl = strcat(string(year(t(i))), "Q", string(quarter(t(i))));
    
    % Offset label
    xLabelPos = x(i) + dx;
    yLabelPos = y(i) + dy;
    
    % Clamp inside axes
    xLabelPos = min(max(xLabelPos, xl(1)), xl(2));
    yLabelPos = min(max(yLabelPos, yl(1)), yl(2));
    
    % Add text
    text(xLabelPos, yLabelPos, lbl, 'FontSize',10, 'FontWeight','bold', 'Color','k');
end
% Compute axis limits
xl = [min(x), max(x)];
yl = [min(y), max(y)];

% Expand limits by 5% on each side
xMargin = 2;
yMargin = 2;
xlim(xl + [-xMargin, xMargin]);
ylim(yl + [-yMargin, yMargin]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
% Save as vector PDF (preferred for LaTeX and journals)
print(gcf, "Main_Figure6_scatter_comparison_oil_DR.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure6_scatter_comparison_oil_DR.eps", '-depsc', '-painters','-loose');



y = zscore(table2array(tab(:,["Var3"])));
x = zscore(table2array(tab(:,["economicActivityShocks"])));
t = datenum(tab.Properties.RowTimes);

figure()
% Fit a line y = a*x + b
p = polyfit(x, y, 1);          % 1 = linear
yFit = polyval(p, x);

% Plot the regression line on top of scatter
hold on;
plot(x, yFit, 'r-', 'LineWidth', 1.5);  % red line, thicker for visibility

scatter(x,y, 60, t, 'filled'); 
colormap(viridis); 
cb=colorbar; 
% Relabel colorbar ticks with datetime strings
cb.Ticks = linspace(min(t), max(t), 5); % choose 5 ticks
cb.TickLabels = datestr(cb.Ticks, 'yyyy-mm-dd');
% Convert tick positions back to datetime
tickDates = datetime(cb.Ticks, 'ConvertFrom', 'datenum');
% Add zero reference lines
xline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
yline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
% Format as "YYYYQn"
tickLabels = strcat(string(year(tickDates)), "Q", string(quarter(tickDates)));
cb.TickLabels = tickLabels;

% Annotate points far from the center
dist = sqrt(x.^2 + y.^2);
%dist = max(abs([x, y]), [], 2);  % max of abs(x) and abs(y)
threshold = prctile(dist, 95); % top 5% farthest points
idxFar = find(dist > threshold);% Get axis limits

[A,B] = sort(dist,'descend');
idxFar = B(1:4);
xl = xlim(gca);
yl = ylim(gca);

% Axis range for offset
dx = 0.02 * diff(xl);  % 2% of x-axis range
dy = 0.02 * diff(yl);  % 2% of y-axis range
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
for i = idxFar
    lbl = strcat(string(year(t(i))), "Q", string(quarter(t(i))));
    
    % Offset label
    xLabelPos = x(i) + dx;
    yLabelPos = y(i) + dy;
    
    % Clamp inside axes
    xLabelPos = min(max(xLabelPos, xl(1)), xl(2));
    yLabelPos = min(max(yLabelPos, yl(1)), yl(2));
    
    % Add text
    text(xLabelPos, yLabelPos, lbl, 'FontSize',10, 'FontWeight','bold', 'Color','k');
end
% Compute axis limits
xl = [min(x), max(x)];
yl = [min(y), max(y)];

% Expand limits by 5% on each side
xMargin = 2;
yMargin = 2;
xlim(xl + [-xMargin, xMargin]);
ylim(yl + [-yMargin, yMargin]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
% Save as vector PDF (preferred for LaTeX and journals)
print(gcf, "Main_Figure6_scatter_comparison_activity_DR.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure6_scatter_comparison_activity_DR.eps", '-depsc', '-painters','-loose');



tab = rmmissing(synchronize(synchronize(synchronize(array2timetable(BRs(:,:,2)',"RowTimes",dates_(T_thres:end)),supply),demand),oil));
corr2 = corr(table2array(tab(:,["Var4","Var3","oilSupplyShocks","economicActivityShocks","oilConsumptionDemandShocks","Oil"])));


figure()
y = zscore(table2array(tab(:,["Var4"])));
x = -1.*zscore(table2array(tab(:,["oilSupplyShocks"])));
t = datenum(tab.Properties.RowTimes);

% Fit a line y = a*x + b
p = polyfit(x, y, 1);          % 1 = linear
yFit = polyval(p, x);

% Plot the regression line on top of scatter
hold on;
plot(x, yFit, 'r-', 'LineWidth', 1.5);  % red line, thicker for visibility

scatter(x,y, 60, t, 'filled'); 
colormap(viridis); 
cb=colorbar; 
% Relabel colorbar ticks with datetime strings
cb.Ticks = linspace(min(t), max(t), 5); % choose 5 ticks
cb.TickLabels = datestr(cb.Ticks, 'yyyy-mm-dd');
% Convert tick positions back to datetime
tickDates = datetime(cb.Ticks, 'ConvertFrom', 'datenum');
% Add zero reference lines
xline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
yline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
% Format as "YYYYQn"
tickLabels = strcat(string(year(tickDates)), "Q", string(quarter(tickDates)));
cb.TickLabels = tickLabels;

% Annotate points far from the center
dist = sqrt(x.^2 + y.^2);
%dist = max(abs([x, y]), [], 2);  % max of abs(x) and abs(y)
threshold = prctile(dist, 95); % top 5% farthest points
idxFar = find(dist > threshold);% Get axis limits

[A,B] = sort(dist,'descend');
idxFar = B(1:4);
xl = xlim(gca);
yl = ylim(gca);

% Axis range for offset
dx = 0.02 * diff(xl);  % 2% of x-axis range
dy = 0.02 * diff(yl);  % 2% of y-axis range
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
for i = idxFar
    lbl = strcat(string(year(t(i))), "Q", string(quarter(t(i))));
    
    % Offset label
    xLabelPos = x(i) + dx;
    yLabelPos = y(i) + dy;
    
    % Clamp inside axes
    xLabelPos = min(max(xLabelPos, xl(1)), xl(2));
    yLabelPos = min(max(yLabelPos, yl(1)), yl(2));
    
    % Add text
    text(xLabelPos, yLabelPos, lbl, 'FontSize',10, 'FontWeight','bold', 'Color','k');
end
% Compute axis limits
xl = [min(x), max(x)];
yl = [min(y), max(y)];

% Expand limits by 5% on each side
xMargin = 2;
yMargin = 2;
xlim(xl + [-xMargin, xMargin]);
ylim(yl + [-yMargin, yMargin]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
% Save as vector PDF (preferred for LaTeX and journals)
print(gcf, "Main_Figure6_scatter_comparison_oil_QR.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure6_scatter_comparison_oil_QR.eps", '-depsc', '-painters','-loose');



y = zscore(table2array(tab(:,["Var3"])));
x = zscore(table2array(tab(:,["economicActivityShocks"])));
t = datenum(tab.Properties.RowTimes);

figure()
% Fit a line y = a*x + b
p = polyfit(x, y, 1);          % 1 = linear
yFit = polyval(p, x);

% Plot the regression line on top of scatter
hold on;
plot(x, yFit, 'r-', 'LineWidth', 1.5);  % red line, thicker for visibility

scatter(x,y, 60, t, 'filled'); 
colormap(viridis); 
cb=colorbar; 
% Relabel colorbar ticks with datetime strings
cb.Ticks = linspace(min(t), max(t), 5); % choose 5 ticks
cb.TickLabels = datestr(cb.Ticks, 'yyyy-mm-dd');
% Convert tick positions back to datetime
tickDates = datetime(cb.Ticks, 'ConvertFrom', 'datenum');
% Add zero reference lines
xline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
yline(0, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
% Format as "YYYYQn"
tickLabels = strcat(string(year(tickDates)), "Q", string(quarter(tickDates)));
cb.TickLabels = tickLabels;

% Annotate points far from the center
dist = sqrt(x.^2 + y.^2);
%dist = max(abs([x, y]), [], 2);  % max of abs(x) and abs(y)
threshold = prctile(dist, 95); % top 5% farthest points
idxFar = find(dist > threshold);% Get axis limits

[A,B] = sort(dist,'descend');
idxFar = B(1:4);
xl = xlim(gca);
yl = ylim(gca);

% Axis range for offset
dx = 0.02 * diff(xl);  % 2% of x-axis range
dy = 0.02 * diff(yl);  % 2% of y-axis range
% Resize figure for publication (e.g., 16 cm wide × 8 cm tall)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 16, 8]);
for i = idxFar
    lbl = strcat(string(year(t(i))), "Q", string(quarter(t(i))));
    
    % Offset label
    xLabelPos = x(i) + dx;
    yLabelPos = y(i) + dy;
    
    % Clamp inside axes
    xLabelPos = min(max(xLabelPos, xl(1)), xl(2));
    yLabelPos = min(max(yLabelPos, yl(1)), yl(2));
    
    % Add text
    text(xLabelPos, yLabelPos, lbl, 'FontSize',10, 'FontWeight','bold', 'Color','k');
end
% Compute axis limits
xl = [min(x), max(x)];
yl = [min(y), max(y)];

% Expand limits by 5% on each side
xMargin = 2;
yMargin = 2;
xlim(xl + [-xMargin, xMargin]);
ylim(yl + [-yMargin, yMargin]);

% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
% Save as vector PDF (preferred for LaTeX and journals)
print(gcf, "Main_Figure6_scatter_comparison_activity_QR.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure6_scatter_comparison_activity_QR.eps", '-depsc', '-painters','-loose');

