%% Model Seasonal Lag Effects Using Indicator Variables
% This example shows how to estimate a seasonal ARIMA model:
% Model the seasonal effects using a multiplicative seasonal model.
% Use indicator variables as a regression component for the seasonal effects, called seasonal dummies.
% Subsequently, their forecasts show that the methods produce similar results. The time series is monthly international airline passenger numbers from 1949 to 1960.
%% Step 1. Load the data.
% Load the data set Data_Airline, and plot the natural log of the monthly passenger totals counts.

load('Data_Airline.mat')
dat = log(Data); % Transform to logarithmic scale 
T = size(dat,1);
y = dat(1:103);  % Estimation sample

% y is the part of dat used for estimation, and the rest of dat is the holdout sample to compare the two models' forecasts.
%
%% Step 2. Define and fit the model specifying seasonal lags.
% Create an ARIMA(0,1,1)(0,1,1)_{12} model
% (1 - L)(1 - {L^{12}}){y_t} = (1 + {\theta _1}L)(1 + {\Theta _{12}}{L^{12}}){\varepsilon _t},
% where \varepsilon_t is an independent and identically distributed normally distributed series with mean 0 and variance \sigma^2. Use estimate to fit Mdl1 to y.

Mdl1 = arima('Constant',0,'MALags',1,'D',1,...
    'SMALags',12,'Seasonality',12);
EstMdl1 = estimate(Mdl1,y);

% The fitted model is
% (1 - L)(1 - {L^{12}}){y_t} = (1 - 0.357L)(1 - 0.615{L^{12}}){\varepsilon _t},
% where \varepsilon_t is an iid normally distributed series with mean 0 and variance 0.0013.
%
%% Step 3. Define and fit the model using seasonal dummies.
% Create an ARIMAX(0,1,1) model with period 12 seasonal differencing and a regression component,
% (1 - L)(1 - {L^{12}}){y_t} = (1 - 0.357L)(1 - 0.615{L^{12}}){\varepsilon _t},
% \{x_t; t = 1,...,T\} is a series of T column vectors having length 12 that indicate in which month observation t was measured. A 1 in row i of x_t indicates that the observation was measured in month i, the rest of the elements are 0s.
% Note that if you include an additive constant in the model, then the T rows of the design matrix X are composed of the row vectors  \left[ {\begin{array}{*{20}{c}}1&{x_t^\prime }\end{array}} \right]. Therefore, X is rank deficient, and one regression coefficient is not identifiable. A constant is left out of this example to avoid distraction from the main purpose. Format the in-sample X matrix

X  = dummyvar(repmat((1:12)',12,1));

% Format the presample X matrix
X0 = [zeros(1,11) 1; dummyvar((1:12)')];
Mdl2 = arima('Constant',0,'MALags',1,'D',1,...
    'Seasonality',12);
EstMdl2   = estimate(Mdl2,y,'X',[X0; X]);

% The fitted model is
% (1 - L)(1 - {L^{12}}){y_t} = X_t \hat \beta  + (1 - 0.407L){\varepsilon _t},
% where \varepsilon_t is an iid normally distributed series with mean 0 and variance 0.0017 and \hat\beta is a column vector with the values Beta1 - Beta12. Note that the estimates MA{1} and Variance between Mdl1 and Mdl2 are not equal.
%
%% Step 4. Forecast using both models.
% Use forecast to forecast both models 41 periods into the future from July 1957. Plot the holdout sample using these forecasts.

yF1 = forecast(EstMdl1,41,y);
yF2 = forecast(EstMdl2,41,y,'X0',X(1:103,:),'XF',X(104:end,:));
l1 = plot(100:T,dat(100:end),'k','LineWidth',3);
hold on
l2 = plot(104:144,yF1,'-r','LineWidth',2);
l3 = plot(104:144,yF2,'-b','LineWidth',2);
hold off
title('Passenger Data: Actual vs. Forecasts')
xlabel('Month')
ylabel('Logarithm of Monthly Passenger Data')
legend({'Observations','Polynomial Forecast',...
 			'Regression Forecast'},'Location','NorthWest')

% Though they overpredict the holdout observations, the forecasts of both models are almost equivalent. One main difference between the models is that EstMdl1 is more parsimonious than EstMdl2.
% References:
% Box, G. E. P., G. M. Jenkins, and G. C. Reinsel. Time Series Analysis: Forecasting and Control. 3rd ed. Englewood Cliffs, NJ: Prentice Hall, 1994.
% Copyright 2015 The MathWorks, Inc.