%% Quantile regressions and density matching
% This script estimates quantile regressions of future GDP on current
% values of GDP and the NFCI, and matches skewed t-distributions to the 
% estimated quantiles. It produces the charts presented in the following
% figures:
% - Figure 1. Distribution of GDP growth over time.
% - Figure 5. Predicted Distributions.
% - Figure 6. Median, Interquartile Range, and 5 Percent Quantile of the
%             Predicted Distribution.
% - Figure 7. The Conditional Quantiles and the Skewed t-Distribution.
% - Figure 8. Probability Densities.
% - Figure 9. Growth Entropy and Expected Shortfall over Time.
% - Figure 13. Mean, Volatility, and Entropy.
% - Figure 14. Median, Interquartile Range, and 5 Percent Quantile.
% - Figure A.1. Predictive Distribution of GDP Growth: Location and Scale
%               Parameters over Time.
% - Figure A.2. Predictive Distribution of GDP Growth: Shape and Degrees of
%               Freedom over Time.
% - Figures A.3 and A.4. Predictive Distribution of GDP Growth: Moments
%                        over Time.
% 
% This script requires MATLAB's Optimization Toolbox.
% 
% From the replication files for:
% Tobias Adrian, Nina Boyarchenko, and Domenico Giannone (2018):
% "Vulnerable Growth," American Economic Review.
%% Clear workspace, set file paths and graphics/estimation settings
clear
close all
clc

addpath('azzalini')

set(0, 'defaultAxesFontName', 'Times');
set(0, 'DefaultAxesFontSize',15)
set(0, 'defaultAxesLineStyleOrder', '-|--|:', 'defaultLineLineWidth', 1.5)
setappdata(0, 'defaultAxesXTickFontSize', 1)
setappdata(0, 'defaultAxesYTickFontSize', 1)

% Folder to store figures
FigSubFolder = 'FigMainResults';
if ~exist(FigSubFolder,'dir')
    mkdir(FigSubFolder);
end

% Should saved density matching results be loaded?
% Or should density matching be reperformed?
loadSavedResults = true;

%% Load data and fix forecast settings
load('DataVulnerability.mat', 'X', 'Time', 'Mnem')
% Use 1973Q1-2015Q4 subsample
jtFirst = find(year(Time) == 1973 & month(Time) ==  1);
jtLast  = find(year(Time) == 2015 & month(Time) == 10);
Time = Time(jtFirst:jtLast);
X = X(jtFirst:jtLast, :);
clear('jtFirst', 'jtLast')
[T, n] = size(X);

