clear;
clc;
close all;

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

%% Load SPF microdata for Core PCE
T = readtable('SPFmicrodata.xlsx','Sheet','COREPCE');
T.Timestamp = datetime(T.YEAR, 3*T.QUARTER, 1); 

% select 4 quarter ahead qoq growth 
var = 'COREPCE6';

% extract ids and timestamps
allIDs   = unique(T.ID);
allTimes = unique(T.Timestamp);

% pivot the table for IDs in columns and time in rows
pivotTbl = unstack(T(:, {'Timestamp','ID',var}), var, 'ID');
pivotTbl = outerjoin(table(allTimes,'VariableNames',{'Timestamp'}),...
pivotTbl, 'Keys','Timestamp', 'MergeKeys',true);
for j = 2:width(pivotTbl)
    % Convert empty cells to 'NaN' and all entries to string
    pivotTbl.(j) = string(pivotTbl.(j));      % converts empty cells to ""
    pivotTbl.(j)(pivotTbl.(j)=="") = '#N/A';  % replace empty strings with "NaN"
end
pivotTbl(:,2:end) = varfun(@(x) str2double(strrep(string(x),'#N/A','NaN')), pivotTbl(:,2:end));

% convert data columns to numeric
for j = 2:width(pivotTbl)   % start at 2, skip the first column
    col = pivotTbl.(j);

    if iscell(col)  % often comes in as cell array of strings
        col = str2double(strrep(col, '#N/A', 'NaN'));
    elseif isstring(col) || ischar(col)
        col = str2double(strrep(string(col), '#N/A', 'NaN'));
    end

    pivotTbl.(j) = col;
end

% convert to timetable
SPFdataPCE = table2timetable(pivotTbl);


%% Load SPF microdata for CPI
T = readtable('SPFmicrodata.xlsx','Sheet','CPI');
T.Timestamp = datetime(T.YEAR, 3*T.QUARTER, 1); 

% select 4 quarter ahead qoq growth 
var = 'CPI6';

% extract ids and timestamps
allIDs   = unique(T.ID);
allTimes = unique(T.Timestamp);

% pivot the table for IDs in columns and time in rows
pivotTbl = unstack(T(:, {'Timestamp','ID',var}), var, 'ID');
pivotTbl = outerjoin(table(allTimes,'VariableNames',{'Timestamp'}),...
pivotTbl, 'Keys','Timestamp', 'MergeKeys',true);
for j = 2:width(pivotTbl)
    % Convert empty cells to 'NaN' and all entries to string
    pivotTbl.(j) = string(pivotTbl.(j));      % converts empty cells to ""
    pivotTbl.(j)(pivotTbl.(j)=="") = '#N/A';  % replace empty strings with "NaN"
end
pivotTbl(:,2:end) = varfun(@(x) str2double(strrep(string(x),'#N/A','NaN')), pivotTbl(:,2:end));

% convert data columns to numeric
for j = 2:width(pivotTbl)   % start at 2, skip the first column
    col = pivotTbl.(j);

    if iscell(col)  % often comes in as cell array of strings
        col = str2double(strrep(col, '#N/A', 'NaN'));
    elseif isstring(col) || ischar(col)
        col = str2double(strrep(string(col), '#N/A', 'NaN'));
    end

    pivotTbl.(j) = col;
end

% convert to timetable
SPFdataCPI = table2timetable(pivotTbl);

%% compute metrics

timestamps = SPFdataPCE.Properties.RowTimes;        % keep time labels
XPCE = table2array(SPFdataPCE);             % numeric forecast panel
XCPI = table2array(SPFdataCPI);             % numeric forecast panel

nT = size(XPCE,1);

% Preallocate result arrays
IQR          = nan(nT,2);
MAD          = nan(nT,2);
STD          = nan(nT,2);
Skew_moment  = nan(nT,2);
BowleySkew   = nan(nT,2);
Kurt_moment  = nan(nT,2);
MoorsKurt    = nan(nT,2);

for t = 1:nT
    xPCE = XPCE(t,:);  % forecasts at time t
    xPCE = xPCE(~isnan(xPCE)); 

    xCPI = XCPI(t,:);  % forecasts at time t
    xCPI = xCPI(~isnan(xCPI)); 
    if all(isnan(xPCE)) && all(isnan(xCPI))
        continue  % skip if no data
    end
    
    % ---- Dispersion ----
    q25 = quantile(xPCE,0.25);
    q75 = quantile(xPCE,0.75);
    IQR(t,1) = q75 - q25;
    
    med = median(xPCE);
    MAD(t,1) = median(abs(xPCE - med));
    
    STD(t,1) = std(xPCE);
    
    % ---- Skewness ----
    Skew_moment(t,1) = skewness(xPCE);
    
    q50 = quantile(xPCE,0.50);
    BowleySkew(t,1) = (q75 + q25 - 2*q50) / (q75 - q25);
    
    % ---- Kurtosis ----
    Kurt_moment(t,1) = kurtosis(xPCE);
    
    P12_5 = quantile(xPCE,0.125);
    P37_5 = quantile(xPCE,0.375);
    P62_5 = quantile(xPCE,0.625);
    P87_5 = quantile(xPCE,0.875);
    
    MoorsKurt(t,1) = ((P87_5 - P62_5) + (P37_5 - P12_5)) / (q75 - q25);
    
    % ---- Dispersion ----
    q25 = quantile(xCPI,0.25);
    q75 = quantile(xCPI,0.75);
    IQR(t,2) = q75 - q25;
    
    med = median(xCPI);
    MAD(t,2) = median(abs(xCPI - med));
    
    STD(t,2) = std(xCPI);
    
    % ---- Skewness ----
    Skew_moment(t,2) = skewness(xCPI);
    
    q50 = quantile(xCPI,0.50);
    BowleySkew(t,2) = (q75 + q25 - 2*q50) / (q75 - q25);
    
    % ---- Kurtosis ----
    Kurt_moment(t,2) = kurtosis(xCPI);
    
    P12_5 = quantile(xPCE,0.125);
    P37_5 = quantile(xPCE,0.375);
    P62_5 = quantile(xPCE,0.625);
    P87_5 = quantile(xPCE,0.875);
    
    MoorsKurt(t,2) = ((P87_5 - P62_5) + (P37_5 - P12_5)) / (q75 - q25);
