%% Estimate quantile regression coefficients and bootstrap bands
% This script produces charts of estimated quantile regression coefficients
% with bootstrapped confidence bands for the following figure:
% - Figure 4. Estimated Quantile Regression Coefficients.
% 
% 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 = 'FigInference';
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

% 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

% Bootstrap settings
nDrawsBoot = 1000; % number of bootstrap replications
nLagsBootVAR = 4;  % number of lags to use in approximating VAR

%% Estimate quantile regressions, generate bootstrapped confidence bands
% Quantile regressions include both GDP and NFCI
Res = QRboot(X(:, strcmp(Mnem, 'NFCI')), y, H, 1, QQ,...
             nDrawsBoot, nLagsBootVAR);

% Plot results for coefficients on NFCI and GDP
for j = 2:3
    % j denotes the column index within the matrix of regressors for the
    % regressor of interest (either NFCI or GDP)
    switch j
        case 2
            regressorName = 'NFCI';
            ylims = [-3.0, 1.5];
        case 3
            regressorName = 'GDP';
            ylims = [-0.2, 0.5];
    end

    for h = H
        BQ  = Res.BQ(j, :, h);
        B2  = Res.B2(j, h);
        bBQ = squeeze(Res.bBQ(j, :, h, :));
        bB2 = squeeze(Res.bB2(j, h, :));
        QQ = Res.QQ;
        filename = fullfile(FigSubFolder, ['GDPv', regressorName, '_coefs_Q', num2str(h), '.pdf']);
        PlotQRbands(BQ, bBQ, B2, bB2, QQ,...
                    'ylims', ylims, 'filename', filename);
    end
end