%% Forecast Unconditional Response Series from VAR(4) Model
% Fit a VAR(4) model to the consumer price index (CPI) and unemployment rate data. Then, forecast unconditional MMSE responses from the estimated model.
% Load the Data_USEconModel data set.

load Data_USEconModel

%
% Plot the two series on separate plots.

figure;
plot(DataTable.Time,DataTable.CPIAUCSL);
title('Consumer Price Index');
ylabel('Index');
xlabel('Date');
figure;
plot(DataTable.Time,DataTable.UNRATE);
title('Unemployment Rate');
ylabel('Percent');
xlabel('Date');

%
% Stabilize the CPI by converting it to a series of growth rates. Synchronize the two series by removing the first observation from the unemployment rate series.

rcpi = price2ret(DataTable.CPIAUCSL);
unrate = DataTable.UNRATE(2:end);
Data = array2timetable([rcpi unrate],'RowTimes',DataTable.Time(2:end),...
    'VariableNames',{'rcpi','unrate'});

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

Mdl = varm(2,4);

%
% Estimate the model using the entire data set.

EstMdl = estimate(Mdl,Data.Variables);

% EstMdl is a fully specified, estimated varm model object.
%
% Forecast responses from the estimated model over a three-year horizon. Specify the entire data set as presample observations.

numperiods = 12;
Y0 = Data.Variables;
Y = forecast(EstMdl,numperiods,Y0);

% Y is a 12-by-2 matrix of forecasted responses. The first and second columns contain the simulated CPI growth rate and unemployment rate, respectively.
%
% Plot the forecasted responses and the last 50 true responses.

fh = dateshift(Data.Time(end),'end','quarter',1:12);

figure;
h1 = plot(Data.Time((end-49):end),Data.rcpi((end-49):end));
hold on;
h2 = plot(fh,Y(:,1));
title('CPI Growth Rate');
ylabel('Growth rate');
xlabel('Date');
h = gca;
fill([Data.Time(end) fh([end end]) Data.Time(end)],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2],'True','Forecast')
hold off;

figure;
h1 = plot(Data.Time((end-49):end),Data.unrate((end-49):end));
hold on;
h2 = plot(fh,Y(:,2));
title('Unemployment Rate');
ylabel('Percent');
xlabel('Date');
h = gca;
fill([Data.Time(end) fh([end end]) Data.Time(end)],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2],'True','Forecast','Location','northwest')
hold off;

% Copyright 2012 The MathWorks, Inc.