Этот пример показывает, как оценить несезонные и сезонные компоненты тренда с помощью параметрических моделей. Временные ряды являются ежемесячными смертями от несчастного случая в США от 1 973 до 1978 (Броквелл и Дэвис, 2002).
Загрузите набор данных смертей от несчастного случая.
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.
Соответствуйте полиному
к наблюдаемому ряду.
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
Вычтите подходящую квадратичную строку из исходных данных.
xt = y - tH;
Создайте индикатор (макет) переменные в течение каждого месяца. Первый индикатор равен одному для наблюдений в январе и нуля в противном случае. Второй индикатор равен одному для наблюдений в феврале и нуля в противном случае. В общей сложности 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 индикаторов и термин прерывания).
Вычтите предполагаемый сезонный компонент из исходного ряда.
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);
Квадратичный тренд намного более ясен с сезонным удаленным компонентом.
Вычтите тренд и сезонные оценки от исходного ряда. Остаток является оценкой неправильного компонента.
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);
Можно опционально смоделировать неправильный компонент с помощью модели стохастического процесса.
Ссылки:
Поле, G. E. P. Г. М. Дженкинс и Г. К. Рейнсель. Анализ timeseries: Прогнозирование и Управление. 3-й редактор Englewood Cliffs, NJ: Prentice Hall, 1994.