%% Plot raw data, create scatterplots for univariate quantile regressions
% This script produces charts for the descriptive analysis presented in the
% following figures:
% - Figure 2. Raw Data.
% - Figure 3. Quantile Regressions.
% The charts are saved in the folder 'FigExploratory'.
% 
% 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 settings
clear
close all
clc

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 = 'FigExploratory';
if ~exist(FigSubFolder, 'dir')
    mkdir(FigSubFolder);
end

%% 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')

% 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));
[~, jq50] = min(abs(QQ - 0.50));
[~, jq95] = min(abs(QQ - 0.95));

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

%% Figure 2: Raw Data (GDP and NFCI)
f = figure;
ax = gca;

yyaxis left
plot(Time, X(:, strcmp(Mnem, 'A191RL1Q225SBEA')))
ax.YColor = 'k';
ylabel('Real GDP growth (annualized)');
datetick('x', 'yyyy')

yyaxis right
plot(Time, X(:, strcmp(Mnem, 'NFCI')), '--')
ax.YColor = 'k';
ylabel('NFCI (index)');

legend('GDP', 'NFCI', 'Location', 'NorthEast')
filename = fullfile(FigSubFolder, 'GDP_NFCI.pdf');
printpdf(f, filename);
clear('f', 'filename')

%% Figure 3: Quantile Regressions (univariate)
% y-axis limits for scatter plots (these are horizon-dependent)
ylimsScatter = [
    -15,  20;
    NaN, NaN;
    NaN, NaN;
    -15,  10;
];

for h = H
    % GDP on GDP
    Z = [ones(size(y)), y];
    yh = Yh((h + 1):end, h);
    Zh = Z(1:(end - h), :);
    B2 = Zh \ yh;
    Y2 = [ones(h, 1) * NaN; Zh * B2];
    for jq = 1:length(QQ)
        BQ(:, jq) = rq(Zh, yh, QQ(jq));
        YQ(:, jq) = [ones(h, 1) * NaN; Zh * BQ(:, jq)];
    end
    
    f=figure;
    scatter(Zh(:, 2), yh, 50);
    box on
    hold on
    p1 = plot(Zh(:, 2), YQ((h + 1):end, jq05), '-k',...
        Zh(:, 2), YQ((h + 1):end, jq50), '-r',...
        Zh(:, 2), YQ((h + 1):end, jq95), '-g',...
        Zh(:, 2), Y2((h + 1):end), '-b');
    hold off
    axis tight
    ylim(ylimsScatter(h, :))
    legend(p1, 'Q5', 'Q50', 'Q95', 'OLS', 'Location', 'SouthEast');
    if h == 1
        ylabel('GDP growth 1 quarter ahead')
    else
        ylabel(['GDP growth ', num2str(h), ' quarters ahead'])
    end
    xlabel('Current quarter GDP growth')
    filename = fullfile(FigSubFolder, ['GDPvGDP_Q', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
    
    % GDP on NFCI
    Z = [ones(size(y)), X(:, strcmp(Mnem, 'NFCI'))];
    yh = Yh((h + 1):end, h);
    Zh = Z(1:(end - h), :);
    B2 = Zh \ yh;
    Y2 = [ones(h, 1) * NaN; Zh * B2];
    for jq = 1:length(QQ)
        BQ(:, jq) = rq(Zh, yh, QQ(jq));
        YQ(:, jq) = [ones(h, 1) * NaN; Zh * BQ(:, jq)];
    end
    
    f=figure;
    scatter(Zh(:, 2), yh, 50);
    box on
    hold on
    p1 = plot(Zh(:, 2), YQ((h + 1):end, jq05), '-k',...
        Zh(:, 2), YQ((h + 1):end, jq50), '-r',...
        Zh(:, 2), YQ((h + 1):end, jq95), '-g',...
        Zh(:, 2), Y2((h + 1):end), '-b');
    hold off
    axis tight
    ylim(ylimsScatter(h, :))
    legend(p1, 'Q5', 'Q50', 'Q95', 'OLS', 'Location', 'SouthWest');
    if h == 1
        ylabel('GDP growth 1 quarter ahead')
    else
        ylabel(['GDP growth ', num2str(h), ' quarters ahead'])
    end
    xlabel('NFCI')
    filename = fullfile(FigSubFolder, ['GDPvNFCI_Q', num2str(h), '.pdf']);
    printpdf(f, filename);
    clear('f', 'filename')
end