Симулируйте модели регрессии с мультипликативными сезонными ошибками

Симулируйте модель регрессии со стационарными мультипликативными сезонными ошибками

В этом примере показано, как симулировать демонстрационные пути из модели регрессии с мультипликативными сезонными ошибками ARIMA simulate. Временные ряды являются ежемесячными международными числами авиапассажира от 1 949 до 1960.

Загрузите наборы данных рецессий и авиакомпания.

load('Data_Airline.mat')
load Data_Recessions

Преобразуйте данные авиакомпании путем применения логарифма и 1-х и 12-х различий.

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

Создайте предиктор (X), который определяет, была ли страна в рецессии в произведенный период. 0 последовательно t означает, что страна не была в рецессии в месяце t, и 1 последовательно t означает, что это было в рецессии в месяце 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);

Задайте наборы индекса тот раздел данные в выборки прогноза и оценку.

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

Оцените модель регрессии с мультипликативными сезонными ошибками:

$y_{t} = c + X_{t} β + u_{t}$

$u_{t} = (1 + B L) (1 + B_{12} L^{12}) ε_{t} .$

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

 
    Regression with ARMA(0,1) Error Model with Seasonal MA(12) (Gaussian Distribution):
 
                   Value      StandardError    TStatistic      PValue  
                 _________    _____________    __________    __________

    Intercept    0.0042146      0.0015333        2.7486       0.0059843
    MA{1}         -0.47713        0.12087       -3.9473      7.9044e-05
    SMA{12}       -0.74115        0.12042       -6.1549      7.5116e-10
    Beta(1)      -0.018912      0.0075665       -2.4994        0.012439
    Variance     0.0016695     0.00031169        5.3564      8.4901e-08

Используйте предполагаемые коэффициенты модели (содержавшийся в EstMdl) симулировать 25 реализации авиапассажира рассчитывает по 60-месячному горизонту. Выведите остаточные значения и используйте их в качестве предварительной выборки.

[~,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

Figure contains an axes object. The axes object with title blank T r a n s f o r m e d , blank S i m u l a t e d blank M o n t h l y blank P a s s e n g e r blank T o t a l s contains 54 objects of type line. These objects represent Observations, Simulation Mean.

Модель регрессии с ошибками SMA, кажется, предсказывает ряд хорошо.

Проверьте прогнозирующую эффективность модели:

Варьируясь размер периода прогноза
Оценка среднеквадратичной погрешности предсказания (PMSE)
Выбор модели с самым низким PMSE

Неназванный

Скрипт Open Live Script

Симулируйте модель регрессии с неустановившимися мультипликативными сезонными ошибками

Загрузите наборы данных рецессий и авиакомпания. Преобразуйте ответ.

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

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

Задайте наборы индекса тот раздел данные в выборки прогноза и оценку.

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

Оцените модель регрессии с мультипликативными сезонными ошибками:

$y_{t} = X_{t} β + u_{t}$

$(1 - L) (1 - L^{12}) u_{t} = (1 + B L) (1 + B_{12} L^{12}) ε_{t} .$

Установите точку пересечения модели регрессии на 0, поскольку это не идентифицируется в интегрированной модели.

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

 
    Regression with ARIMA(0,1,1) Error Model Seasonally Integrated with Seasonal MA(12) (Gaussian Distribution):
 
                   Value      StandardError    TStatistic      PValue  
                 _________    _____________    __________    __________

    Intercept            0              0           NaN             NaN
    MA{1}         -0.35662        0.10393       -3.4312      0.00060088
    SMA{12}       -0.67729        0.11294       -5.9972      2.0077e-09
    Beta(1)      0.0015098       0.020533       0.07353         0.94138
    Variance     0.0015198     0.00021411        7.0983      1.2631e-12

Используйте предполагаемые коэффициенты модели (содержавшийся в EstMdl), чтобы симулировать количества авиапассажира по 60-месячному горизонту. Выведите остаточные значения и используйте их в качестве предварительной выборки.

[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);

figure
h1 = plot(dates(estInds),y(estInds));
title('{\bf Simulated Monthly Passenger Totals}')
hold on
plot(dates(foreInds),ySim,'Color',[.85,.85,.85])
h2 = plot(dates(foreInds),meanYSim,'k.-','LineWidth',2);
plot([dates(estInds(end)),dates(foreInds(1))],...
    [repmat(y(estInds(end)),numPaths,1),ySim(1,:)'],...
    'Color',[.85,.85,.85])
plot([dates(estInds(end)),dates(foreInds(1))],...
    [y(estInds(end)),meanYSim(1)],'k.-','LineWidth',2)
plot(dates(foreInds),y(foreInds))
datetick
legend([h1,h2],'Observations','Monte Carlo Forecasts',...
    'Location','NorthWest')
axis tight
hold off

Figure contains an axes object. The axes object with title blank S i m u l a t e d blank M o n t h l y blank P a s s e n g e r blank T o t a l s contains 1004 objects of type line. These objects represent Observations, Monte Carlo Forecasts.

Симулированные прогнозы показывают рост и сезонную периодичность, похожую на наблюдаемый ряд. Модель регрессии с ошибками SMA, кажется, предсказывает ряд хорошо, хотя немного переоценивая.

Проверьте прогнозирующую эффективность модели:

Варьируясь размер периода прогноза
Оценка среднеквадратичной погрешности предсказания (PMSE)
Выбор модели с самым низким PMSE

Документация