%% Forecast Responses Using Posterior Predictive Distribution
% Consider the linear regression model in Estimate Marginal Posterior Distributions.
% Create a prior model for the regression coefficients and disturbance variance, then estimate the marginal posterior distributions. Hold out the last 10 periods of data from estimation so you can use them to forecast real GNP. Turn the estimation display off.

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

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',VarNames);

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws,'Display',false);

%
% Forecast responses 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;
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','Location','NW')
title('Real Gross National Product');
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.