%% Simulate a Regression Model with Stationary Multiplicative Seasonal Errors
% This example shows how to simulate sample paths from a regression model with multiplicative seasonal ARIMA errors using simulate. The time series is monthly international airline passenger numbers from 1949 to 1960.
% Load the airline and recessions data sets.

load('Data_Airline.mat')
load Data_Recessions

%
% Transform the airline data by applying the logarithm, and the 1st and 12th differences.

y = Data;
logY = log(y);
DiffPoly = LagOp([1 -1]);
SDiffPoly = LagOp([1 -1],'Lags',[0, 12]);
dLogY = filter(DiffPoly*SDiffPoly,logY);

%
% 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
X = X(14:end); % Remove the first 14 observations for consistency
dates = dates(14:end);

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

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

%
% Estimate the regression model with multiplicative seasonal errors:
% y_t = c + X_t\beta + u_t
% u_t = (1 + BL)(1 + B_{12}L^{12})\varepsilon_t.

Mdl = regARIMA('MALags',1,'SMALags',12);
EstMdl = estimate(Mdl,dLogY(estInds),'X',X(estInds));

%
% Use the estimated coefficients of the model (contained in EstMdl) to simulate 25 realizations of airline passenger counts over the 60-month horizon. Infer the residuals, and use them as a presample.

[~,u0] = infer(EstMdl,dLogY(estInds),'X',X(estInds));
rng(5);
numPaths = 25;
dLogYSim = simulate(EstMdl,60,'numPaths',numPaths,'U0',u0,'X',X(foreInds));
meanDLogYSim = mean(dLogYSim,2);

figure
h1 = plot(dates(estInds),dLogY(estInds));
title('{\bf Transformed, Simulated Monthly Passenger Totals}')
hold on
plot(dates(foreInds),dLogYSim,'Color',[.85,.85,.85])
h2 = plot(dates(foreInds),meanDLogYSim,'k.-','LineWidth',2);
plot([dates(estInds(end)),dates(foreInds(1))],...
    [repmat(dLogY(estInds(end)),numPaths,1),dLogYSim(1,:)'],...
    'Color',[.85,.85,.85])
plot([dates(estInds(end)),dates(foreInds(1))],...
    [dLogY(estInds(end)),meanDLogYSim(1)],'k.-','LineWidth',2)
plot(dates(foreInds),dLogY(foreInds))
datetick
legend([h1,h2],'Observations','Simulation Mean','Location','NorthWest')
axis tight
hold off

% The regression model with SMA errors seems to forecast the series well.
% Check the predictive performance of the model by:
% Varying the size of the forecast period
% Estimating the prediction mean square error (PMSE)
% Choosing the model with the lowest PMSE
% Copyright 2015 The MathWorks, Inc.