%% Forecast Unit Root Nonstationary Process Using Monte Carlo Simulations
% Regress the unit root nonstationary, quarterly log GDP onto the CPI using a regression model with ARIMA(1,1,1) errors with known intercept. Forecast log GDP using Monte Carlo simulation.
% Load the US 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);
frstHzn =  T+1:numObs;                % Forecast horizon
hoCPI = DataTable.CPIAUCSL(frstHzn);  % Holdout sample

%
% Fit a regression model with ARIMA(1,1,1). The intercept is not identifiable in a model with integrated errors, so fix its value before estimation.

intercept = 5.867;
Mdl = regARIMA('ARLags',1,'MALags',1,'D',1,'Intercept',intercept);
EstMdl = estimate(Mdl,logGDP,'X',cpi);

%
% Infer unconditional disturbances.

[~,u0] = infer(EstMdl,logGDP,'X',cpi);

%
% Simulate 1000 paths with 15 observations each. Use the inferred unconditional disturbances as presample data.

rng(1); % For reproducibility
GDPF = simulate(EstMdl,15,'NumPaths',1000,...
    'U0',u0,'X',hoCPI);

%
% Plot the simulation mean forecast and approximate 95% forecast intervals.

lower = prctile(GDPF,2.5,2);
upper = prctile(GDPF,97.5,2);
mn = mean(GDPF,2);

figure
plot(dates(end-65:end),log(DataTable.GDP(end-65:end)),'Color',...
   [.7,.7,.7])
datetick
hold on
h1 = plot(dates(frstHzn),lower,'r:','LineWidth',2);
plot(dates(frstHzn),upper,'r:','LineWidth',2)
h2 = plot(dates(frstHzn),mn,'k','LineWidth',2);
h = gca;
ph = patch([repmat(dates(frstHzn(1)),1,2) repmat(dates(frstHzn(end)),1,2)],...
    [h.YLim fliplr(h.YLim)],...
    [0 0 0 0],'b');
ph.FaceAlpha = 0.1;
legend([h1 h2],{'95% Interval','Simulation Mean'},'Location','NorthWest',...
    'AutoUpdate','off')
axis tight
title('{\bf log GDP Forecast - 15 Quarter Horizon}')
hold off

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