%% Out-of-sample forecasting exercise for alternative models
% This script performs an out-of-sample forecasting exercise assessing the
% performance of two alternative econometric approaches (a conditionally
% Gaussian model and a flexible nonparametric approach) and produces charts
% for the following figure:
% - Figure A.8. Out-of-Sample Performance of Alternative Approaches.
% 
% NOTE: this script should be run separately for one quarter and four
% quarter ahead forecast horizons by manually changing the variable 'h'.
% 
% 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; 
clear
close all
clc

%% Set forecast horizon (run script separately for h = 1 and h = 4)
h = 1;

%% Set file paths and graphics/estimation settings
addpath('azzalini')
addpath('Tatevik')
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 = 'FigOOS';
if ~exist(FigSubFolder,'dir')
    mkdir(FigSubFolder);
end

% Should saved recursive estimation results and out-of-sample forecasts be
% loaded?
% Or should estimation and out-of-sample forecasting 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 (NOTE: h is set in the section above)
QQ = 0.05:0.05:0.95;   % Quantiles to estimate in quantile regressions
deltaYY = 0.1;         % Step size and grid to use for computing CDFs/PDFs
YY = (-20):deltaYY:20;
% Index of start date for out-of-sample forecasting
jtFirstOOS = find(year(Time) == 1993 & month(Time) == 1); 
% 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 rate and matrices of regressors
y = X(:, strcmp(Mnem, 'A191RL1Q225SBEA'));
yh = filter(ones(h, 1) / h, 1, y);
yh(1:(h - 1)) = NaN;
Z = [ones(size(y)), X(:, strcmp(Mnem, 'NFCI')), y];
ZGDPonly = [ones(size(y)), y];

%% In-sample and out-of-sample estimation/forecasting OR load saved results
if loadSavedResults
    filename = ['ResAlternativeOOS_H', num2str(h), '.mat'];
    disp(['Loading saved estimation and forecasting results from file ', filename])
    load(filename)
    clear('filename')
