%% Forecast Stationary Process Using Monte Carlo Simulations
% Regress the stationary, quarterly log GDP onto the CPI using a regression model with ARMA(1,1) errors, and forecast log GDP using Monte Carlo simulation.
% Load the US Macroeconomic data set and preprocess the data.

load Data_USEconModel;
logGDP = log(DataTable.GDP);
dlogGDP = diff(logGDP);          % For stationarity
dCPI = diff(DataTable.CPIAUCSL); % For stationarity
numObs = length(dlogGDP);
gdp = dlogGDP(1:end-15);   % Estimation sample
cpi = dCPI(1:end-15);
T = length(gdp);        % Effective sample size
frstHzn =  T+1:numObs;  % Forecast horizon
hoCPI = dCPI(frstHzn);  % Holdout sample
dts = dates(2:end);     % Date nummbers

%
% Fit a regression model with ARMA(1,1) errors.

Mdl = regARIMA('ARLags',1,'MALags',1);
EstMdl = estimate(Mdl,gdp,'X',cpi);

%
% Infer unconditional disturbances.

[~,u0] = infer(EstMdl,gdp,'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(dts(end-65:end),dlogGDP(end-65:end),'Color',[.7,.7,.7])
datetick
hold on
h1 = plot(dts(frstHzn),lower,'r:','LineWidth',2);
plot(dts(frstHzn),upper,'r:','LineWidth',2)
h2 = plot(dts(frstHzn),mn,'k','LineWidth',2);
h = gca;
ph = patch([repmat(dts(frstHzn(1)),1,2) repmat(dts(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

% Copyright 2012 The MathWorks, Inc.