%% Forecast Responses from Posterior Predictive Distribution
% Consider the 3-D VAR(4) model of docid:econ_ug.mw_93ab52c2-232a-4877-8152-03c57b96825f.
%% Load and Preprocess Data
% Load the US macroeconomic data set.  Compute the inflation rate, stabilize the unemployment and federal funds rates, and remove missing values (the data includes leading missing values only).

load Data_USEconModel
seriesnames = ["INFL" "UNRATE" "FEDFUNDS"];
DataTable.INFL = 100*[NaN; price2ret(DataTable.CPIAUCSL)];
DataTable.DUNRATE = [NaN; diff(DataTable.UNRATE)];
DataTable.DFEDFUNDS = [NaN; diff(DataTable.FEDFUNDS)];
seriesnames(2:3) = "D" + seriesnames(2:3);
rmDataTable = rmmissing(DataTable); 

%
%% Create Prior Model
% Create a weak conjugate prior model. Specify the response series names.

numseries = numel(seriesnames);
numlags = 4;

PriorMdl = conjugatebvarm(numseries,numlags,'SeriesNames',seriesnames);
numcoeffseqn = size(PriorMdl.V,1);
PriorMdl.V = 1e4*eye(numcoeffseqn);

%
%% Forecast Responses Using forecast
% Directly forecast two years (eight quarters) of observations from the posterior predictive distribution. forecast estimates the posterior distribution of the parameters, and then forms the posterior predictive distribution.

rng(1); % For reproducibility
numperiods = 8;
YF = forecast(PriorMdl,numperiods,rmDataTable{:,seriesnames});

% YF is an 8-by-3 matrix of forecasted responses.
% Plot the forecasted responses.

fh = rmDataTable.Time(end) + calquarters(1:8);
for j = 1:PriorMdl.NumSeries
    subplot(3,1,j)
    plot(rmDataTable.Time(end - 20:end),rmDataTable{end - 20:end,seriesnames(j)},'r',...
        [rmDataTable.Time(end) fh],[rmDataTable{end,seriesnames(j)}; YF(:,j)],'b');
    legend("Observed","Forecasted",'Location','NorthWest')
    title(seriesnames(j))
end

%
%% Forecast Responses Using simsmooth
% Configure the data set for obtaining forecasts from simsmooth by concatenating a numperiods-by-numseries timetable of missing values to the end of the set.

fTT = array2timetable(NaN(numperiods,numseries),'RowTimes',fh,...
    'VariableNames',seriesnames);
frmDataTable = [rmDataTable(:,seriesnames); fTT];
tail(frmDataTable)

% Forecast the VAR model using simsmooth. Like forecast, simsmooth estimates the posterior distribution, so it requires a prior model and data. Specify the data containing the trailing NaNs.

[~,~,~,YMean] = simsmooth(PriorMdl,frmDataTable.Variables);

% YMean is almost equal to frmDataTable, with these exceptions:
% YMean excludes rows corresponding to the presample period frmDataTable(1:4,:).
% If frmDataTable(j,k) is NaN, YMean(j,k) is the mean of the posterior predictive distribution of response k at time j.
% Create a timetable from YMean.

YMeanTT = array2timetable(YMean,'RowTimes',frmDataTable.Time((PriorMdl.P + 1):end),...
    'VariableNames',seriesnames);

% Plot the forecasted responses.

for j = 1:PriorMdl.NumSeries
    subplot(3,1,j)
    plot(YMeanTT.Time((end - 20 - numperiods):(end - numperiods)),YMeanTT{(end - 20 - numperiods):(end - numperiods),j},'r',...
        YMeanTT.Time((end - numperiods):end),YMeanTT{(end - numperiods):end,j},'b');
        legend("Observed","Forecasted",'Location','NorthWest')
    title(seriesnames(j))
end

% Copyright 2020 The MathWorks, Inc.
%
%