%% Estimate Standard Deviations of Posterior Predictive Distribution
% Consider the 3-D VAR(4) model of docid:econ_ug.mw_06d1f680-c634-4382-9b19-b10ba3b3f1f0.
% Load the US macroeconomic data set.  Compute the inflation rate, stabilize the unemployment and federal funds rates, and remove missing values.

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 a diffuse prior model. Specify the response series names.

numseries = numel(seriesnames);
numlags = 4;

PriorMdl = bayesvarm(numseries,numlags,'SeriesNames',seriesnames);

% Directly forecast two years (eight quarters) of response observations from the posterior predictive distribution. Return the posterior standard deviations.

numperiods = 8;
[YF,YStd] = forecast(PriorMdl,numperiods,rmDataTable{:,seriesnames});

% YF and YStd are 8-by-3 matrices of forecasted responses and corresponding standard deviations, respectively.
% Plot the forecasted responses and approximate 95% credible intervals.

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

% Copyright 2020 The MathWorks, Inc.
%
%