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

Этот пример показывает, как применяться $S_{n \times m}$ сезонные фильтры к deseasonalize временные ряды (использующий мультипликативное разложение). Временные ряды являются ежемесячными международными количествами авиапассажира от 1 949 до 1960.

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

Загрузите набор данных авиакомпании.

load(fullfile(matlabroot,'examples','econ','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 терминами.

Прежде, чем оценить сезонный компонент, оцените и удалите линейный тренд. Примените симметричное скользящее среднее значение с 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;

h = plot(yS,'r','LineWidth',2);
legend(h,'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

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

Примените 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';

Заметьте, что сезонный уровень переключает область значений данных. Это иллюстрирует различие между $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('13-Term Henderson Filter')
hold off

Примените 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';

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('Original Series','13-Term Henderson Filter',...
       'Trend and Seasonal Components')
hold off

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

Детрендируйте и 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';

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

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

cellfun | conv | conv2

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