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.