%% Estimate Forecast Intervals
% Analyze forecast accuracy using forecast intervals over a three-year horizon. This example follows from Forecast Unconditional Response Series from VEC Model.
% Load the Data_USEconVECModel data set and preprocess the data.

load Data_USEconVECModel
FRED.GDP = 100*log(FRED.GDP);      
FRED.GDPDEF = 100*log(FRED.GDPDEF);
FRED.COE = 100*log(FRED.COE);       
FRED.HOANBS = 100*log(FRED.HOANBS); 
FRED.PCEC = 100*log(FRED.PCEC);     
FRED.GPDI = 100*log(FRED.GPDI);

%
% Estimate a VEC(1) model. Reserve the last three years of data to assess forecast accuracy. Assume that the appropriate cointegration rank is 4, and the H1 Johansen form is appropriate for the model.

bfh = FRED.Time(end) - years(3);
estIdx = FRED.Time < bfh;
Mdl = vecm(7,4,1);
Mdl.SeriesNames = FRED.Properties.VariableNames;
EstMdl = estimate(Mdl,FRED{estIdx,:});

%
% Forecast responses from the estimated model over a three-year horizon. Specify all in-sample observations as a presample. Return the MSE of the forecasts.

numperiods = 12;
Y0 = FRED{estIdx,:};
[Y,YMSE] = forecast(EstMdl,numperiods,Y0);

% Y is a 12-by-7 matrix of forecasted responses. YMSE is a 12-by-1 cell vector of 7-by-7 matrices corresponding to the MSEs.
%
% Extract the main diagonal elements from the matrices in each cell of YMSE. Apply the square root of the result to obtain standard errors.

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

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

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

YFI(:,:,1) = Y - 2*SE;
YFI(:,:,2) = Y + 2*SE;

% Plot the forecasted responses and the last 40 true responses.
%

figure;
subplot(2,2,1)
h1 = plot(FRED.Time((end-39):end),FRED.GDP((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,1));
h3 = plot(FRED.Time(~estIdx),YFI(:,1,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,1,2),'k--');
title('Gross Domestic Product');
ylabel('Index (scaled)');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off
subplot(2,2,2)
h1 = plot(FRED.Time((end-39):end),FRED.GDPDEF((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,2));
h3 = plot(FRED.Time(~estIdx),YFI(:,2,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,2,2),'k--');
title('GDP Deflator');
ylabel('Index (scaled)');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off
subplot(2,2,3)
h1 = plot(FRED.Time((end-39):end),FRED.COE((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,3));
h3 = plot(FRED.Time(~estIdx),YFI(:,3,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,3,2),'k--');
title('Paid Compensation of Employees');
ylabel('Billions of $ (scaled)');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off
subplot(2,2,4)
h1 = plot(FRED.Time((end-39):end),FRED.HOANBS((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,4));
h3 = plot(FRED.Time(~estIdx),YFI(:,4,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,4,2),'k--');
title('Nonfarm Business Sector Hours');
ylabel('Index (scaled)');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off

figure;
subplot(2,2,1)
h1 = plot(FRED.Time((end-39):end),FRED.FEDFUNDS((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,5));
h3 = plot(FRED.Time(~estIdx),YFI(:,5,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,5,2),'k--');
title('Federal Funds Rate');
ylabel('Percent');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off
subplot(2,2,2)
h1 = plot(FRED.Time((end-39):end),FRED.PCEC((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,6));
h3 = plot(FRED.Time(~estIdx),YFI(:,6,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,6,2),'k--');
title('Consumption Expenditures');
ylabel('Billions of $ (scaled)');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off
subplot(2,2,3)
h1 = plot(FRED.Time((end-39):end),FRED.GPDI((end-39):end));
hold on
h2 = plot(FRED.Time(~estIdx),Y(:,7));
h3 = plot(FRED.Time(~estIdx),YFI(:,7,1),'k--');
plot(FRED.Time(~estIdx),YFI(:,7,2),'k--');
title('Gross Private Domestic Investment');
ylabel('Billions of $ (scaled)');
xlabel('Date');
h = gca;
fill([bfh h.XLim([2 2]) bfh],h.YLim([1 1 2 2]),'k',...
    'FaceAlpha',0.1,'EdgeColor','none');
legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
    'Location','best');
hold off

% Copyright 2012 The MathWorks, Inc.