%% Forecast VARX(4) Model
% Estimate a four-degree vector autoregression model including exogenous predictors (VARX(4)) of the consumer price index (CPI), the unemployment rate, and the gross domestic product (GDP). Include a linear regression component containing the current quarter and the last four quarters of government consumption expenditures and investment (GCE). Forecast a response path from the estimated model.
% Load the Data_USEconModel data set. Compute the real GDP.

load Data_USEconModel
DataTable.RGDP = DataTable.GDP./DataTable.GDPDEF*100;

%
% Plot all variables on separate plots.

figure;
subplot(2,2,1)
plot(DataTable.Time,DataTable.CPIAUCSL);
ylabel('Index');
title('Consumer Price Index');
subplot(2,2,2)
plot(DataTable.Time,DataTable.UNRATE);
ylabel('Percent');
title('Unemployment Rate');
subplot(2,2,3)
plot(DataTable.Time,DataTable.RGDP);
ylabel('Output');
title('Real Gross Domestic Product');
subplot(2,2,4)
plot(DataTable.Time,DataTable.GCE);
ylabel('Billions of $');
title('Government Expenditures');

%
% Stabilize the CPI, GDP, and GCE by converting each to a series of growth rates. Synchronize the unemployment rate series with the others by removing its first observation.

inputVariables = {'CPIAUCSL' 'RGDP' 'GCE'};
Data = varfun(@price2ret,DataTable,'InputVariables',inputVariables);
Data.Properties.VariableNames = inputVariables;
Data.UNRATE = DataTable.UNRATE(2:end);

%
% Expand the GCE rate series to a matrix that includes its current value and up through four lagged values. Remove the GCE variable from Data.

rgcelag4 = lagmatrix(Data.GCE,0:4);
Data.GCE = [];

%
% Create a default VAR(4) model using the shorthand syntax.

Mdl = varm(3,4);
Mdl.SeriesNames = ["rcpi" "unrate" "rgdpg"];

%
% Estimate the model using all but the last three years of data. Specify the GCE matrix as data for the regression component.

bfh = Data.Time(end) - years(3);
estIdx = Data.Time < bfh;
EstMdl = estimate(Mdl,Data(estIdx,:).Variables,'X',rgcelag4(estIdx,:));

%
% Forecast a path of quarterly responses three years into the future.

Y0 = Data(estIdx,:).Variables;
Y = forecast(EstMdl,12,Data(estIdx,:).Variables,'X',rgcelag4(~estIdx,:));

% Y is a 12-by-3 matrix of simulated responses. The columns correspond to the CPI growth rate, unemployment rate, and real GDP growth rate, respectively.
%
% Plot the forecasted responses and the last 50 true responses.

figure;
for j = 1:Mdl.NumSeries
    subplot(2,2,j)
    h1 = plot(Data.Time((end-49):end),Data{(end-49):end,j});
    hold on;
    h2 = plot(Data.Time(~estIdx),Y(:,j));
    title(Mdl.SeriesNames{j});
    h = gca;
    fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
        'FaceAlpha',0.1,'EdgeColor','none');
    hold off;
end

hl = legend([h1 h2],'Data','Forecast');
hl.Position = [0.6 0.25 hl.Position(3:4)];

% Copyright 2012 The MathWorks, Inc.