%% C1_make_paper_outputs.m
clear; clc;

% Paths
root = fileparts(pwd);
inT  = fullfile(root,'Data','thispaper_exports');
inK  = fullfile(root,'Data','kaenzig_exports');
outTables = fullfile(root,'Output','tables');
outFigs   = fullfile(root,'Output','figures');
mkdir(outTables); mkdir(outFigs);

addpath(fullfile(pwd,'C1_lib'));

% Filenames
fnT_shocks = fullfile(inT,'thispaper_baseline_structural_shocks.csv');
fnT_dfse   = fullfile(inT,'thispaper_baseline_ngml_df_and_se.csv');
fnT_corrK  = fullfile(inT,'thispaper_baseline_corr_with_kaenzig_news.csv');

fnT_irf.h    = fullfile(inT,'thispaper_baseline_irf_horizon.csv');
fnT_irf.pt   = fullfile(inT,'thispaper_baseline_irf_point.csv');
fnT_irf.lo68 = fullfile(inT,'thispaper_baseline_irf_lo68.csv');
fnT_irf.hi68 = fullfile(inT,'thispaper_baseline_irf_hi68.csv');
fnT_irf.lo90 = fullfile(inT,'thispaper_baseline_irf_lo90.csv');
fnT_irf.hi90 = fullfile(inT,'thispaper_baseline_irf_hi90.csv');
fnT_hd       = fullfile(inT,'thispaper_baseline_hd_point.csv');

fnT_extcorr  = fullfile(inT,'thispaper_extensions_shock_correlations.csv');

fnK_irf.h    = fullfile(inK,'kaenzig_baseline_irf_horizon.csv');
fnK_irf.pt   = fullfile(inK,'kaenzig_baseline_irf_point.csv');
fnK_irf.lo68 = fullfile(inK,'kaenzig_baseline_irf_lo68.csv');
fnK_irf.hi68 = fullfile(inK,'kaenzig_baseline_irf_hi68.csv');
fnK_irf.lo90 = fullfile(inK,'kaenzig_baseline_irf_lo90.csv');
fnK_irf.hi90 = fullfile(inK,'kaenzig_baseline_irf_hi90.csv');
fnK_varnames = fullfile(inK,'kaenzig_baseline_varnames.csv');
fnK_hd       = fullfile(inK,'kaenzig_baseline_hd.csv');

fnT_B  = fullfile(inT,'thispaper_baseline_B.csv');
fnK_b1 = fullfile(inK,'kaenzig_baseline_b1.csv');

%
% load once (baseline shocks + df/se)
Eraw = readmatrix(fnT_shocks);
E    = Eraw(:, end-5:end);                 % last 6 columns = shocks (works for 6 or 7 cols)

Mdf  = readmatrix(fnT_dfse);
df   = Mdf(:,2);
dfse = Mdf(:,3);

%
% 1) TABLE 1
JB = zeros(6,1); P = zeros(6,1);
for j=1:6
    [JB(j), P(j)] = jb_test_simple(E(:,j));
end
JB = round(JB,3); P = round(P,3); df = round(df,3); dfse = round(dfse,3);

disp(' ');
disp('Table 1: Assessment of Non-Gaussianity of Structural Shocks');
fprintf('%-6s %12s %12s %12s %12s\n','Shock','JB','p-value','df','df_se');
for j=1:6
    fprintf('eps_%d %12.3f %12.3f %12.3f %12.3f\n', j, JB(j), P(j), df(j), dfse(j));
end
writetable(table((1:6)',JB,P,df,dfse,'VariableNames',{'shock','JB','p_value','df','df_se'}), ...
    fullfile(outTables,'Table1_JB_df.csv'));

%
% 2) TABLE 2
T2 = readtable(fnT_corrK);
n = string(T2{:,1});
v = round(T2{:,2},3);
disp(' ');
disp('Table 2: Correlations of Non-Gaussian Structural Shocks with Kaenzigs Shock');
fprintf('%8s','');
fprintf('%10s',n); fprintf('\n');
fprintf('%8s','Corr');
fprintf('%10.3f',v); fprintf('\n');
writetable(array2table(v','VariableNames',cellstr(n)), ...
           fullfile(outTables,'Table2_Corr_with_Kaenzig_news.csv'));

%
% 3) FIGURE 1 (IRFs)
IRF_T = read_irf_bundle(fnT_irf);
IRF_K = read_irf_bundle(fnK_irf);
try
    vt = readtable(fnK_varnames,'ReadVariableNames',false);
    varnames = string(vt{:,1});
catch
    varnames = strings(0,1);
end
plot_irf_compare(IRF_T, IRF_K, varnames, fullfile(outFigs,'Figure1_IRF_compare.pdf'));

%
% 4) FIGURE 2 (HD)

% This paper (header + ISO date)
T = readtable(fnT_hd);
HD_T.t = datenum(datetime(T.date,'InputFormat','yyyy-MM-dd'));
HD_T.y = T.hd_point;

% Kaenzig (decimal year -> first of month)
HD_K = read_hd_series(fnK_hd);
x  = HD_K.t(:);
dt = datetime(floor(x), 1, 1) + calmonths(round((x-floor(x))*12));
HD_K.t = datenum(dt);

plot_hd_compare(HD_T, HD_K, fullfile(outFigs,'Figure2_HD_compare.pdf'));

%
% 5) MATRICES
disp(' ');
disp('=== Baseline B matrix (this paper) ===');
B = readmatrix(fnT_B);
for i = 1:size(B,1)
fprintf('%8.3f ', B(i,:));
fprintf('\n');
end
disp('=== Kaenzig baseline b1'' ===');
b1 = readmatrix(fnK_b1);
fprintf('%8.3f\n', b1(:));

%
% 6) TABLE A1
A1 = readtable(fnT_extcorr);
namesA1 = string(A1{:,1});
valsA1  = round(A1{:,2},3);

disp(' ');
disp('Table A1: Correlations of Non-Gaussian Structural Shocks with Kaenzigs Shocks from Extended Models');
fprintf('%-30s %10s\n','Extension','Corr');
fprintf('%-30s %10s\n',repmat('-',1,9),repmat('-',1,4));
for i=1:numel(valsA1)
    fprintf('%-30s %10.3f\n', namesA1(i), valsA1(i));
end
A1{:,2} = valsA1;
writetable(A1, fullfile(outTables,'TableA1_Extensions_Corr.csv'));

disp(' ');
disp('=== C1 finished ===');
