%% Forecast VAR Model
% This example shows how to use forecast to forecast a VAR model.
% forecast enables you to generate MMSE forecasts with error estimates. forecast requires:
% A fully-specified model (for example , EstMdl in what follows)
% The number of periods for the forecast (for example, numperiods in what follows)
% A presample time series (e.g., Y(end-3:end,:) in what follows)
% forecast optionally takes:
% An exogenous data series
% Extra paths
% Load the Data_USEconModel data set. This example uses two time series: the logarithm of real GDP, and the real 3-month T-bill rate, both differenced to be approximately stationary. Suppose that a VAR(4) model is appropriate to describe the time series.

load Data_USEconModel
DEF = log(DataTable.CPIAUCSL);
GDP = log(DataTable.GDP);
rGDP = diff(GDP - DEF); % Real GDP is GDP - deflation
TB3 = 0.01*DataTable.TB3MS;
dDEF = 4*diff(DEF); % Scaling
rTB3 = TB3(2:end) - dDEF; % Real interest is deflated
Y = [rGDP,rTB3];

%
% Fit a VAR(4) model specification:

Mdl = varm(2,4);
Mdl.SeriesNames = {'Transformed real GDP','Transformed real 3-mo T-bill rate'};
EstMdl = estimate(Mdl,Y);

%
% Predict the evolution of the time series:

numperiods = 21;
FDates = dateshift(DataTable.Time(end),'end','quarter',1:numperiods);
[Forecast,ForecastMSE] = forecast(EstMdl,numperiods,Y(end-3:end,:));

%
% Extract the main diagonal elements from the matrices in each cell of ForecastMSE. Apply the square root of the result to obtain standard errors.

extractMSE = @(x)diag(x)';
MSE = cellfun(extractMSE,ForecastMSE,'UniformOutput',false);
SE = sqrt(cell2mat(MSE));

%
% Estimate approximate 95% forecast intervals for each response series.

ForecastFI = zeros(numperiods,Mdl.NumSeries,2);

ForecastFI(:,:,1) = Forecast - 2*SE;
ForecastFI(:,:,2) = Forecast + 2*SE;

%
% Plot the forecasts and the data.

figure;
h1 = plot(DataTable.Time((end-49):end),Y((end-49):end,1));
hold on;
h2 = plot(FDates,Forecast(:,1));
h3 = plot(FDates,ForecastFI(:,1,1),'k--');
plot(FDates,ForecastFI(:,1,2),'k--');
title('Real GDP Growth Rate');
h = gca;
fill([FDates(1) h.XLim([2 2]) FDates(1)],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','northwest')
hold off;

figure;
h1 = plot(DataTable.Time((end-49):end),Y((end-49):end,2));
hold on;
h2 = plot(FDates,Forecast(:,2));
h3 = plot(FDates,ForecastFI(:,2,1),'k--');
plot(FDates,ForecastFI(:,2,2),'k--');
title('3-month T-bill rate');
h = gca;
fill([FDates(1) h.XLim([2 2]) FDates(1)],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','northwest')
hold off;

% Copyright 2012 The MathWorks, Inc.