Шаг 1: Загрузка данных

Загрузите набор данных о случайных смертях.

load Data_Accidental 
y = Data;
T = length(y);

figure
plot(y/1000)
h1 = gca;
h1.XLim = [0,T];
h1.XTick = 1:12:T;
h1.XTickLabel = datestr(dates(1:12:T),10);
title 'Monthly Accidental Deaths';
ylabel 'Number of Deaths (in thousands)';
hold on

Данные показывают потенциальный квадратичный тренд и сильный сезонный компонент с периодичностью 12.

Шаг 2: подгонка квадратичного тренда

Аппроксимация полинома

к наблюдаемой серии.

t = (1:T)';
X = [ones(T,1) t t.^2];

b = X\y;
tH = X*b;
 
h2 = plot(tH/1000,'r','LineWidth',2);
legend(h2,'Quadratic Trend Estimate')
hold off

Шаг 3. Детренд оригинальной серии.

Вычесть подобранную квадратичную линию из исходных данных.

xt = y - tH;

Шаг 4. Оценка переменных сезонного индикатора

Создайте переменные индикатора (фиктивные) для каждого месяца. Первый показатель равен единице для январских наблюдений и нулю в противном случае. Второй показатель равен единице для февральских наблюдений и нулю в противном случае. Всего за 12 месяцев создается 12 переменных индикатора. Регрессируйте детрендированные ряды по сезонным показателям.

mo = repmat((1:12)',6,1);
sX = dummyvar(mo);
  
bS = sX\xt;
st = sX*bS;

figure
plot(st/1000)
title 'Parametric Estimate of Seasonal Component (Indicators)';
h3 = gca;
h3.XLim = [0,T];
ylabel 'Number of Deaths (in thousands)';
h3.XTick = 1:12:T;
h3.XTickLabel = datestr(dates(1:12:T),10);

В этой регрессии все 12 сезонных показателей включены в матрицу проекта. Чтобы предотвратить коллинеарность, термин точки пересечения не включается (кроме того, можно включать 11 индикаторов и термин точки пересечения).

Шаг 5. Десеасонализация оригинальных серий

Вычесть расчетный сезонный компонент из исходного ряда.

dt = y - st;

figure
plot(dt/1000)
title 'Monthly Accidental Deaths (Deseasonalized)';
h4 = gca;
h4.XLim = [0,T];
ylabel 'Number of Deaths (in thousands)';
h4.XTick = 1:12:T;
h4.XTickLabel = datestr(dates(1:12:T),10);

Квадратичный тренд намного яснее с удаленным сезонным компонентом.

Шаг 6. Оценка нерегулярного компонента

Вычесть тренд и сезонные оценки из исходного ряда. Остальная часть является оценкой нерегулярного компонента.

bt = y - tH - st;

figure
plot(bt/1000)
title('Irregular Component')
h5 = gca;
h5.XLim = [0,T];
ylabel 'Number of Deaths (in thousands)';
h5.XTick = 1:12:T;
h5.XTickLabel = datestr(dates(1:12:T),10);

Вы можете опционально смоделировать нерегулярный компонент с помощью стохастической модели процесса.

Ссылки:

Box, G. E. P., G. M. Jenkins, and G. C. Reinsel. Анализ временных рядов: прогнозирование и управление. 3-й эд. Englewood Cliffs, Нью-Джерси: Prentice Hall, 1994.