%% Forecast a Regression Model with Multiplicative Seasonal ARIMA Errors
% This example shows how to forecast a multiplicative seasonal ARIMA model using forecast. The response series is monthly international airline passenger numbers from 1949 to 1960.
% Load the airline and recessions data sets. Transform the response.

load('Data_Airline.mat')
load Data_Recessions
y = log(Data);

%
% Construct the predictor (X), which determines whether the country was in a recession during the sampled period. A 0 in row t means the country was not in a recession in month t, and a 1 in row t means that it was in a recession in month t.

X = zeros(numel(dates),1); % Preallocation
for j = 1:size(Recessions,1)
    X(dates >= Recessions(j,1) & dates <= Recessions(j,2)) = 1;
end

%
% Define index sets that partition the data into estimation and forecast samples.

nSim = 60; % Forecast period
T = length(y);
estInds = 1:(T-nSim);
foreInds = (T-nSim+1):T;

%
% Estimate the regression model with multiplicative seasonal ARIMA} (0,1,1)\times(0,1,1)_{12} errors:
% y_t = X_t\beta + u_t
% (1 - L)(1 - L^{12})u_t = (1 + BL)(1 + B_{12}L^{12})\varepsilon_t
% Set the regression model intercept to 0 since it is not identifiable in a model with integrated errors.

Mdl = regARIMA('D',1,'Seasonality',12,'MALags',1,'SMALags',12,...
    'Intercept',0);
EstMdl = estimate(Mdl,y(estInds),'X',X(estInds)); 

%
% Use the estimated coefficients of the model (contained in EstMdl), to generate MMSE forecasts and corresponding mean square errors over a 60-month horizon. Use the observed series as presample data. By default, forecast infers presample innovations and unconditional disturbances using the specified model and observations.

[YF,YMSE] = forecast(EstMdl,nSim,'X0',X(estInds),...
    'Y0',y(estInds),'XF',X(foreInds));
ForecastInt = [YF,YF] + 1.96*[-sqrt(YMSE), sqrt(YMSE)];

figure
h1 = plot(dates,y);
title('{\bf Forecasted Monthly Passenger Totals}')
hold on
h2 = plot(dates(foreInds),YF,'Color','r','LineWidth',2);
h3 = plot(dates(foreInds),ForecastInt,'k--','LineWidth',2);
datetick
legend([h1,h2,h3(1)],'Observations','MMSE Forecasts',...
    '95% MMSE Forecast Intervals','Location','NorthWest')
axis tight
hold off

% The regression model with SMA errors seems to forecast the series well, albeit slightly overestimating. Since the error process is nonstationary, the forecast intervals widen as time increases.
%
% Compare the MMSE forecasts to Monte Carlo forecasts by simulating 500 sample paths from EstMdl over the forecast horizon.

[e0,u0] = infer(EstMdl,y(estInds),'X',X(estInds));
rng(5);
numPaths = 500;
ySim = simulate(EstMdl,nSim,'numPaths',numPaths,...
    'E0',e0,'U0',u0,'X',X(foreInds));
meanYSim = mean(ySim,2);
ForecastIntMC = [prctile(ySim,2.5,2),prctile(ySim,97.5,2)];

figure
h1 = plot(dates(foreInds),y(foreInds));
title('{\bf Forecasted Monthly Passenger Totals}')
hold on
h2 = plot(dates(foreInds),YF,'Color',[0.85,0.85,0.85],...
    'LineWidth',4);
h3 = plot(dates(foreInds),ForecastInt,'--','Color',...
    [0.85,0.85,0.85],'LineWidth',4);
h4 = plot(dates(foreInds),meanYSim,'k','LineWidth',2);
h5 = plot(dates(foreInds),ForecastIntMC,'k--','LineWidth',2);
datetick
legend([h1,h2,h3(1),h4,h5(1)],'Observations',...
    'MMSE Forecasts','95% MMSE Forecast Intervals',...
    'Monte Carlo Forecasts','95% Monte Carlo Forecast Intervals',...
    'Location','NorthWest')
axis tight
hold off

% The MMSE forecasts and Monte Carlo mean forecasts are virtually indistinguishable. However, there are slight discrepancies between the theoretical 95% forecast intervals and the simulation-based 95% forecast intervals.
% Copyright 2015 The MathWorks, Inc.