В этом примере показано, как применяться сезонные фильтры к 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
Данные показывают восходящий линейный тренд и сезонный компонент с периодичностью 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.
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(["Passenger counts" "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 ряд с помощью стационарной модели стохастического процесса.