%% Estimate Forecast Intervals
% Consider the regression model in Forecast Responses Using Posterior Predictive Distribution.
% Create a normal-inverse-gamma semiconjugate prior model for the linear regression parameters. Specify the number of predictors p and the names of the regression coefficients.

p = 3;
PriorMdl = bayeslm(p,'ModelType','semiconjugate','VarNames',["IPI" "E" "WR"]);

%
% Load the Nelson-Plosser data set. Create variables for the response and predictor series.

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

%
% Hold out the last 10 periods of data from estimation so you can use them to forecast real GNP. Turn off the estimation display.

fhs = 10; % Forecast horizon size
X = DataTable{1:(end - fhs),PriorMdl.VarNames(2:end)}; 
y = DataTable{1:(end - fhs),'GNPR'};
XF = DataTable{(end - fhs + 1):end,PriorMdl.VarNames(2:end)}; % Future predictor data
yFT = DataTable{(end - fhs + 1):end,'GNPR'};                  % True future responses

%
% Forecast responses and their covariance matrix by using the posterior predictive distribution and the future predictor data XF. Specify the in-sample observations X and y (the observations from which MATLAB® composes the posterior).

[yF,YFCov] = forecast(PriorMdl,XF,X,y);

%
% Because the predictive posterior distribution is not analytical, a reasonable approximation to a set of 95% credible intervals is
% \hat y_i \pm z_{0.975}SE(\hat y_i),
% for all i in the forecast horizon. Estimate 95% credible intervals for the forecasts using this formula.

n = sum(all(~isnan([X y]')));
cil = yF - norminv(0.975)*sqrt(diag(YFCov));
ciu = yF + norminv(0.975)*sqrt(diag(YFCov));

%
% Plot the data, forecasts, and forecast intervals.

figure;
plot(dates(end-30:end),DataTable.GNPR(end-30:end));
hold on
h = gca;
plot(dates((end - fhs + 1):end),yF)
plot(dates((end - fhs + 1):end),[cil ciu],'k--')
hp = patch([dates(end - fhs + 1) dates(end) dates(end) dates(end - fhs + 1)],...
    h.YLim([1,1,2,2]),[0.8 0.8 0.8]);
uistack(hp,'bottom');
legend('Forecast horizon','True GNPR','Forecasted GNPR',...
    'Credible interval','Location','NW')
title('Real Gross National Product: 1909 - 1970');
ylabel('rGNP');
xlabel('Year');
hold off

% Copyright 2012 The MathWorks, Inc.