%% Out-of-sample forecasting exercise
% This script estimates the models presented in the paper in psuedo
% real-time, and compares out-of-sample forecasting performance. It
% produces the charts presented in the following figures:
% - Figure 10. Out-of-Sample Predictions.
% - Figure 11. Out-of-Sample Accuracy.
% 
% 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 Toolbox.
% 
% 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')

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 = ['ResOOS_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));
    YQ_IS         = YQ_NaNs;
    YQ_OOS        = YQ_NaNs;
    YQGDPonly_IS  = YQ_NaNs;
    YQGDPonly_OOS = YQ_NaNs;
    YQunc_IS      = YQ_NaNs;
    YQunc_OOS     = YQ_NaNs;
    clear('YQ_NaNs')
    
    % PDFs (evaluated over grid)
    P_NaNs         = NaN(length(Time),length(YY));
    PST_IS         = P_NaNs;
    PST_OOS        = P_NaNs;
    PSTGDPonly_IS  = P_NaNs;
    PSTGDPonly_OOS = P_NaNs;
    PSTunc_IS      = P_NaNs;
    PSTunc_OOS     = P_NaNs;
    clear('P_NaNs')
    
    % Smoothed quantiles
    Q_NaNs = NaN(length(Time),length(QQ));
    QST_IS         = Q_NaNs;
    QST_OOS        = Q_NaNs;
    QSTGDPonly_IS  = Q_NaNs;
    QSTGDPonly_OOS = Q_NaNs;
    QSTunc_IS      = Q_NaNs;
    QSTunc_OOS     = Q_NaNs;
    clear('QST_NaNs')
    
    % CDFs (evaluated over grid)
    C_NaNs = NaN(length(Time),length(YY));
    CST_IS          = C_NaNs;
    CST_OOS         = C_NaNs;
    CSTGDPonly_IS   = C_NaNs;
    CSTGDPonly_OOS  = C_NaNs;
    CSTunc_IS       = C_NaNs;
    CSTunc_OOS      = C_NaNs;
    clear('C_NaNs')
    
    % Skewed t-distribution parameters
    STpar_NaNs = NaN(length(Time),4);
    STpar_IS         = STpar_NaNs;
    STpar_OOS        = STpar_NaNs;
    STparGDPonly_IS  = STpar_NaNs;
    STparGDPonly_OOS = STpar_NaNs;
    STparunc_IS      = STpar_NaNs;
    STparunc_OOS     = STpar_NaNs;
    clear('STpar_NaNs')
    
    % Predictive scores
    Score_NaNs = NaN(length(Time),1);
    ScoreST_IS         = Score_NaNs;
    ScoreST_OOS        = Score_NaNs;
    ScoreSTGDPonly_IS  = Score_NaNs;
    ScoreSTGDPonly_OOS = Score_NaNs;
    ScoreSTunc_IS      = Score_NaNs;
    ScoreSTunc_OOS     = Score_NaNs;
    clear('Score_NaNs')
    
    % Probability integral transforms
    Pit_NaNs = NaN(length(Time),1);
    PitST_IS         = Pit_NaNs;
    PitST_OOS        = Pit_NaNs;
    PitSTGDPonly_IS  = Pit_NaNs;
    PitSTGDPonly_OOS = Pit_NaNs;
    PitSTunc_IS      = Pit_NaNs;
    PitSTunc_OOS     = Pit_NaNs;
    clear('Pit_NaNs')
    
    % Left entropy
    Entropy_NaNs = NaN(length(Time),1);
    LeftEntropy_IS  = Entropy_NaNs;
    LeftEntropy_OOS = Entropy_NaNs;
    clear('Entropy_NaNs')
    
    %% In-sample estimation of conditional quantiles
    for jq = 1:length(QQ)
        % Quantile regression with both NFCI and GDP
        b = rq(Z(1:(end - h), :), yh((h + 1):end), QQ(jq));
        YQ_IS((h + 1):end, jq) = Z(1:(end-h),:) * b;
        % 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
    
    %% 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 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 NFCI and GDP, in-sample
        qqTarg = YQ_IS(jt + h,:);
        [lc,sc,sh,df] = QuantilesInterpolation(qqTarg, QQ);
        PST_IS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
        QST_IS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
        CST_IS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
        STpar_IS(jt + h, :)   = [lc, sc, sh, df];
        ScoreST_IS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
        PitST_IS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);
        
        Temp = PST_IS(jt + h, :) .* (YY < QST_IS(jt + h, jq50));
        LeftEntropy_IS(jt + h, :) = -sum((log(PSTunc_IS(jt + h, :)) - log(PST_IS(jt + h, :))) .* Temp * deltaYY);
        
        % Fit skewed-t density for quantile regression with GDP only, in-sample
        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);
        
        if jt >= jtFirstOOS
            if month(Time(jt))==1
                disp(['Now computing the real-time predictive densities in ', num2str(year(Time(jt)))])
            end
            
            for jq = 1:length(QQ)
                % Quantile regression with both NFCI and GDP, out-of-sample
                b = rq(Z(1:(jt - h), :), yh(h + 1:jt, :), QQ(jq));
                YQ_OOS(jt + h, jq) = Z(jt, :) * b;
                % Quantile regression with GDP only, out-of-sample
                bGDPonly = rq(ZGDPonly(1:(jt - h), :), yh((h + 1):jt, :), QQ(jq));
                YQGDPonly_OOS(jt + h, jq) = ZGDPonly(jt, :) * bGDPonly;
                % Unconditional quantiles, out-of-sample
                bunc = rq(ones(size(yh((h + 1):jt, :))), yh((h + 1):jt, :), QQ(jq));
                YQunc_OOS((jt + h), jq) = bunc;
            end
            
            % Fit skewed-t distribution for quantile regression with NFCI and
            % GDP, out-of-sample
            qqTarg = YQ_OOS(jt + h, :);
            [lc, sc, sh, df] = QuantilesInterpolation(qqTarg, QQ);
            PST_OOS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
            QST_OOS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
            CST_OOS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
            STpar_OOS(jt + h, :)   = [lc, sc, sh, df];
            ScoreST_OOS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
            PitST_OOS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);
            
            % 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 skewed t-distribution for unconditional quantiles,
            % out-of-sample
            qqTarg = YQunc_OOS(jt + h, :);
            [lc, sc, sh, df]  = QuantilesInterpolation(qqTarg, QQ);
            PSTunc_OOS(jt + h, :)     = dskt(YY, lc, sc, sh, df);
            QSTunc_OOS(jt + h, :)     = qskt(QQ, lc, sc, sh, df);
            CSTunc_OOS(jt + h, :)     = pskt(YY, lc, sc, sh, df);
            STparunc_OOS(jt + h, :)   = [lc, sc, sh, df];
            ScoreSTunc_OOS(jt + h, :) = dskt(yhRealized, lc, sc, sh, df);
            PitSTunc_OOS(jt + h, :)   = pskt(yhRealized, lc, sc, sh, df);
            
            % Compute entropy for skewed t-distribution from quantile
            % regression with GDP and NFCI, out-of-sample
            Temp = PST_OOS(jt + h, :) .* (YY < QST_OOS(jt + h, jq50));
            LeftEntropy_OOS(jt + h, :) = -sum((log(PSTunc_OOS(jt + h, :)) - log(PST_OOS(jt + h, :))) .* Temp * deltaYY);
        end
    end
    
    % Save results for quick generation of figures
    filename = ['ResOOS_H', num2str(h), '.mat'];
    disp(['Saving results to file ', filename])
    save(filename,...
        'YQ_IS',      'YQ_OOS',      'YQGDPonly_IS',      'YQGDPonly_OOS',      'YQunc_IS',      'YQunc_OOS',...
        'PST_IS',     'PST_OOS',     'PSTGDPonly_IS',     'PSTGDPonly_OOS',     'PSTunc_IS',     'PSTunc_OOS',...
        'QST_IS',     'QST_OOS',     'QSTGDPonly_IS',     'QSTGDPonly_OOS',     'QSTunc_IS',     'QSTunc_OOS',...
        'CST_IS',     'CST_OOS',     'CSTGDPonly_IS',     'CSTGDPonly_OOS',     'CSTunc_IS',     'CSTunc_OOS',...
        'STpar_IS',   'STpar_OOS',   'STparGDPonly_IS',   'STparGDPonly_OOS',   'STparunc_IS',   'STparunc_OOS',...
        'ScoreST_IS', 'ScoreST_OOS', 'ScoreSTGDPonly_IS', 'ScoreSTGDPonly_OOS', 'ScoreSTunc_IS', 'ScoreSTunc_OOS',...
        'PitST_IS',   'PitST_OOS',   'PitSTGDPonly_IS',   'PitSTGDPonly_OOS',   'PitSTunc_IS',   'PitSTunc_OOS',...
        'LeftEntropy_IS', 'LeftEntropy_OOS')
    clear('filename')
