В этом примере показано, как применить сезонные фильтры для десеасонализации временных рядов (с помощью мультипликативного разложения). Временные ряды - ежемесячные международные счетчики пассажиров авиакомпаний с 1949 по 1960 год.
Загрузите набор данных авиакомпании.
load('Data_Airline.mat') y = Data; T = length(y); figure plot(y) h1 = gca; h1.XLim = [0,T]; h1.XTick = 1:12:T; h1.XTickLabel = datestr(dates(1:12:T),10); title 'Airline Passenger Counts'; hold on
Данные показывают восходящий линейный тренд и сезонный компонент с периодичностью 12.
Перед оценкой сезонного компонента оцените и удалите линейный тренд. Примените 13-терминовое симметричное скользящее среднее значение, повторяя первое и последнее наблюдения шесть раз, чтобы предотвратить потерю данных. Используйте вес 1/24 для первого и последнего терминов скользящего среднего значения и вес 1/12 для всех интерьерных терминов.
Разделите исходный ряд на сглаженные ряды, чтобы удалить данные. Добавьте оценку скользящего среднего тренда к наблюдаемому графику временных рядов.
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; plot(yS,'r','LineWidth',2) legend(["Passenger counts" "13-Term Moving Average"]) hold off
Создайте массив ячеек, sidx, для хранения индексов, соответствующих каждому периоду. Данные являются ежемесячными, с периодичностью 12, поэтому первый элемент sidx является вектором с элементами 1, 13, 25,..., 133 (соответствующим январским наблюдениям). Второй элемент sidx - вектор с элементами 2, 14, 16,..., 134 (соответствующий февральским наблюдениям). Это повторяется в течение всех 12 месяцев .
s = 12; sidx = cell(s,1); % Preallocation for i = 1:s sidx{i,1} = i:s:T; end sidx{1:2}
ans = 1×12
1 13 25 37 49 61 73 85 97 109 121 133
ans = 1×12
2 14 26 38 50 62 74 86 98 110 122 134
Использование массива ячеек для хранения индексов допускает возможность того, что каждый период не происходит одинаковое количество раз в пределах диапазона наблюдаемых рядов.
Применить 5-термин сезонное скользящее среднее значение к детрендированному ряду xt. То есть применить скользящее среднее значение к значениям января (по индексам 1, 13, 25,..., 133), а затем применить скользящее среднее значение к серии февраля (по индексам 2, 14, 26,..., 134) и так далее на оставшиеся месяцы.
Используйте асимметричные веса в концах скользящего среднего значения (используя conv2). Поместите сглаженные значения назад в один вектор.
Чтобы центрировать сезонный компонент вокруг единицы, оцените, а затем разделите на 13-кратное скользящее среднее значение расчетного сезонного компонента.
% S3x3 seasonal filter % Symmetric weights sW3 = [1/9;2/9;1/3;2/9;1/9]; % Asymmetric weights for end of series aW3 = [.259 .407;.37 .407;.259 .185;.111 0]; % Apply filter to each month shat = NaN*y; for i = 1:s ns = length(sidx{i}); first = 1:4; last = ns - 3:ns; dat = xt(sidx{i}); sd = conv(dat,sW3,'same'); sd(1:2) = conv2(dat(first),1,rot90(aW3,2),'valid'); sd(ns -1:ns) = conv2(dat(last),1,aW3,'valid'); shat(sidx{i}) = sd; end % 13-term moving average of filtered series sW13 = [1/24;repmat(1/12,11,1);1/24]; sb = conv(shat,sW13,'same'); sb(1:6) = sb(s+1:s+6); sb(T-5:T) = sb(T-s-5:T-s); % Center to get final estimate s33 = shat./sb; figure plot(s33) h2 = gca; h2.XLim = [0,T]; h2.XTick = 1:12:T; h2.XTickLabel = datestr(dates(1:12:T),10); title 'Estimated Seasonal Component';
Заметьте, что сезонный уровень изменяется в области значений данных. Это иллюстрирует различие между сезонный фильтр и стабильный сезонный фильтр. Стабильный сезонный фильтр принимает, что сезонный уровень является постоянным в области значений данных.
Чтобы получить улучшенную оценку компонента тренда, примените 13-терминный фильтр Хендерсона к сезонно скорректированному ряду. Необходимые симметричные и асимметричные веса приведены в следующем коде.
% Deseasonalize series dt = y./s33; % Henderson filter weights sWH = [-0.019;-0.028;0;.066;.147;.214; .24;.214;.147;.066;0;-0.028;-0.019]; % Asymmetric weights for end of series aWH = [-.034 -.017 .045 .148 .279 .421; -.005 .051 .130 .215 .292 .353; .061 .135 .201 .241 .254 .244; .144 .205 .230 .216 .174 .120; .211 .233 .208 .149 .080 .012; .238 .210 .144 .068 .002 -.058; .213 .146 .066 .003 -.039 -.092; .147 .066 .004 -.025 -.042 0 ; .066 .003 -.020 -.016 0 0 ; .001 -.022 -.008 0 0 0 ; -.026 -.011 0 0 0 0 ; -.016 0 0 0 0 0 ]; % Apply 13-term Henderson filter first = 1:12; last = T-11:T; h13 = conv(dt,sWH,'same'); h13(T-5:end) = conv2(dt(last),1,aWH,'valid'); h13(1:6) = conv2(dt(first),1,rot90(aWH,2),'valid'); % New detrended series xt = y./h13; figure plot(y) h3 = gca; h3.XLim = [0,T]; h3.XTick = 1:12:T; h3.XTickLabel = datestr(dates(1:12:T),10); title 'Airline Passenger Counts'; hold on plot(h13,'r','LineWidth',2); legend(["Passenger counts" "13-Term Henderson Filter"]) hold off
Чтобы получить 6. улучшенная оценка сезонного компонента, применить 7-термин сезонное скользящее среднее значение к вновь детрендированной серии. Симметричные и асимметричные веса приведены в следующем коде. Центрируйте сезонную оценку, чтобы колебаться около 1.
Десеасонализуйте исходный ряд путем деления его на центрированную сезонную оценку.
% S3x5 seasonal filter % Symmetric weights sW5 = [1/15;2/15;repmat(1/5,3,1);2/15;1/15]; % Asymmetric weights for end of series aW5 = [.150 .250 .293; .217 .250 .283; .217 .250 .283; .217 .183 .150; .133 .067 0; .067 0 0]; % Apply filter to each month shat = NaN*y; for i = 1:s ns = length(sidx{i}); first = 1:6; last = ns-5:ns; dat = xt(sidx{i}); sd = conv(dat,sW5,'same'); sd(1:3) = conv2(dat(first),1,rot90(aW5,2),'valid'); sd(ns-2:ns) = conv2(dat(last),1,aW5,'valid'); shat(sidx{i}) = sd; end % 13-term moving average of filtered series sW13 = [1/24;repmat(1/12,11,1);1/24]; sb = conv(shat,sW13,'same'); sb(1:6) = sb(s+1:s+6); sb(T-5:T) = sb(T-s-5:T-s); % Center to get final estimate s35 = shat./sb; % Deseasonalized series dt = y./s35; figure plot(dt) h4 = gca; h4.XLim = [0,T]; h4.XTick = 1:12:T; h4.XTickLabel = datestr(dates(1:12:T),10); title 'Deseasonalized Airline Passenger Counts';
Десеасонализированный ряд состоит из долгосрочного тренда и нерегулярных компонентов. С удаленного сезонного компонента легче увидеть переломные точки в тренде.
Сравните исходную серию с серией, восстановленной с помощью оценок компонентов.
figure plot(y,'Color',[.85,.85,.85],'LineWidth',4) h5 = gca; h5.XLim = [0,T]; h5.XTick = 1:12:T; h5.XTickLabel = datestr(dates(1:12:T),10); title 'Airline Passenger Counts'; hold on plot(h13,'r','LineWidth',2) plot(h13.*s35,'k--','LineWidth',1.5) legend(["Passenger counts" "13-Term Henderson Filter" ... "Trend and Seasonal Components"]) hold off
Детренд и десеасонализация оригинальной серии. Постройте график оставшейся оценки неправильного компонента.
Irr = dt./h13;
figure
plot(Irr)
h6 = gca;
h6.XLim = [0,T];
h6.XTick = 1:12:T;
h6.XTickLabel = datestr(dates(1:12:T),10);
title 'Airline Passenger Counts Irregular Component';
Можно опционально смоделировать детрендированный и десесонализированный ряд с помощью стационарной стохастической модели процесса.