% Forecast settings
H = [1, 4];          % Horizons to forecast (# of quarters ahead)
QQ = 0.05:0.05:0.95; % Quantiles to estimate in quantile regressions
% Indices for selected quantiles
[~, jq05] = min(abs(QQ - 0.05));
[~, jq25] = min(abs(QQ - 0.25));
[~, jq50] = min(abs(QQ - 0.50));
[~, jq75] = min(abs(QQ - 0.75));
[~, jq95] = min(abs(QQ - 0.95));

% Construct average growth rates
y = X(:, strcmp(Mnem, 'A191RL1Q225SBEA'));
for h = 1:4
    Yh(:, h) = filter(ones(1, h) / h, 1, y);
    Yh(1:(h - 1), h) = NaN;
end

%Dates for which you want to plot the density
JT = [
    find(year(Time) == 2006 & month(Time) ==  4);
    find(year(Time) == 2008 & month(Time) == 10);
    find(year(Time) == 2014 & month(Time) == 10);
    ];

% axis limits for various plots (these are horizon-dependent)
ylimsPredictedDistribution = [
    -15,  20;
	NaN, NaN;
	NaN, NaN;
	-15,  10;
    ];

%% Main results
% Estimate quantile regressions (for all quantiles and forecast horizons)
ResMain    = QRboot(X(:, strcmp(Mnem, 'NFCI')), y, H, 1, QQ);
ResGDPonly = QRboot([], y, H, 1, QQ);
ResUnc     = QRboot([], y, H, 0, QQ);

% Loop over forecast horizons
for h = H
    % Get quantiles of predictive distribution
    YQunc     = ResUnc.YQ(end, :, h);
    YQ        = ResMain.YQ(:, :, h);
    Y2        = ResMain.Y2(:, h);
    YQGDPonly = ResGDPonly.YQ(:, :, h);
    QQ        = ResMain.QQ;
    
    %% Figure 5: Predicted Distributions
    filename = fullfile(FigSubFolder, ['GDP_Q', num2str(h), '_TS.pdf']);
    if h == 1
        legendLocation = 'NorthEast';
    elseif h == 4
        legendLocation = 'SouthEast';
    end
    PlotPredictiveTS(YQ, QQ, Time, Yh(:, h),...
                     'ylims', ylimsPredictedDistribution(h, :),...
                     'filename', filename,...
                     'legendLocation', legendLocation);
    clear('filename', 'legendLocation')
    
    %% Figure 6: Median, Interquartile Range, and 5 Percent Quantile
    % (a)/(b) Median and Interquartile Range
    f = figure;
    scatter(YQ(:, jq75) - YQ(:, jq25), YQ(:, jq50), 75, [20, 91, 142] / 255)
    xlim([1, 8])
    ylim([-3, 5])
    box on;
    xlabel('Interquartile range')
    ylabel('Median')
    filename = fullfile(FigSubFolder, ['IQRvM_Q', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    % (c)/(d) Q5 and Interquartile Range
    f = figure;
    scatter(YQ(:, jq75) - YQ(:, jq25), YQ(:, jq05), 75, [20, 91, 142] / 255)
    box on;
    xlabel('Interquartile range')
    ylabel('5th quantile')
    xlim([1, 8])
    ylim([-15, 5])
    filename = fullfile(FigSubFolder, ['IQRvQ5_Q', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')

    %% Fit skewed-t densities to estimated quantiles OR load saved results
    if loadSavedResults
        filename = ['ResMatch_H', num2str(h), '.mat'];
        disp(['Loading saved density matching results from file ', filename])
        load(filename)
        clear('filename')
    else
        disp('Matching densities for quantile regression using GDP and NFCI.')
        ResMatch = Step2match(YQ, YQunc, QQ);
        disp('Matching densities for quantile regression using GDP only.')
        ResMatchGDPonly = Step2match(YQGDPonly, YQunc, QQ);

        % Save results to .mat file
        filename = ['ResMatch_H', num2str(h), '.mat'];
        disp(['Saving results to file ', filename])
        save(filename, 'ResMatch', 'ResMatchGDPonly')
        clear('filename')
    end
    
    %% Figure 1: 3D density plot
    f = figure;
    meshc(Time, ResMatch.YY', ResMatch.PST')
    datetick('x', 'yyyy')
    view(85, 50)
    set(gca, 'YLim', [-15, 10])
    set(gca, 'XLim', [Time(1), Time(end)])
    xlabel('Time')
    filename = fullfile(FigSubFolder, ['Dens3D_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    %% Figure A.1: Location and Scale Parameters over Time
    % (a)/(b) Location
    f = figure;
    plot(Time, ResMatch.STpar(:, 1))
    ylabel('Location')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([-5, 5])
    filename = fullfile(FigSubFolder, ['Location_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    % (b)/(c) Scale
    f = figure;
    plot(Time, ResMatch.STpar(:, 2))
    ylabel('Scale')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([1, 7])
    filename = fullfile(FigSubFolder, ['Scale_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename');
    
    %% Figure A.2: Shape and Degrees of Freedom over Time
    % (a)/(b) Shape
    f = figure;
    plot(Time, ResMatch.STpar(:, 3))
    ylabel('Shape')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([-5, 20])
    filename = fullfile(FigSubFolder, ['Shape_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    % (c)/(d) Degrees of Freedom
    f = figure;
    plot(Time, ResMatch.STpar(:, 4))
    ylabel('Degrees of Freedom')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([0, 30])
    filename = fullfile(FigSubFolder, ['DF_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    %% Figure A.3: Moments over Time (Mean and Variance)
    % (a)/(b) Mean
    f = figure;
    plot(Time, ResMatch.MeanST)
    ylabel('Mean')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([-4, 6])
    filename = fullfile(FigSubFolder, ['Mean_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    % (c)/(d) Variance
    f = figure;
    plot(Time, ResMatch.VarianceST)
    ylabel('Variance')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([0, 30])
    filename = fullfile(FigSubFolder, ['Variance_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    %% Figure A.4: Moments over Time (Skewness and Kurtosis)
    % (a)/(b) Skewness
    f = figure;
    plot(Time, ResMatch.SkewnessST)
    ylabel('Skewness')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([-2.5, 1.5])
    hline = refline(0, 0);
    hline.Color = ones(1,3) * 0.25;
    hline.LineWidth = 0.5;
    uistack(hline, 'bottom')
    filename = fullfile(FigSubFolder, ['Skewness_H', num2str(h),'.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    % (c)/(d) Kurtosis
    f = figure;
    plot(Time, ResMatch.KurtosisST)
    ylabel('Kurtosis')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([2, 20])
    hline = refline(0, 3);
    hline.Color = ones(1, 3) * 0.25;
    hline.LineWidth = 0.5;
    uistack(hline, 'bottom')
    filename = fullfile(FigSubFolder, ['Kurtosis_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    %% Plot densities for specific dates (Figures 7 and 8)
    for jt = JT'
        % Figure 8: Probability Densities
        f = figure;
        plot(ResMatch.YY, ResMatch.PST(jt + h, :),...
             ResMatch.YY, ResMatchGDPonly.PST(jt + h, :), '--');
        axis([-15, 15, 0, 0.4])
        ylabel('PDF')
        legend('GDP and NFCI', 'GDP only', 'Location', 'NorthWest')
        grid on
        filename = fullfile(FigSubFolder, ['PDF_H', num2str(h), '_', num2str(year(Time(jt))), 'Q', num2str(ceil(month(Time(jt)) / 3)), '.pdf']);
        printpdf(f, filename);
        clear('f', 'filename');
        
        % Figure 9: Conditional Quantiles and the Skewed t-Distribution
        f = figure;
        plot(QQ, ResMatch.QST(jt + h, :),...
             QQ, ResMatchGDPonly.QST(jt + h, :), '--',...
             QQ, YQ(jt + h, :), '-*');
        xlabel('\tau')
        ylabel('GDP growth')
        ylim([-15, 10])
        grid on
        legend('GDP and NFCI', 'GDP only', 'Raw', 'Location', 'NorthWest')
        grid on
        filename = fullfile(FigSubFolder, ['InverseCDF_H', num2str(h), '_', num2str(year(Time(jt))), 'Q', num2str(ceil(month(Time(jt)) / 3)), '.pdf']);
        printpdf(f, filename);
        clear('f', 'filename')
    end

    %% Figure 9: Growth Entropy and Expected Shortfall over Time
    % (a)/(b) Entropy
    f = figure;
    hold on
    plot(Time, ResMatch.LeftEntropyST, '-');
    plot(Time, ResMatch.RightEntropyST, '--');
    hold off
    legend('Downside', 'Upside')
    ylabel('Relative Entropy')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([-0.5, 1.5])
    box on
    filename = fullfile(FigSubFolder, ['Entropy_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')

    % (c)/(d) Expected Shortfall and Longrise
    f = figure;
    hold on
    plot(Time, ResMatch.EsfST, '-')
    plot(Time, ResMatch.EljST, '--')
    hold off
    legend('Shortfall', 'Longrise')
    datetick('x', 'yyyy')
    xlim([Time(1), Time(end)])
    ylim([-20, 20])
    box on
    filename = fullfile(FigSubFolder, ['ExpectedShortfall_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    %% Figure 13: Mean, Volatility, and Entropy
    % (a)/(b) Mean and Volatility
    f=figure;
    scatter(ResMatch.MeanST, sqrt(ResMatch.VarianceST), 75, [20, 91, 142] / 255)
    box on;
    ylabel('Standard Deviation')
    xlabel('Mean')
    ylim([1, 6])
    yticks(1:6)
    xlim([-5, 6])
    xticks((-5):6)
    filename = fullfile(FigSubFolder, ['MEANvSTD_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')

    % (c)/(d) Mean and Entropy
    f = figure;
    scatter(ResMatch.MeanST, ResMatch.LeftEntropyST, 75, [20, 91, 142] / 255)
    box on;
    ylabel('Downside Entropy')
    xlabel('Mean')
    ylim([-0.5, 1.5])
    xlim([-5, 6])
    xticks((-5):6)
    filename = fullfile(FigSubFolder, ['MEANvENT_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')

    %% Figure 14: Median, Interquartile Range, and 5 Percent Quantile
    % (a)/(b) Cross-Correlogram Mean and Volatility
    nLags = 16;
    [XFC, Lags, Bounds] = crosscorr(ResMatch.MeanST((h + 1):end), sqrt(ResMatch.VarianceST((h + 1):end)), nLags);
    % Determine significant/insignificant correlations
    ix = ((-nLags):nLags) + nLags + 1;
    ix_sig = ix((XFC >= Bounds(1)) | (XFC <= Bounds(2))); % Correlations outside of significance bounds
    ix_ins = ix((Bounds(1) > XFC) & (Bounds(2) < XFC));
    
    f = figure;
    hold on
    stem(Lags(ix_sig), XFC(ix_sig), 'filled', 'Color', [20, 91, 142] / 255, 'LineWidth', 2, 'MarkerSize', 10)
    stem(Lags(ix_ins), XFC(ix_ins), 'Color',[20,91,142]/255, 'LineWidth',2, 'MarkerSize',10)
    hold off
    box on
    xlim([-nLags, nLags])
    xticks(-nLags:4:nLags)
    ylim([-1.0, 0.2])
    xlabel('Volatility Lag')
    filename = fullfile(FigSubFolder, ['MeanStdxcorr_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')

    % (c)/(d) Cross-Correlogram Mean and Entropy
    nLags=16;
    [XFC, Lags, Bounds] = crosscorr(ResMatch.MeanST((h + 1):end), ResMatch.LeftEntropyST((h + 1):end), nLags);
    % Determine significant/insignificant correlations
    ix = ((-nLags):nLags) + nLags + 1;
    ix_sig = ix((XFC >= Bounds(1)) | (XFC <= Bounds(2)));
    ix_ins = ix((Bounds(1) > XFC) & (Bounds(2) < XFC));
    
    f = figure;
    hold on
    stem(Lags(ix_sig), XFC(ix_sig), 'filled', 'Color', [20, 91, 142] / 255, 'LineWidth', 2, 'MarkerSize', 10)
    stem(Lags(ix_ins), XFC(ix_ins), 'Color', [20, 91, 142] / 255, 'LineWidth', 2, 'MarkerSize', 10)
    hold off
    box on
    xlim([-nLags, nLags])
    xticks(-nLags:4:nLags)
    ylim([-1.0, 0.2])
    xlabel('Entropy Lag')
    filename = fullfile(FigSubFolder, ['MeanEntrxcorr_H', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
end
    