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

Скрипт Open Live Script

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

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

% 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 object. The axes object with title Deseasonalized Airline Passenger Counts contains an object of type line.

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

Постройте компоненты и исходный ряд.

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

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 object. The axes object with title Airline Passenger Counts contains 3 objects of type line. These objects represent Passenger counts, 13-Term Henderson Filter, Trend and Seasonal Components.

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

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

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 object. The axes object with title Airline Passenger Counts Irregular Component contains an object of type line.

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

Смотрите также

conv | conv2 | cellfun

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