%% Forecast Responses Using Posterior Predictive Distribution
% Consider the multiple linear regression model that predicts the US real gross national product (GNPR) using a linear combination of industrial production index (IPI), total employment (E), and real wages (WR).
% \texttt{GNPR}_t=\beta_0+\beta_1\texttt{IPI}_t+\beta_2\texttt{E}_t+\beta_3\texttt{WR}_t+\varepsilon_t.
% For all t, \varepsilon_t is a series of independent Gaussian disturbances with a mean of 0 and variance \sigma^2.
% Assume these prior distributions:
% \beta\vert\sigma^2 \sim N_4\left({\rm M},\sigma^2 V\right). {\rm M} is a 4-by-1 vector of means, and V is a scaled 4-by-4 positive definite covariance matrix.
% \sigma^2 \sim IG(A,B). A and B are the shape and scale, respectively, of an inverse gamma distribution.
% These assumptions and the data likelihood imply a normal-inverse-gamma conjugate model.
% Create a normal-inverse-gamma conjugate prior model for the linear regression parameters. Specify the number of predictors p and the variable names.

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

% Mdl is a conjugateblm Bayesian linear regression model object representing the prior distribution of the regression coefficients and disturbance variance.
%
% Load the Nelson-Plosser data set. Create variables for the predictor and response data. Hold out the last 10 periods of data from estimation so you can use them to  forecast real GNP.

load Data_NelsonPlosser
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

%
% Estimate the marginal posterior distributions. Turn off the estimation display.

PosteriorMdl = estimate(PriorMdl,X,y,'Display',false);

% PosteriorMdl is a conjugateblm model object that contains the posterior distributions of \beta and \sigma^2.
%
% Forecast responses by using the posterior predictive distribution and the future predictor data XF. Plot the true values of the response and the forecasted values.

yF = forecast(PosteriorMdl,XF);

figure;
plot(dates,DataTable.GNPR);
hold on
plot(dates((end - fhs + 1):end),yF)
h = gca;
p = 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(p,'bottom');
legend('Forecast Horizon','True GNPR','Forecasted GNPR','Location','NW')
title('Real Gross National Product: 1909 - 1970');
ylabel('rGNP');
xlabel('Year');
hold off

% yF is a 10-by-1 vector of future values of real GNP corresponding to the future predictor data.
%
% Estimate the forecast root mean squared error (RMSE).

frmse = sqrt(mean((yF - yFT).^2))

% The forecast RMSE is a relative measure of forecast accuracy. Specifically, you estimate several models using different assumptions. The model with the lowest forecast RMSE is the best-performing model of the ones being compared.
% Copyright 2012 The MathWorks, Inc.