%% Estimate conditional quantiles using alternative econometric approaches
% This script estimates conditional quantiles of GDP growth using two
% alternative econometric approaches (a conditionally Gaussian model and 
% a flexible nonparametric approach) and produces charts for the following
% figure:
% - Figure 12. Alternative Econometric Approaches: Predicted Distributions.
% 
% This script requires MATLAB's Optimization and Global Optimization
% Toolboxes.
% 
% 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('cgm')
addpath('nonpar')

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

%% Figure 12. Alternative Econometric Approaches: Predicted Distributions.
% Use current quarter NFCI and GDP growth as predictors
Z = [ones(length(Time), 1), X(:, strcmp(Mnem, 'NFCI')), y];

for h = H
    yh = Yh((h + 1):end, h);
    Zh = Z(1:(end - h), :);
    
    % Conditionally Gaussian model (MLE)
    [gammahat, deltahat] = EstimateCGM(Zh, yh);
    YQ_ml = NaN(length(Time), length(QQ));
    YQ_ml((h + 1):end, :) = ComputeCGMCondQuantiles(Zh, gammahat, deltahat, QQ);
    
    % Nonparametric
    YQ_np = NaN(length(Time), length(QQ));
    YQ_np((h + 1):end, :) = ComputeNonparCondQuantiles(Zh(:, 2:end), yh, QQ);
        
    % x- and y-axis limits for plots
    xlims = [Time(1), Time(end)];
    if h == 1
        ylims = [-15, 20];
    elseif h == 4
        ylims = [-15, 10];
    end
    
    % (a) MLE
    filename = fullfile(FigSubFolder, ['ml_Q', num2str(h), '.pdf']);
    PlotPredictiveTS(YQ_ml, QQ, Time, Yh(:, h),...
                     'ylims', ylims, 'filename', filename);
    % (b) Non-parametric
    filename = fullfile(FigSubFolder, ['npqreg_Q', num2str(h), '.pdf']);
    PlotPredictiveTS(YQ_np, QQ, Time, Yh(:, h),...
                     'ylims', ylims, 'filename', filename);
end