else
    %% Initialize matrices to store forecasts
    % Raw quantiles
    YQ_NaNs = NaN(length(Time),length(QQ));
    YQGDPonly_IS  = YQ_NaNs;
    YQGDPonly_OOS = YQ_NaNs;
    YQunc_IS    = YQ_NaNs;
    YQunc_OOS   = YQ_NaNs;
    YQml_IS       = YQ_NaNs;
    YQml_OOS      = YQ_NaNs;
    YQnp_IS       = YQ_NaNs;
    YQnp_OOS      = YQ_NaNs;
    clear('YQ_NaNs')
    
    % PDFs (evaluated over grid)
    P_NaNs = NaN(length(Time),length(YY));
    PSTGDPonly_IS  = P_NaNs;
    PSTGDPonly_OOS = P_NaNs;
    PSTunc_IS    = P_NaNs;
    PSTunc_OOS   = P_NaNs;
    Pml_IS         = P_NaNs;
    Pml_OOS        = P_NaNs;
    PSTnp_IS       = P_NaNs;
    PSTnp_OOS      = P_NaNs;
    clear('P_NaNs')

    % Smoothed quantiles
    Q_NaNs = NaN(length(Time),length(QQ));
    QSTGDPonly_IS  = Q_NaNs;
    QSTGDPonly_OOS = Q_NaNs;
    QSTunc_IS    = Q_NaNs;
    QSTunc_OOS   = Q_NaNs;
    Qml_IS         = Q_NaNs;
    Qml_OOS        = Q_NaNs;
    QSTnp_IS       = Q_NaNs;
    QSTnp_OOS      = Q_NaNs;
    clear('Q_NaNs')

    % CDFs (evaluated over grid)
    C_NaNs = NaN(length(Time),length(YY));
    CSTGDPonly_IS  = C_NaNs;
    CSTGDPonly_OOS = C_NaNs;
    CSTunc_IS    = C_NaNs;
    CSTunc_OOS   = C_NaNs;
    Cml_IS         = C_NaNs;
    Cml_OOS        = C_NaNs;
    CSTnp_IS       = C_NaNs;
    CSTnp_OOS      = C_NaNs;
    clear('C_NaNs')

    % Skewed t-distribution parameters
    STpar_NaNs = NaN(length(Time),4);
    STparGDPonly_IS  = STpar_NaNs;
    STparGDPonly_OOS = STpar_NaNs;
    STparunc_IS    = STpar_NaNs;
    STparunc_OOS   = STpar_NaNs;
    STparnp_IS       = STpar_NaNs;
    STparnp_OOS      = STpar_NaNs;
    clear('STpar_NaNs')

    % Predictive scores
    Score_NaNs = NaN(length(Time),1);
    ScoreSTGDPonly_IS  = Score_NaNs;
    ScoreSTGDPonly_OOS = Score_NaNs;
    ScoreSTunc_IS    = Score_NaNs;
    ScoreSTunc_OOS   = Score_NaNs;
    Scoreml_IS         = Score_NaNs;
    Scoreml_OOS        = Score_NaNs;
    ScoreSTnp_IS       = Score_NaNs;
    ScoreSTnp_OOS      = Score_NaNs;
    clear('Score_NaNs')

    % Probability integral transforms
    Pit_NaNs = NaN(length(Time),1);
    PitSTGDPonly_IS  = Pit_NaNs;
    PitSTGDPonly_OOS = Pit_NaNs;
    PitSTunc_IS    = Pit_NaNs;
    PitSTunc_OOS   = Pit_NaNs;
    Pitml_IS         = Pit_NaNs;
    Pitml_OOS        = Pit_NaNs;
    PitSTnp_IS       = Pit_NaNs;
    PitSTnp_OOS      = Pit_NaNs;
    clear('Pit_NaNs')

    %% In-sample estimation of conditional quantiles
    for jq = 1:length(QQ)
        % Quantile regression with GDP only
        bGDPonly = rq(ZGDPonly(1:(end - h), :), yh((h + 1):end), QQ(jq));
        YQGDPonly_IS((h + 1):end, jq) = ZGDPonly(1:(end - h), :) * bGDPonly;
        % Unconditional quantiles (quantile regression on constant)
        bunc = rq(ones(size(yh((h + 1):end))), yh((h + 1):end), QQ(jq));
        YQunc_IS((h + 1):end, jq) = ones(length(Time) - h, 1) * bunc;
    end
    % Conditionally Gaussian model
    [gamma_IS, delta_IS] = EstimateCGM(Z(1:(end - h), :), yh((h + 1):end));
    YQml_IS((h + 1):end, :) = ComputeCGMCondQuantiles(Z(1:(end - h), :), gamma_IS, delta_IS, QQ);
    % Nonparametric model (need to exclude constant)
    YQnp_IS((h + 1):end, :) = ComputeNonparCondQuantiles(Z(1:(end - h), 2:end), yh((h + 1):end), QQ);

    %% Fit skewed-t distribution for in-sample unconditional quantiles
    qqTarg = YQunc_IS(end, :);
    [lc, sc, sh, df] = QuantilesInterpolation(qqTarg, QQ);
    PSTunc_IS((h + 1):end, :)     = repmat(dskt(YY, lc, sc, sh, df), length(Time) - h, 1);
    QSTunc_IS((h + 1):end, :)     = repmat(qskt(QQ, lc, sc, sh, df), length(Time) - h, 1);
    CSTunc_IS((h + 1):end, :)     = repmat(pskt(YY, lc, sc, sh, df), length(Time) - h, 1);
    STparunc_IS((h + 1):end, :)   = repmat([lc, sc, sh, df], length(Time) - h, 1);
    ScoreSTunc_IS((h + 1):end, :) = dskt(yh((h + 1):end), lc, sc, sh, df);
    PitSTunc_IS((h + 1):end, :)   = pskt(yh((h + 1):end), lc, sc, sh, df);

    %% Fit Gaussian distributions in-sample
    Pml_IS((h + 1):end, :)     = ComputeCGMCondDensity(Z(1:(end - h),:), YY, gamma_IS, delta_IS);
    Qml_IS                     = YQml_IS; % no skewed-t fitting of quantiles
    Cml_IS((h + 1):end, :)     = ComputeCGMCondCDF(Z(1:(end - h), :), YY, gamma_IS, delta_IS);
    Scoreml_IS(1:(end - h), :) = ComputeCGMCondDensity(Z(1:(end - h), :), yh((h + 1):end), gamma_IS, delta_IS);
    Pitml_IS((h + 1):end, :)   = ComputeCGMCondCDF(Z(1:(end - h), :), yh((h + 1):end), gamma_IS, delta_IS);

    %% Fit t-densities for in-sample and out-of-sample estimation
    for jt = 1:(length(Time) - h)
        if (month(Time(jt)) == 1) && (jt >= jtFirstOOS)
            disp(['Computing in-sample and out-of-sample predictive densities in ', num2str(year(Time(jt)))])
        else
            disp(['Computing in-sample predictive densities in ', num2str(year(Time(jt)))])
        end
        
        yhRealized = yh(jt + h, :); % realized value for average GDP growth

        % Fit skewed-t density for quantile regression with GDP only, in-sample fit
        qqTarg_GDPonly = YQGDPonly_IS(jt + h, :);
        [lc, sc, sh, df] = QuantilesInterpolation(qqTarg_GDPonly, QQ);
        PSTGDPonly_IS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
        QSTGDPonly_IS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
        CSTGDPonly_IS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
        STparGDPonly_IS(jt + h, :)   = [lc, sc, sh, df];
        ScoreSTGDPonly_IS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
        PitSTGDPonly_IS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);

        % Fit skewed-t density for nonparametric model with GDP and NFCI
        qqTarg_np = YQnp_IS(jt + h, :);
        [lc, sc, sh, df] = QuantilesInterpolation(qqTarg_np, QQ);
        PSTnp_IS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
        QSTnp_IS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
        CSTnp_IS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
        STparnp_IS(jt + h, :)   = [lc, sc, sh, df];
        ScoreSTnp_IS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
        PitSTnp_IS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);

        if jt >= jtFirstOOS
            for jq = 1:length(QQ)
                % Quantile regression with GDP only (+constant), out-of-sample
                bGDPonly = rq(ZGDPonly(1:(jt - h), :), yh((h + 1):jt, :), QQ(jq));
                YQGDPonly_OOS(jt + h, jq) = ZGDPonly(jt, :) * bGDPonly;
            end
            
            % Conditionally Gaussian model, out-of-sample
            if jt == jtFirstOOS
                % First OOS forecast: use default initial conditions
                [gamma_OOS, delta_OOS] = EstimateCGM(Z(1:(jt - h), :), yh((h + 1):jt));
            else
                % Subsequent OOS forecasts: use last period estimates as initial conditions
                [gamma_OOS, delta_OOS] = EstimateCGM(Z(1:(jt - h), :),yh((h + 1):jt), gamma_OOS, delta_OOS);
            end
            YQml_OOS(jt + h, :) = ComputeCGMCondQuantiles(Z(jt, :), gamma_OOS, delta_OOS, QQ);
            
            % Nonparametric model, out-of-sample
            % (no constant is used, so the first column of Z containing the
            % constant is dropped)
            YQnp_OOS(jt+h,:) = ComputeNonparCondQuantiles(Z(1:(jt - h), 2:end), yh((h + 1):jt), QQ, Z(jt ,2:end));

            % Fit skewed-t distribution for quantile regression with GDP only,
            % out-of-sample
            qqTarg = YQGDPonly_OOS(jt + h, :);
            [lc, sc, sh, df] = QuantilesInterpolation(qqTarg, QQ);
            PSTGDPonly_OOS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
            QSTGDPonly_OOS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
            CSTGDPonly_OOS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
            STparGDPonly_OOS(jt + h, :)   = [lc, sc, sh, df];
            ScoreSTGDPonly_OOS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
            PitSTGDPonly_OOS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);

            % Fit conditionally Gaussian model, out-of-sample
            Pml_OOS(jt + h, :)     = ComputeCGMCondDensity(Z(jt,:),YY,gamma_OOS,delta_OOS);
            Qml_OOS(jt + h, :)     = YQml_OOS(jt+h,:); % no skewed-t fitting of quantiles for ML
            Cml_OOS(jt + h, :)     = ComputeCGMCondCDF(Z(jt, :), YY, gamma_OOS, delta_OOS);
            Scoreml_OOS(jt + h, :) = ComputeCGMCondDensity(Z(jt, :), yhRealized, gamma_OOS, delta_OOS);
            Pitml_OOS(jt + h, :)   = ComputeCGMCondCDF(Z(jt, :), yhRealized, gamma_OOS, delta_OOS);
            
            % Fit skewed t-distribution for nonparametric model, out-of-sample
            qqTarg = YQnp_OOS(jt + h, :);
            [lc, sc, sh, df] = QuantilesInterpolation(qqTarg, QQ);
            PSTnp_OOS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
            QSTnp_OOS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
            CSTnp_OOS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
            STparnp_OOS(jt + h, :)   = [lc, sc, sh, df];
            ScoreSTnp_OOS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
            PitSTnp_OOS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);
        end
    end

    % Save results for quick generation of figures
    filename = ['ResAlternativeOOS_H', num2str(h), '.mat'];
    disp(['Saving results to file ', filename])
    save(filename,...
         'YQGDPonly_IS',      'YQGDPonly_OOS',      'YQml_IS',    'YQml_OOS',    'YQnp_IS',      'YQnp_OOS',...
         'PSTGDPonly_IS',     'PSTGDPonly_OOS',     'Pml_IS',     'Pml_OOS',     'PSTnp_IS',     'PSTnp_OOS',...
         'QSTGDPonly_IS',     'QSTGDPonly_OOS',     'Qml_IS',     'Qml_OOS',     'QSTnp_IS',     'QSTnp_OOS',...
         'CSTGDPonly_IS',     'CSTGDPonly_OOS',     'Cml_IS',     'Cml_OOS',     'CSTnp_IS',     'CSTnp_OOS',...
         'STparGDPonly_IS',   'STparGDPonly_OOS',   'STparnp_IS', 'STparnp_OOS',...
         'ScoreSTGDPonly_IS', 'ScoreSTGDPonly_OOS', 'Scoreml_IS', 'Scoreml_OOS', 'ScoreSTnp_IS', 'ScoreSTnp_OOS',...
         'PitSTGDPonly_IS',   'PitSTGDPonly_OOS',   'Pitml_IS',   'Pitml_OOS',   'PitSTnp_IS',   'PitSTnp_OOS')
    clear('filename')
