Шаг 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);

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

Ссылки:

Бокс, Г. Э. П., Г. М. Дженкинс и Г. К. Рейнсель. Анализ временных рядов: прогнозирование и контроль. 3-й ред. Энглвуд Клиффс, Нью-Джерси: Прентис Холл, 1994.