Сезонная регулировка с использованием S (n, m) сезонных фильтров

Открыть Live Script

В этом примере показано, как применить $S_{n \times m}$ сезонные фильтры для десеасонализации временных рядов (с помощью мультипликативного разложения). Временные ряды - ежемесячные международные счетчики пассажиров авиакомпаний с 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

Figure contains an axes. The axes with title Airline Passenger Counts contains an object of type line.

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

Детрендируйте данные с помощью 13-срочного скользящего среднего значения.

Перед оценкой сезонного компонента оцените и удалите линейный тренд. Примените 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

Figure contains an axes. The axes with title Airline Passenger Counts contains 2 objects of type line. These objects represent Passenger counts, 13-Term Moving Average.

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

Создайте массив ячеек, 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

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

Применить фильтр S (3,3

Применить 5-термин $S_{3 \times 3}$ сезонное скользящее среднее значение к детрендированному ряду 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';

Figure contains an axes. The axes with title Estimated Seasonal Component contains an object of type line.

Заметьте, что сезонный уровень изменяется в области значений данных. Это иллюстрирует различие между $S_{n \times m}$ сезонный фильтр и стабильный сезонный фильтр. Стабильный сезонный фильтр принимает, что сезонный уровень является постоянным в области значений данных.

Примените 13-терминный фильтр Хендерсона.

Чтобы получить улучшенную оценку компонента тренда, примените 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

Figure contains an axes. The axes with title Airline Passenger Counts contains 2 objects of type line. These objects represent Passenger counts, 13-Term Henderson Filter.

Применить S (3,5) сезонный фильтр

Чтобы получить 6. улучшенная оценка сезонного компонента, применить 7-термин $S_{3 \times 5}$ сезонное скользящее среднее значение к вновь детрендированной серии. Симметричные и асимметричные веса приведены в следующем коде. Центрируйте сезонную оценку, чтобы колебаться около 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 contains an axes. The axes with title Deseasonalized Airline Passenger Counts contains an object of type line.

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

Постройте график компонентов и исходного ряда.

Сравните исходную серию с серией, восстановленной с помощью оценок компонентов.

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

Figure contains an axes. The axes with title Airline Passenger Counts contains 3 objects of type line. These objects represent Passenger counts, 13-Term Henderson Filter, Trend and Seasonal Components.

Оцените неправильный компонент.

Детренд и десеасонализация оригинальной серии. Постройте график оставшейся оценки неправильного компонента.

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';

Figure contains an axes. The axes with title Airline Passenger Counts Irregular Component contains an object of type line.

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

См. также

cellfun | conv | conv2

Документация