%% Forecast Using a Regression Model with ARIMA Errors and a Known Intercept
% Forecast unit root nonstationary, log GDP using a regression model with ARIMA(1,1,1) errors, including CPI as a predictor and a known intercept.
% Load the U.S. Macroeconomic data set and preprocess the data.

load Data_USEconModel;
numObs = length(DataTable.GDP);
logGDP = log(DataTable.GDP(1:end-15));
cpi = DataTable.CPIAUCSL(1:end-15);
T = length(logGDP);                  % Effective sample size
frstHzn =  T+1:numObs;               % Forecast horizon
hoCPI = DataTable.CPIAUCSL(frstHzn); % Holdout sample

%
% Specify the model for the estimation period.

Mdl = regARIMA('ARLags',1,'MALags',1,'D',1);

%
% The intercept is not identifiable in a model with integrated errors, so fix its value before estimation. One way to do this is to estimate the intercept using simple linear regression.

Reg4Int = [ones(T,1), cpi]\logGDP;
intercept = Reg4Int(1);

% Consider performing a sensitivity analysis by using a grid of intercepts.
%
% Set the intercept and fit the regression model with ARIMA(1,1,1) errors.

Mdl.Intercept = intercept;
EstMdl = estimate(Mdl,logGDP,'X',cpi,...
    'Display','off')

%
% Forecast GDP over a 15-quarter horizon. Use the estimation sample as a presample for the forecast.

[gdpF,gdpMSE] = forecast(EstMdl,15,'Y0',logGDP,...
    'X0',cpi,'XF',hoCPI);

%
% Plot the forecasts and 95% forecast intervals.

figure
h1 = plot(dates(end-65:end),log(DataTable.GDP(end-65:end)),...
    'Color',[.7,.7,.7]);
datetick
hold on
h2 = plot(dates(frstHzn),gdpF,'b','LineWidth',2);
h3 = plot(dates(frstHzn),gdpF+1.96*sqrt(gdpMSE),'r:',...
		'LineWidth',2);
plot(dates(frstHzn),gdpF-1.96*sqrt(gdpMSE),'r:',...
    'LineWidth',2);
ha = gca;

title(['{\bf Forecasts and Approximate 95% }'...
			'{\bf Forecast Intervals for log GDP}']);
ph = patch([repmat(dates(frstHzn(1)),1,2) repmat(dates(frstHzn(end)),1,2)],...
        [ha.YLim fliplr(ha.YLim)],...
        [0 0 0 0],'b');
ph.FaceAlpha = 0.1;
legend([h1 h2 h3],{'Observed GDP','Forecasted GDP',...
		'95% Forecast Interval'},'Location','Best','AutoUpdate','off');
axis tight
hold off

% The unconditional disturbances, u_{t}, are nonstationary, therefore the widths of the forecast intervals grow with time.
% Copyright 2012 The MathWorks, Inc.