www.gusucode.com > econ 案例源码程序 matlab代码 > econ/ForecastRegressionModelWithMultiplicativeSeasonalARIMAErExample.m

    %% 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.

% Copyright 2015 The MathWorks, Inc.

load(fullfile(matlabroot,'examples','econ','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 forcast 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.