end

%% Figure 10. Out-of-sample Predictions.
% (a)/(b) Quantiles
f = figure;
plot(Time, YQ_OOS(:, [jq05, jq50, jq95]), '-b',...
     Time, YQ_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_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% (c)/(d) Downside Entropy
f = figure;
plot(Time, LeftEntropy_OOS, '-b', Time, LeftEntropy_IS, ':k')
datetick('x', 'yyyy')
legend('out of sample', 'in sample')
axis([Time(1), Time(end), -0.5, 1.5])
filename = fullfile(FigSubFolder, ['EntropyOOS_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

%% Figure 11. Out-of-sample Accuracy.
% (a)/(b) Predictive scores
f = figure;
plot(Time, ScoreST_OOS, Time, ScoreSTGDPonly_OOS, '--');
datetick('x', 'yyyy')
axis([Time(1), Time(end), 0.0, 0.4])
ylabel('Scores')
legend('GDP and NFCI', 'GDP only', 'Location', 'NorthWest')
filename = fullfile(FigSubFolder, ['ScoresOOS_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')

% (c)/(d): PITs
% 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;
zST_ecdf = PITtest(PitST_OOS, rvec);
zSTGDPonly_ecdf = PITtest(PitSTGDPonly_OOS, rvec);
if h == 1
    % Use asymptotic 5% critical value from Rossi and Sekhposyan (2017): 1.34
    kappa        = 1.34;
    kappaGDPonly = 1.34;
elseif h == 4
    % Compute bootstrapped 5% critical values
    PITs = [PitST_OOS, PitSTGDPonly_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        = testcritvalues(1, 2, 1);
    kappaGDPonly = testcritvalues(1, 2, 2);
end

% Plot PIT for full quantile regression vs. quantile regression with GDP only
f = figure;
plot(rvec, zST_ecdf, '-b', rvec, zSTGDPonly_ecdf, '--r')
hold on
P = sum(isnan(PitST_OOS) == 0);
plot(rvec, [rvec', rvec' - (kappa / sqrt(P)), rvec' + (kappa / sqrt(P))], '-b');
plot(rvec, [rvec', rvec' - (kappaGDPonly / sqrt(P)), rvec' + (kappaGDPonly / sqrt(P))], '--r');
hold off
xlim([0, 1])
ylim([0, 1])
xlabel('\tau')
ylabel('Empirical CDF')
grid on
legend('GDP and NFCI', 'GDP only', 'Theoretical and 5% Critical Values',...
       'Location', 'SouthEast');
filename = fullfile(FigSubFolder, ['PITS_H', num2str(h), '.pdf']);
printpdf(f, filename);
clear('f', 'filename')