end

%% Figure A.8. Out-of-Sample Performance of Alternative Approaches.
% (a) Quantiles: MLE
f = figure;
plot(Time, YQml_OOS(:, [jq05, jq50, jq95]), '-b',...
     Time, YQml_IS(:, [jq05, jq50, jq95]),'-.k');
datetick('x','yyyy');
if h == 1
    axis([Time(1), Time(end), -15, 20])
elseif h == 4
    axis([Time(1), Time(end), -15, 10])
end
filename = fullfile(FigSubFolder, ['QuantOOS_ml_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% (b) Quantiles: Non-parametric
f = figure;
plot(Time, YQnp_OOS(:, [jq05, jq50, jq95]), '-b',...
     Time, YQnp_IS(:, [jq05, jq50, jq95]), '-.k');
datetick('x', 'yyyy');
if h == 1
    axis([Time(1), Time(end), -15, 20])
elseif h == 4
    axis([Time(1), Time(end), -15, 10])
end
filename = fullfile(FigSubFolder, ['QuantOOS_np_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% (c) Scores: MLE
f = figure;
plot(Time, Scoreml_OOS, Time, ScoreSTGDPonly_OOS, '--');
datetick('x', 'yyyy')
axis([Time(1), Time(end), 0.0, 0.4])
ylabel('Scores')
legend('Maximum likelihood', 'Quantile regression, GDP only', 'Location', 'NorthWest')
filename = fullfile(FigSubFolder, ['ScoresOOS_ml_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% (d) Scores: Non-parametric
f = figure;
plot(Time, ScoreSTnp_OOS, Time, ScoreSTGDPonly_OOS, '--');
datetick('x', 'yyyy')
axis([Time(1), Time(end), 0.0, 0.4])
ylabel('Scores')
legend('Nonparametric', 'Quantile regression, GDP only', 'Location', 'NorthWest')
filename = fullfile(FigSubFolder, ['ScoresOOS_np_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% (e)/(f): PITs (MLE and Non-parametric)
% The code below was modified from files provided by Barbara Rossi and
% Tatevik Sekhposyan implementing the specification tests for predictive
% densities described in Rossi and Sekhposyan (2017).
rvec = 0:0.001:1;
zSTGDPonly_ecdf = PITtest(PitSTGDPonly_OOS, rvec);
zSTnp_ecdf = PITtest(PitSTnp_OOS, rvec);
zml_ecdf = PITtest(Pitml_OOS, rvec);
if h == 1
    % Use asymptotic 5% critical value from Rossi and Sekhposyan (2017): 1.34
    kappa_GDPonly = 1.34;
    kappa_ml      = 1.34;
    kappa_np      = 1.34;
elseif h == 4
    % Compute bootstrapped 5% critical values
    PITs = [PitSTGDPonly_OOS, Pitml_OOS, PitSTnp_OOS];
    PITs = PITs((jtFirstOOS + h):end, :); % remove NaNs
    for i = 1:size(PITs, 2)
        [~, critvalues] = rstestboot(PITs(:, i));
        testcritvalues(:, :, i) = round(critvalues .* 100) ./ 100;
    end
    kappa_GDPonly = testcritvalues(1, 2, 1);
    kappa_ml      = testcritvalues(1, 2, 2);
    kappa_np      = testcritvalues(1, 2, 3);
end

% Plot PIT for conditionally Gaussian model vs. quantile regression with GDP only
f = figure;
plot(rvec, zml_ecdf, '-b', rvec, zSTGDPonly_ecdf, '--r')
hold on
P = sum(isnan(Pitml_OOS) == 0);
plot(rvec,[rvec', rvec' - (kappa_ml / sqrt(P)), rvec' + (kappa_ml / sqrt(P))], '-b');
plot(rvec,[rvec', rvec' - (kappa_GDPonly / sqrt(P)), rvec' + (kappa_GDPonly / sqrt(P))], '--r');
hold off
xlim([0, 1])
ylim([0, 1])
xlabel('\tau')
ylabel('Empirical CDF')
grid on
legend('Maximum likelihood', 'Quantile regression, GDP only', 'Theoretical and 5% Critical Values',...
       'Location','NorthWest');
filename = fullfile(FigSubFolder, ['PITS_ml_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% Plot PIT for nonparametric model vs. quantile regression with GDP only
f = figure;
plot(rvec, zSTnp_ecdf, '-b', rvec, zSTGDPonly_ecdf, '--r')
hold on
P = sum(isnan(PitSTnp_OOS)==0);
plot(rvec,[rvec', rvec' - (kappa_np / sqrt(P)), rvec' + (kappa_np / sqrt(P))], '-b');
plot(rvec,[rvec', rvec' - (kappa_GDPonly / sqrt(P)), rvec' + (kappa_GDPonly / sqrt(P))], '--r');
hold off
xlim([0, 1])
ylim([0, 1])
xlabel('\tau')
ylabel('Empirical CDF')
grid on
legend('Nonparametric', 'Quantile regression, GDP only', 'Theoretical and 5% Critical Values',...
       'Location','NorthWest');
filename = fullfile(FigSubFolder, ['PITS_np_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')