%% Forecast VECX Model
% Consider the model and data in 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);

%
% The Data_Recessions data set contains the beginning and ending serial dates of recessions. Load the data set. Convert the matrix of date serial numbers to a datetime array.

load Data_Recessions
dtrec = datetime(Recessions,'ConvertFrom','datenum');

%
% Create a dummy variable that identifies periods in which the U.S. was in a recession or worse. Specifically, the variable should be 1 if FRED.Time occurs during a recession, and 0 otherwise.

isin = @(x)(any(dtrec(:,1) <= x & x <= dtrec(:,2)));
isrecession = double(arrayfun(isin,FRED.Time));

%
% Create a VEC(1) model using the shorthand syntax. Assume that the appropriate cointegration rank is 4. You do not have to specify the presence of a regression component when creating the model. Specify the variable names.

Mdl = vecm(7,4,1);
Mdl.SeriesNames = FRED.Properties.VariableNames;

%
% Estimate the model using all but the last three years of data. Specify the predictor identifying whether the observation was measured during a recession.

bfh = FRED.Time(end) - years(3);
estIdx = FRED.Time < bfh;
EstMdl = estimate(Mdl,FRED{estIdx,:},'X',isrecession(estIdx));

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

Y0 = FRED{estIdx,:};
Y = forecast(EstMdl,12,Y0,'X',isrecession(~estIdx));

% Y is a 12-by-7 matrix of simulated responses. Rows correspond to the forecast horizon, and columns correspond to the variables in EstMdl.SeriesNames.
%
% 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));
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],'True','Forecast','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));
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],'True','Forecast','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));
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],'True','Forecast','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));
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],'True','Forecast','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));
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],'True','Forecast','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));
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],'True','Forecast','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));
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],'True','Forecast','Location','Best')
hold off

% Copyright 2012 The MathWorks, Inc.