end

%% load forecasting results

load("PCEPILFE_H_5_priorM_2_priorK_0_cluster_.mat");


dates_(T_thres+12:end);
startid = find(timestamps-calmonths(2)==dates_(T_thres+12));
endid   = find(timestamps-calmonths(2)==dates_(end));

zscor_xnan = @(x) bsxfun(@rdivide, bsxfun(@minus, x, mean(x,'omitnan')), std(x, 'omitnan'));

meas = IQR(startid:endid,:);
% backcast PCE with CPI
aux1(:,1) = rescale(movmean(meas(:,1),[1,1]),0,1);
aux1(:,2) = rescale(movmean(meas(:,2),[1,1]),0,1);
aux2 = rescale(movvar(y(T_thres+12:end),[1,1]),0,1);

% --- Step 1: Find first valid row of the series to extend ---
firstRow = find(~isnan(aux1(:,1)), 1, 'first');
validGuide = aux1(1:firstRow,2);           % only from first valid row
changes = flipud(diff(flipud(validGuide)));                % changes between consecutive periods

% --- Step 3: Backward extension ---
% Reverse and negate the changes for backward reconstruction
backward_changes = flipud(cumsum(flipud(changes)));

% Starting value for extension
startValue = aux1(firstRow,1);

% Reconstruct backward values
backward_values = startValue + (backward_changes);
aux1(:,1) = [backward_values; aux1(firstRow:end,1)];


colors = [0 0.4470 0.7410;...	        %Blue
          0.8500 0.3250 0.0980;...	    %Orange
          0.9290 0.6940 0.1250];	        %Yellow

qs = squeeze(QScore10(:,4,[1,2,3]));

figure()
yyaxis right
plot(dates_(T_thres+12:end),cumsum(qs),"Color",colors(1,:))
ax = gca;
ax.YColor = colors(1,:);                  % Left axis color
yyaxis left
hold on
area(dates_(T_thres+12:end),aux1(:,1),'FaceColor',[0.4980, 0.4980, 0.4980],'EdgeColor','none','FaceAlpha',0.2)
ylim([0 1])
plot(dates_(T_thres+12:end),aux2,"--",'Color','black')
hold off
ax = gca;
ax.YColor = 'black';                  % Left axis color
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 11, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "Main_Figure2_CQS_10" + "_.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure2_CQS_10" + "_.eps", '-depsc', '-painters','-loose');


qs = squeeze(QScore50(:,4,[1,2,3]));

figure()
yyaxis right
plot(dates_(T_thres+12:end),cumsum(qs),"Color",colors(2,:))
ax = gca;
ax.YColor = colors(2,:);                  % Left axis color
yyaxis left
hold on
area(dates_(T_thres+12:end),aux1(:,1),'FaceColor',[0.4980, 0.4980, 0.4980],'EdgeColor','none','FaceAlpha',0.2)
ylim([0 1])
plot(dates_(T_thres+12:end),aux2,"--",'Color','black')
hold off
ax = gca;
ax.YColor = 'black';                  % Left axis color
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 11, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "Main_Figure2_CQS_50" + "_.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure2_CQS_50" + "_.eps", '-depsc', '-painters','-loose');


qs = squeeze(QScore90(:,4,[1,2,3]));

figure()
yyaxis right
plot(dates_(T_thres+12:end),cumsum(qs),"Color",colors(3,:))
ax = gca;
ax.YColor = colors(3,:);                  % Left axis color
yyaxis left
hold on
area(dates_(T_thres+12:end),aux1(:,1),'FaceColor',[0.4980, 0.4980, 0.4980],'EdgeColor','none','FaceAlpha',0.2)
ylim([0 1])
plot(dates_(T_thres+12:end),aux2,"--",'Color','black')
hold off
ax = gca;
ax.YColor = 'black';                  % Left axis color
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Units', 'centimeters', 'Position', [2, 2, 11, 8]);
% Set renderer to 'painters' for vector output (PDF/EPS)
set(gcf, 'Renderer', 'painters');
print(gcf, "Main_Figure2_CQS_90" + "_.pdf", '-dpdf', '-painters');
print(gcf, "Main_Figure2_CQS_90" + "_.eps", '-depsc', '-painters','-loose');


