%% Work with Prior and Posterior Distributions
% Consider the 3-D VAR(4) model of docid:econ_ug#mw_22b40c74-722c-4ed0-bd9d-6de2110b22c7. Estimate the posterior distribution, and generate forecasts from the corresponding posterior predictive distribution.
%
%% Load and Preprocess Data
% Load the US macroeconomic data set.  Compute the inflation rate.  Plot all response series.

load Data_USEconModel
seriesnames = ["INFL" "UNRATE" "FEDFUNDS"];
DataTable.INFL = 100*[NaN; price2ret(DataTable.CPIAUCSL)];

figure
plot(DataTable.Time,DataTable{:,seriesnames})
legend(seriesnames)

% Stabilize the unemployment and federal funds rates by applying the first difference to each series.

DataTable.DUNRATE = [NaN; diff(DataTable.UNRATE)];
DataTable.DFEDFUNDS = [NaN; diff(DataTable.FEDFUNDS)];
seriesnames(2:3) = "D" + seriesnames(2:3);

% Remove all missing values from the data.

rmDataTable = rmmissing(DataTable);

%
%% Create Prior Model
% Create a conjugate Bayesian VAR(4) prior model for the three response series. Specify the response variable names.

numseries = numel(seriesnames);
numlags = 4;

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

%
%% Estimate Posterior Distribution
% Estimate the posterior distribution by passing the prior model and entire data series to estimate.

PosteriorMdl = estimate(PriorMdl,rmDataTable{:,seriesnames},'Display','equation');

% Because the prior is conjugate for the data likelihood, the posterior is a conjugatebvarm object. By default, estimate uses the first four observations as a presample to initialize the model.
%
%% Generate Forecasts from Posterior Predictive Distribution
% From the posterior predictive distribution, generate forecasts over a two-year horizon.  Because sampling from the posterior predictive distribution requires the entire data set, specify the prior model in forecast instead of the posterior.

fh = 8;
FY = forecast(PriorMdl,fh,rmDataTable{:,seriesnames});

% FY is an 8-by-3 matrix of forecasts.
% Plot the end of the data set and the forecasts.

fp = rmDataTable.Time(end) + calquarters(1:fh);
figure
plotdata = [rmDataTable{end - 10:end,seriesnames}; FY];
plot([rmDataTable.Time(end - 10:end); fp'],plotdata)
hold on
plot([fp(1) fp(1)],ylim,'k-.')
legend(seriesnames)
title('Data and Forecasts')
hold off

%
%% Compute Impulse Responses
% Plot impulse response functions by passing posterior estimations to armairf.

armairf(PosteriorMdl.AR,[],'InnovCov',PosteriorMdl.Covariance)

% Copyright 2020 The MathWorks, Inc.