%% Forecast VARX Model
% Consider the 2-D VARX(1) model for the US real GDP (RGDP) and investment (GCE) rates that treats the personal consumption (PCEC) rate  as exogenous:
% \left[ {\begin{array}{*{20}{l}}
% {\texttt{RGDP}_t}\\
% {\texttt{GCE}_t}
% \end{array}} \right] = c + \Phi\left[ {\begin{array}{*{20}{l}}
% {\texttt{RGDP}_{t-1}\\
% {\texttt{GCE}_{t-1}
% \end{array}} \right] + \texttt{PCEC}_t\beta+\varepsilon_t.
% For all t, \varepsilon_t is a series of independent 2-D normal innovations with a mean of 0 and covariance \Sigma. Assume the following prior distributions:
% {\left\lbrack \begin{array}{ccc}
% \Phi  & c & \beta
% \end{array}\right\rbrack }^{\prime } \left|\Sigma \sim N_{4\times 2} \left(M,V,\Sigma \right)\right., where M is a 4-by-2 matrix of means and V is the 4-by-4 among-coefficient scale matrix. Equivalently, \textrm{vec}\left({\left\lbrack \begin{array}{ccc}
% \Phi  & c & \beta
% \end{array}\right\rbrack }^{\prime } \right)\left|\Sigma \sim N_8 \left(\textrm{vec}\left(M\right),\Sigma \otimes \;V\right)\right..
% \Sigma\sim{\rm Inverse\,Wishart}\left(\Omega,\nu\right),where Ω is the 2-by-2 scale matrix and \nu is the degrees of freedom.
% Load the US macroeconomic data set.  Compute the real GDP, investment, and personal consumption rate series.  Remove all missing values from the resulting series.

load Data_USEconModel
DataTable.RGDP = DataTable.GDP./DataTable.GDPDEF;
seriesnames = ["PCEC"; "RGDP"; "GCE"];
rates = varfun(@price2ret,DataTable,'InputVariables',seriesnames);
rates = rmmissing(rates);
rates.Properties.VariableNames = seriesnames;

% Create a conjugate prior model for the 2-D VARX(1) model parameters.

numseries = 2;
numlags = 1;
numpredictors = 1;
PriorMdl = conjugatebvarm(numseries,numlags,'NumPredictors',numpredictors,...
    'SeriesNames',seriesnames(2:end));

% Create index sets that partition the data into estimation and forecast samples. Specify a forecast horizon of two years.

T = size(rates,1);
numperiods = 8;
idxest = 1:(T - numperiods); % Includes presample
idxf = (T - numperiods + 1):T;
idxtot = [idxest idxf];

% The simulation smoother forecasts by imputing missing values. Therefore, create a data set that contains missing values for the responses in the forecast horizon.

missingrates = rates;
missingrates{idxf,PriorMdl.SeriesNames} = nan(numperiods,PriorMdl.NumSeries);

% Forecast responses in the forecast horizon. Specify presample observations and the exogenous predictor data. Return standard deviations of the posterior predictive distribution.

rng(1) % For reproducibility
[~,~,~,YMean,YStd] = simsmooth(PriorMdl,missingrates{:,PriorMdl.SeriesNames},...
    'X',missingrates{:,seriesnames(1)});

% Create a timetable from YMean.

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

% Plot the forecasted responses.

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

% Copyright 2020 The MathWorks, Inc.