Загрузите данные.

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

load('Data_Accidental.mat')
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 (thousands)';
hold on

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

Применить 13-терминовое скользящее среднее значение.

Сглаживайте данные с помощью 13-срочного скользящего среднего значения. Чтобы предотвратить потерю наблюдения, повторите первое и последнее сглаженные значения шесть раз. Вычесть сглаженный ряд из исходного ряда, чтобы искажать данные. Добавьте оценку скользящего среднего тренда к наблюдаемому графику временных рядов.

sW13 = [1/24;repmat(1/12,11,1);1/24];
yS = conv(y,sW13,'same');
yS(1:6) = yS(7); yS(T-5:T) = yS(T-6);

xt = y-yS;

h = plot(yS/1000,'r','LineWidth',2);
legend(h,'13-Term Moving Average')
hold off

Временные ряды без тренда xt.

Использование параметра shape 'same' при вызове conv возвращает сглаженный ряд той же длины, что и исходный ряд.

Шаг 3. Создайте сезонные индексы.

Создайте массив ячеек, sidx, для хранения индексов, соответствующих каждому периоду. Данные являются ежемесячными, с периодичностью 12, поэтому первый элемент sidx является вектором с элементами 1, 13, 25,..., 61 (соответствующим январским наблюдениям). Второй элемент sidx - вектор с элементами 2, 14, 16,..., 62 (соответствующий февральским наблюдениям). Это повторяется в течение всех 12 месяцев .

s = 12;
sidx = cell(s,1);
for i = 1:s
 sidx{i,1} = i:s:T;
end

sidx{1:2}

ans = 1×6

     1    13    25    37    49    61

ans = 1×6

     2    14    26    38    50    62

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

Шаг 4. Применить стабильный сезонный фильтр.

Примените стабильный сезонный фильтр к детрендированной серии xt. Используя индексы, созданные на этапе 3, усредняйте данные без тренда, соответствующие каждому периоду. То есть, среднее значение все значения января (по индексам 1, 13, 25,..., 61), а затем среднее значение все значения февраля (по индексам 2, 14, 26,..., 62) и так далее для остальных месяцев. Поместите сглаженные значения назад в один вектор.

Центрируйте сезонную оценку, чтобы колебаться около нуля.

sst = cellfun(@(x) mean(xt(x)),sidx);

% Put smoothed values back into a vector of length N
nc = floor(T/s); % no. complete years
rm = mod(T,s); % no. extra months
sst = [repmat(sst,nc,1);sst(1:rm)];

% Center the seasonal estimate (additive)
sBar = mean(sst); % for centering
sst = sst-sBar;

figure
plot(sst/1000)
title 'Stable Seasonal Component';
h2 = gca;
h2.XLim = [0 T];
ylabel 'Number of deaths (thousands)';
h2.XTick = 1:12:T;
h2.XTickLabel = datestr(dates(1:12:T),10);

Стабильный сезонный компонент имеет постоянную амплитуду в рядах. Сезонная оценка центрирована и колеблется около нуля.

Шаг 5. Десеасонализируйте серию.

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

dt = y - sst;

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

Десеасонализированный ряд состоит из долгосрочного тренда и нерегулярных компонентов. Широкомасштабный квадратичный тренд в количестве случайных смертей очевиден при удалении сезонного компонента.

Ссылки:

Броквелл, П. Дж. и Р. А. Дэвис. Введение во временные ряды и прогнозирование. 2nd ed. New York, NY: Springer, 2002.