Z Skrypty dla studentów Ekonofizyki UPGOW

Aktualna wersja na dzień 11:56, 10 lut 2011

Analiza Szeregów Czasowych
<<< Procesy stochastyczne | Modelowanie szeregów czasowych >>>

Indeks Dow Jones Industrial Average w latach 1928-1930

Właściwa analiza dowolnego szeregu czasowego zaczyna się od wizualizacji danych. Znaczy to, że na samym początku powinniśmy stworzyć wykres badanych danych. Daje nam to możliwość oszacowania jak konkretny zbiór danych zachowuje się w czasie. Zazwyczaj możemy ocenić istnienie trendu czy sezonowości. W trudniejszych przypadkach zadecydować, gdy nie podzielić danego szeregu na dwie czy więcej części i badać je osobno. Przykładowo można rozpatrywać zmienność indeksów giełdowych przed krachem i po nim (np: przed i po Czarnym czwartku), opisując serię danych dwoma różnymi modelami. Dodatkowo mamy możliwość decyzji czy przypadkiem część danych nie została uzyskanych przypadkowo (np: błąd pomiaru, przypadkowy pomiar innych wielkości).

Klasyczna dekompozycja sygnału

W szczególności oględziny wykresu analizowanych danych dają możliwość (po pewnej wprawie) na decyzję, czy możemy daną serię danych reprezentować jako realizację procesu

\( X_t = m_t + s_t + Y_t, \!\)

gdzie

• \(X_t\) to dane pomiarowe,
• \(m_t\) to wielkość opisująca wolno zmienną funkcję, czyli trend,
• \(s_t\) to periodycznie zmienna funkcja zwana sezonowością,
• \(Y_t\) to komponent opisujący (stacjonarny) proces losowy, często nazywany szumem.

Celem dekompozycji szeregu czasowego jest oszacowanie i ekstrakcja deterministycznych części szeregu - trendu \(m_t\) oraz sezonowości \(s_t\) w nadziei, że pozostałe dane, czyli teoretycznie zmienna losowa \(Y_t\) okaże się stacjonarnym procesem losowym. W przypadku, kiedy okaże się to prawdą, tj. reszty \(Y_t\) mogą być opisane stacjonarnym procesem losowym \(\{Y_t\}\), możemy przystąpić do przewidywania przyszłego zachowania się szeregu, wykorzystując oczywiście wszystkie posiadane wiadomości: trend, okres oraz zidentyfikowany z pewną dokładnością proces losowy.

Metoda Boxa-Jenkinsa

Inną, alternatywną metodą do opisanej wyżej, jest metoda która dopasowywuje modele ARMA i ARIMA do istniejących danych. Metoda ta została nazwana po nazwiskach dwóch statystyków Georgea Boxa oraz Gwilyma Jenkinsa, którzy rozwinęli tą metodę w latach 70-tych.

Algorytm Boxa-Jenkinsa składa się z trzech kroków:

1. Identyfikacja oraz wybór modelu: pierwszym krokiem jest upewnienie się, że analizujemy dane stacjonarne; następnie identyfikujemy sezonowość i usuwamy ją z danych aby w końcu wykorzystując wykresy funkcji autokorelacji oraz częściowej autokorelacji zdecydować jakie komponenty AR (autoregresji), I (scałkowane) lub MA (średniej ruchomej) wykorzystać do budowy modelu.
2. Znalezienie parametrów wybranego modelu za pomocą wybranych metod (numerycznych) tak, aby dopasowanie danych do modelu było najlepsze. Najczęstszymi metodami wykorzystywanymi w praktyce są: maximum likelihood estimation lub (nieliniowa) metoda najmniejszych kwadratów.
3. Sprawdzenie poprawności wyboru danego modelu. W szczególności należy sprawdzić czy proces jest stacjonarny - reszty muszą być od siebie niezależne, oraz ich średnia i wariancja musi być stała w czasie. Można
   • narysować wykres średniej, wariancji oraz reszt versus czas (indeks) i przeprowadzić na nich test Ljunga-Boxa,
   • narysować wykresy funkcji autokorelacji i częściowej autokorelacji reszt.

Jeżeli oszacowanie jest niedobre, należy wrócić na początek algorytmu i postarać się polepszyć istniejący model lub zbudować nowy.

Eliminacja Trendu

Jeżeli analizowane dane nie wykazują sezonowości, ogólny model zastąpić możemy poprzez uproszczony

\( X_t = m_t + Y_t, \!\)

\( t = 1, 2, ..., n\). Zakładamy dodatkowo, bez straty ogólności, że \( E Y_t = 0.\)

Rysunek 4.2. Liczba ludności USA w latach 1790 - 1980, w odstępach dziesięcioletnich.

Metoda 1 (metoda minimum sumy kwadratów błędów)

Potoczna nazwa (i obecnie w pełni przyjęta) to metoda najmniejszych kwadratów (MNK). Służy do estymacji i wyznaczania linii trendu na podstawie zbioru danych w postaci par liczb. Najczęściej jest stosowana przy regresji liniowej, ale może też być stosowana do statystycznego wyznaczania parametrów nieliniowych linii trendu.

Przykład 4.1: Liczba ludności USA w latach 1790 - 1980, w odstępach dziesięcioletnich.

Po wstępnej analizie wykresu danych populacji USA z lat 1790 - 1980 (patrz wstęp) zakładamy, że nasze dane mogą być dopasowane funkcją potęgową drugiego rzędu

\( m_t = a_2 t^2 + a_1 t + a_0,\)

poprzez odpowiedni dobór parametrów modelu. Dobór parametrów następuje poprzez minimalizację funkcji

\(\chi^2(a_2,a_1,a_0) = \sum_{t=1790}^{1980} (X_t - m_t)^2 = \sum_{t=1790}^{1980} (X_t - a_2 t^2 - a_1 t - a_0)^2. \)

Rysunek 4.1. Liczba ludności USA w latach 1790 - 1980, w odstępach dziesięcioletnich.

Prosty rachunek (w Octave) daje nam szacowanie parametrów:

\( a_0 = 2.0979 \cdot 10^{10}, \)

\( a_1 = -2.3350 \cdot 10^{7}, \)

\( a_2 = 6.4986 \cdot 10^{3}. \)

Na rysunku 4.1 przedstawiono dane wraz z dofitowanymi metodą najmniejszych kwadratów wielomianami rzędu 1 i 2. Znając wartości parametrów modelu, a więc mając określony wielomian możemy pokusić się o dalszą analizę szeregu, mianowicie o ocenę reszt modelu. Tym zagadnieniem zajmiemy się nieco później. Na razie wystarczy, jeżeli narysujemy reszty, \(Y_t = X_t - m_t\) i poddamy je ocenie wzrokowej (możemy też łatwo policzyć średnią).

Przykład 

Używając metody najmniejszych kwadratów, dopasujemy teraz wielomian 1 i 2-go rzędu do danych "Liczba ludności USA w latach 1790 - 1980".

close all
clear all
y = csvread('PopulacjaUSA1790-1980.csv');
X = y(:, 1);
Y = y(:, 2) / 10^6;
 
plot(X, Y, "sr; Dane;", 
         'MarkerSize', 14, 'markeredgecolor', 'black');
 
x = X(1) : X(length(X));
txt = "gbk"
 
for (deg = 1 : 2)
  printf("NMK Stopień: %d\n", deg);
  printf("Współczynniki (od najwyższej potęgi):\n");
  a = polyfit(X, Y, deg)
  hold on;
  plot(x, polyval(a,x), "linewidth", 2, 
         sprintf("--%s; MNK rzedu %d;", txt(deg), deg));
endfor;
 
xlabel('t');
ylabel('X_t (w milionach)');
grid on;
title('Populacja USA w latach 1790-1980');
legend('Location', 'SouthEast');
print('cw41.png', '-dpng', '-FTimes-Roman:24')
 
%%%
% Obliczanie reszt
close all;
Z = Y - polyval(a, X);
mZ = mean(Z);
plot(X, Z, "linewidth", 1, "--sb; reszty;");
hold on;
 
%%%
% srednia
plot(X, 0*X + mZ, "linewidth", 2, sprintf("-k; srednia = %.4f;",mZ));
 
xlabel('t');
ylabel('x_t - m_t (w milionach)');
grid on;
legend('Location', 'SouthWest');
print('cw41reszty.png', '-dpng', '-FTimes-Roman:24')

Metoda 2 (wygładzanie krzywych)

Metoda średniej kroczącej

Przykład wygładzania 7-punktową i 11-punktową średnią kroczącą

Niech q będzie nieujemną liczbą całkowitą. Rozważmy dwustronną średnią kroczącą procesu \( \{ X_t = m_t + Y_t\} \ \), zdefiniowaną jako

\( W_t = (2q - 1)^{-1} \sum_{j=-q}^q X_{t+j}. \ \)

Zakładając, że na wybranym przedziale \([t-q,t+q] \ \), trend \(m_t \ \) jest w dobrym przybliżeniu liniowy oraz średnia reszt modelu \(Y_t \ \) jest bliska zeru, wówczas, dla \(q + 1 \le t \le n - q \ \) powyższe równanie możemy zapisać jako

\( W_t = (2q - 1)^{-1} \sum_{j=-q}^q m_{t+j} + (2q - 1)^{-1} \sum_{j=-q}^q Y_{t+j} \simeq m_t. \ \)

Metoda średniej kroczącej pozwala nam oszacować trend

\( \hat{m}_t = (2q-1)^{-1} \sum_{j=-q}^q X_{t+j}, ~~~q + 1 \le t \le n - q. \ \)

Ponieważ \( \{ X_t \} \ \) reprezentuje szereg czasowy nie posiadający rejestru (danych) dla indeksu czasowego mniejszego niż 0 oraz dla indeksu czasowego większego niż rozmiar próbki (n), dlatego metody tej nie da się zastosować dla początkowego i końcowego fragmentu szeregu, tj. gdy \( t \le q \ \) oraz \( t > n -q \ \).

function [z] = asc_moving_average(x, c)
 
  if ( isempty (x) || (!isvector(x)) ) 
    error ("movingAverage: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if (mod(c,2) == 0 || c < 1 || c > length(x))
    error("s=%d, Podaj nieparzyste s >= 1\n",c);
    return;
  endif
 
 q = (c-1)/2;
 for i = (q+1):(length(x)-q)
   z(i-q) = mean(x(i-q:i+q));
 endfor
 
endfunction

Wygładzanie wykładnicze

Metoda ta pozwala zmniejszyć wariancję oryginalnego szeregu czasowego za pomocą ważonej średniej ruchomej z przeszłych wartości, o wagach malejących wykładniczo wraz z odległością w czasie.

Proste wygładzanie wykładnicze (model Browna)

Dla dowolnej liczby \( a \in [0,1] \ \) możemy policzyć jednostronną średnią kroczącą \( \hat{m}_t, t = 1,\dots,n \ \) zdefiniowaną poprzez rekurencje (czasami zwane rekursjami)

\( \begin{align} \hat{m}_t &= a X_t + (1 - a) \hat{m}_{t-1}, ~~~t=2,\dots,n, \\ \hat{m}_1 &= X_1. \end{align} \)

Łatwo można policzyć, że dla \( t \ge 2 \ \) poszczególne \( \hat{m}_t \) są

\( \hat{m}_t = \sum_{j=0}^{t-2} a (1 - a)^j x_{t-j} + (1 - a)^{t-1} X_1, \ \)

czyli, że są średnimi kroczącymi z wagami zanikającymi wykładniczo (za wyjątkiem ostatniego wyrazu).

Wybór parametru a odbywa się zazwyczaj metodą prób i błędów. Poniżej czytelnik znajdzie program w języku Matlab/Octave realizujący wygładzanie wykładnicze.

function [m] = asc_exponential_smoothing(x, a = 0.5)
 
  if ( isempty (x) || (!isvector(x)) ) 
    error ("movingAverage: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if (a < 0 || a > 1)
    error("a=%d, Podaj 0 < a < 1\n",a);
    return;
  endif
 
  m(1) = x(1);
  for t = 2 : length(x)
     m(t) = a * x(t) + (1 - a) * m(t - 1);
  endfor
 
endfunction

Model liniowy Holta

Model ten stosowany jest do wygładzania szeregów czasowych z trendem. Model opisujemy dwoma równaniami:

\(F_{t} = \alpha X_{t} + \left( {1 - \alpha } \right)\left( {F_{t - 1} + S_{t - 1} } \right)\)

\(S_{t} = \beta \left( {F_{t} - F_{t - 1} } \right) + \left( {1 - \beta } \right)S_{t - 1}\)

• \(F_{t}\,\) - wygładzona wartość zmiennej w chwili \(t\)
• \(S_{t}\,\) - wygładzona wartość przyrostu (trendu) w chwili \(t\)
• \(\alpha, \beta\,\) - parametry modelu o wartościach z przedziału [0,1]
  

Wartości \(\alpha, \beta\,\), podobnie jak parametr a poprzedniej metody znajdujemy metodą prób i błędów. Wartości początkowe \(F_1\,\) i \(S_1\,\) można uzyskać na dwa sposoby:

• z użyciem pierwszych wyrazów wektora danych
  • \(F_1\,\) = \(X_1\,\),
  • \(S_1\,\) = \(X_2-X_1\,\),
• z użyciem dopasowania funkcji liniowej do trendu (szeregu)
  • \(F_1\,\) wyrazu wolnego liniowej funkcji trendu (po uprzednim dopasowaniu),
  • \(S_1\,\) współczynnik kierunkowy tej samej funkcji.

function [F, S] = asc_exponential_smoothing_Holt(x, a = 0.5, b = 0.5, initial = "vec")
 
  length_x = length(x);
 
  if ( isempty (x) || (!isvector(x)) || length_x < 2) 
    error ("asc_exponential_smoothing_Holt: \
    pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if (a < 0 || a > 1)
    error("asc_exponential_smoothing_Holt: \
    a=%d, Podaj 0 < a < 1\n",a);
    return;
  endif
 
  if (b < 0 || b > 1)
    error("asc_exponential_smoothing_Holt: \
    b=%d, Podaj 0 < b < 1\n",b);
    return;
  endif
 
  if initial == "vec"
    F(1) = x(1);
    S(1) = x(2) - x(1);
  elseif initial == "lin"
    [p1, p2, p3] = polyfit([1 : length_x]', x, 1);
    F(1) = p1(1);
    S(1) = p1(2);
  else
    error("expecting 'vec' or 'lin' for initial values, '%s' given", initial);
    return;
  endif
 
  for t = 2 : length_x
    F(t) = a * x(t) + (1 - a) * (F(t - 1) + S(t - 1));
    S(t) = b * (F(t) - F(t - 1)) + (1 - b) * S(t - 1);
  endfor
 
endfunction

Filtr liniowy

Użytecznie jest myśleć o procesie \( \{ m_t \} \ \) jako o procesie potomnym otrzymanym z procesu pierwotnego \( \{ X_t \} \ \) poprzez zastosowanie operatora liniowego bądź też filtra liniowego

\( \hat{m}_t = \sum_{j=-\infty}^{\infty} a_j X_{t+j}, \ \)

z wagami

\( \begin{align} a_j &= (2q + 1)^{-1}, ~~~ -q \le j \le q, \\ a_j &= 0, ~~~ |j| > q. \end{align} \)

Ten szczególny typ filtra jest filtrem dolnoprzepustowym, pozwala on bowiem usunąć elementy szybkozmienne (elementy o wysokiej częstotliwości) pozostawiając elementy wolno-zmienne (o częstotliwości niskiej) - takie jak trend (czy sezonowość).

Filtr zdefiniowany w paragrafie traktującym o średniej kroczącej to tylko jeden z wielu możliwych filtrów stosowanych do wygładzania danych. Dla dostatecznie dużych q czyli takich, by średnia szumu na badanym okresie była bliska zeru

\( (2q-1)^{-1} \sum_{j=-q}^{q} Y_{t+j} \simeq 0 \ \)

taki filtr przede wszystkim usunie szum z danych, ale za razem przepuści sygnał

\( m_t = at + b \ \)

który będzie liniową funkcją "czasu". Cały kłopot polega na takim dobraniu parametru q by obliczona właśnie średnia dalej była liniowa. Jeżeli q będzie za duże może się okazać, że przefiltrowany proces pomimo, że będzie ładnie wygładzony, nie będzie dobrym oszacowaniem procesu \( m_t \ \). Z drugiej strony - dobrze zaprojektowany filtr (czyli dobrze dobrany parametr q oraz wagi \( \{ a_i \} \ \)) pozwoli nie tylko na pozbycie się części losowej, ale również "przepuści" szeroką klasę funkcji trendu.

15 punktowa średnia krocząca Spencera

Przykład wygładzania 15-punktową średnia krocząca Spencera

Eliminacja wysokich częstości.

\( \begin{align} a_i &= 0, ~~~|i| > 7, \\ a_i &= a_{-i}, ~~~|i| \le 7 \\ \end{align} \)

a wagi zdefiniowane są następująco

\( [a_0, a_1, \dots, a_7] = \frac{1}{320} [74,67,46,21,3,-5,-6,-3]. \ \)

function [z] = asc_spencer(x)
 
  length_x = length(x);
 
  if ( isempty (x) || (!isvector(x)) || length_x < 18) 
    error ("\n\nasc_spencer:\npierwszy agrument musi byc niepustym \
    wektorem\no długości większej lub równej 18\npodano wektor od \
    długości: %d \n\n",length(x));
    return;
  endif
 
  %%%
  % wagi spencera
  spencer = [74, 67, 46, 21, 3, -5, -6, -3] / 320;
  length_spencer = length(spencer);  
 
  z = zeros(1, length_x - 2 * length_spencer);
  for i = length_spencer + 1 : length_x - length_spencer - 1
    for j = -length_spencer + 1 : length_spencer - 1
      z(i - length_spencer) += spencer(abs(j) + 1) * x(i + j);
    endfor
  endfor
 
endfunction

Metoda 3 (różnicowanie)

Rys.M3.1: szereg czasowy wykazujący trend

Rys.M3.2: reszty modelu \( Y_i \)

Metoda ta jest poniekąd odwrotnością metody 2. Zamiast usuwać szum z danych wejściowych usuniemy z nich trend i uzyskamy w ten sposób wektor szumu. Zastosujemy metodę różnicowania do danych stacjonarnych.

Na początku zdefiniujemy operator różnicowania jako

\( \nabla X_t = X_t - X_{t-1} = (1 - B) X_t, \ \)

gdzie \(B \ \) nazywamy operatorem przesunięcia w tył i definiujemy jako

\( B X_t = X_{t-1}. \ \)

Podamy teraz kilka właściwości operatora przesunięcia w tył \[ \begin{align} (i)~~ &B^j X_t = X_{t-j}, \\ (ii)~~ &\nabla^j X_t = \nabla(\nabla^{j-1} X_t),~~~ j \ge 1, \\ (iii)~~ &\nabla^0 X_t = X_t. \end{align} \]

Potęgi operatora \(\nabla \ \) zachowuja się tak samo jak wielomiany

\( \nabla^2 X_t = \nabla (\nabla (X_t)) = (1-B)(1-B) X_t = (1 - 2B + B^2) X_t = X_t - 2 X_{t-1} + X_{t-2} \ \)

Jeżeli operatow \(\nabla \ \) zaaplikujemy do funkcji liniowej

\(m_t = at + b \ \)

wynikiem będzie funkcja stała

\(\nabla m_t = a. \ \)

W ten sposób dowolny wielomian rzędu \( k \) możemy zredukować do stałej poprzez k-krotne użycie operatora \( \nabla \ \) (\( \nabla^k \ \)).

Jeżeli rozpatrujemy model

\( X_t = m_t + Y_t \ \)

gdzie

\( m_t = \sum_{i=0}^{k} a_j t^j, \)

a \( Y_t \ \) jest procesem stacjonarnym o średniej zero, to

\( D_t = \nabla^k X_t = k! a_k + \nabla^k Y_t \)

jest procesem stacjonarnym o średniej

\( \langle D_t \rangle = k! a_k. \)

Możemy zatem, mając dowolny proces \( \{ x_t \} \ \) w którym trend dany jest wielomianem niewiadomego rzędu zróżnicować go, powiedzmy m-razy, aż znajdziemy taki zbiór danych \( \{ \nabla^m x_t \} \), który będzie realizacją stochastycznego procesu stacjonarnego.

W rzeczywistych realizacjach tej metody zwyczajowo wystarczy 1-2 krotne różnicowanie szeregu.

Uwaga

Jeżeli amplituda danych otrzymanych metodą 3 wykazuje wzrost wraz z wzrostem t, to można temu zapobiec logarytmując na początku dane wejściowe

\( \{ y_t = \ln (x_t) \} \ \)

a dopiero później różnicując potomny proces \( \{ y_t \} \ \).

function [z] = asc_diff(x, k = 2)
 
  if ( isempty (x) || (!isvector(x)) ) 
    error ("asc_diff: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if ( k < 1 )
    error("k > 0, podałeś k = %d", k);
    return
  endif
 
  for i = 1 : k
    z = x(2 : length(x)) - x(1 : length(x)-1);
    x = z;
  endfor
 
endfunction

Eliminacja Trendu oraz Sezonowości

przykładowe dane wykazujące sezonowość oraz wykładniczy wzrost (trend)

Poprzednie metody w naturalny sposób dają się rozszerzyć i mogą zostać zaadaptowane do szeregów czasowych wykazujących zarówno trend jak i sezonowość

\( X_t = m_t + s_t + Y_t \ \)

gdzie

\( \begin{align} (i)~~ &E Y_t = 0,\\ (ii)~~ &s_{t+d} = s_t, \ \\ (iii)~~ &\sum_{j=1}^d S_j = 0. \end{align} \)

Metoda S1: Metoda małego trendu

Rys.MS1.1: sezonowy szereg czasowy o małym trendzie

Rys.MS1.2: trend \( \hat{m}_j \)

Rys.MS1.3: czynnik okresowy \( \hat{s}_k \)

Rys.MS1.4: reszty modelu \( Y_{jk} \)

Zgodnie z powyższymi założeniami mamy

\( \sum_{j=1}^d S_j = 0. ~~~ \ \) (MS1.1)

Dodatkowo, jeżeli o naszych danych wejściowych możemy powiedzieć, że dla danego okresu d trend jest stały, to w j-tym okresie

\( \hat{m}_j = \frac{1}{s} \sum_{k=1}^s X_{jk}, \ \) (MS1.2)

czyli

\( \hat{s}_k = \frac{1}{N} \sum_{j=1}^N (X_{jk} - \hat{m}_j), \ \), N - ilość okresów. (MS1.3)

Powyższy wzór spełnia zależność (M1). W rezultacie obliczamy wektor reszt.

\( Y_{jk} = X_{jk} - \hat{m}_j - \hat{s}_k,\)

\( ~~~~~j = 1,\dots,N, k=1,\dots,d\ \). (MS1.4)

\( Y_{jk} \ \) to reszta (błąd) dla k-tej wartości w <j>-tym okresie.

Przykład 

Rozpatrzymy przykładowy (sztuczny) szereg czasowy wykazujący zarówno sezonowość, jak i trend (Rys.MS1.1). Trend, zgodnie z założeniami powinien być dostatecznie mały, aby na danym okresie mógł być uznany za stały. Aby zilustrować działanie metody 1, obliczymy kolejno \( \hat{m}_j \) (Rys.MS1.2), \( \hat{s}_k \) (Rys.MS1.3) i ostatecznie oszacujemy błąd (Rys.MS1.4) korzystając z formuły (MS1.4). Poniżej czytelnik znajdzie propozycję skryptu (m-plik) zawierającego funkcję pozwalającą oszacować wszystkie powyższe wielkości.

function [m_hat, s_hat, Y] = asc_small_trend(x, d)
 
  if ( isempty (x) || (!isvector(x)) ) 
    error ("asc_diff: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if ( d < 1 )
    error("d > 0, podałeś d = %d", d);
    return;
  endif  
 
  ilosc_okresow = floor(length(x) / d);
 
  %%%
  % średnia na okres
  m_hat = zeros(1,ilosc_okresow);
  for j = 1 : ilosc_okresow
    for k = 1 : d
      m_hat(j) += x((j-1)*d + k);
    endfor
  endfor
  m_hat /= d;
 
  %%%
  % sezonowość
  s_hat = zeros(1,d);
  for k = 1 : d
    for j = 1 : ilosc_okresow
       s_hat(k) += x((j-1)*d + k) - m_hat(j);
    endfor
  endfor
  s_hat /= ilosc_okresow;
 
  %%%
  % wektor reszt
  for j = 1 : ilosc_okresow
    for k = 1 : d
      Y((j-1)*d + k) = x((j-1)*d + k) - m_hat(j) - s_hat(k);
    endfor
  endfor
 
endfunction

Metoda S2: Szacowanie metodą średniej kroczącej

Rys.MS2.1: sezonowy szereg czasowy o silnym (eksponencjalnym) trendzie

Rys.MS2.2: wstępny trend \( \hat{m}_t \)

Rys.MS2.3: czynnik okresowy \( \hat{s}_k \)

Technika ta ma lekką przewagą nad swoją poprzedniczką, bowiem nie musimy zakładać, że średnia liczona na okresie zmian jest praktycznie stała.

Przepis

1. Zakładamy, że mamy wektor obserwacji \( \{ X_t \}, ~~~t=1,\dots,N \ \).
2. Na początku szacujemy trend poprzez zastosowanie specjalnego filtra, tak, by wygładzić część odpowiedzialną za sezonowość i by zmniejszyć szum.
   1. Jeżeli okres jest liczbą parzystą (d = 2q), to stosujemy następujący filtr \( \hat{m}_t = (\frac{1}{2} X_{t-q} + X_{t-q+1} + \ \) \( \dots + X_{t+q-1} + X_{t+q})/d, ~~~q < t \le N-q. \ \)
   2. Jeżeli okres jest liczbą nieparzystą, to stosujemy zwykłą średnią kroczącą.
3. Kolejny krok to próba oszacowania sezonowości. Dla każdej wartości \( k = 1, \dots, d \) obliczamy średnią odchyleń

\( w_k = \{ (x_{k + jd} - \hat{m}_{k + jd}): q < k + jd \le N-q \}. \)

Suma właśnie obliczonych odchyleń niekoniecznie daje zero. Komponent sezonowy lepiej jest zatem liczyć ze wzoru

\( \hat{s}_k = w_k - \frac{1}{d} \sum_{i=1}^{d} w_i, ~~~k=1,\dots,d \ \),

mając dodatkowo na uwadze, że \( \hat{s}_k = \hat{s}_{k-d}, ~~~k>d. \ \) Dane pozbawione sezonowości można zapisać jako

\( d_t = X_t - \hat{s}_t \ \).

4. Na sam koniec szacujemy trend posługując się najbardziej odpowiednią metodą, np. filtrem MA, czy też próbując dopasować wielomian do danych \( \{ d_t \} \ \). Dzięki temu otrzymujemy ostatecznie wektor trendu \( \hat{m}_t \ \).
5. Ostatecznie obliczamy wektor reszt, standardowo odejmując komponenty sezonowości i trendu od danych wejściowych

\( Y_t = X_t - \hat{m}_t - \hat{s}_t, ~~~t=1,\dots,N. \ \)

Przykład 

Rozpatrzymy inny (sztuczny) szereg czasowy wykazujący zarówno sezonowość, jak i trend (Rys.MS2.1). Trend tym razem może być dowolnie duży, tu zastosujemy trend eksponencjalny. Aby zilustrować działanie metody 1, obliczymy kolejno wstępny trend (metodą odpowiedniej średniej kroczącej) \( \hat{m}_j \) (Rys.MS2.2), średni okres \( \hat{s}_k \) (Rys.MS2.3), prawdziwy trend \( d_t \) (Rys.MS2.4) i ostatecznie oszacujemy błąd (Rys.MS1.5). Poniżej czytelnik znajdzie propozycję skryptu (m-plik) zawierającego funkcję pozwalającą oszacować wszystkie powyższe wielkości.

Rys.MS2.4: dane bez sezonowości \( d_t \).

Rys.MS2.5: reszty modelu \( Y_t \).

function [m_hat, s_hat, d, Y] = asc_seasonal_moving_average(x, d = -1, stopien_MA = 17)
 
  if ( isempty (x) || (!isvector(x)) ) 
    error ("asc_seasonalMA: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if ( d == -1 )
    error("Jeszcze nie zakodowane...");
  elseif (d == 0)
    error("asc_seasonalMA: d > 0, podałeś d = %d", d);
    return;
  endif  
 
  length_x = length(x);
  ilosc_okresow = floor(length_x / d);
 
  %%%
  % wstępne wygładzanie danych
  q = floor(d / 2);
  switch (mod(d, 2))
    case 0  % even
      for i = q + 1 : length_x - q
        m_hat(i - q) = sum(x(i - q:i + q)) - 0.5 * (x(i - q) + x(i + q));
      endfor
      m_hat /= d;
    case 1  % odd
      m_hat = asc_moving_average(x, d);
  endswitch
 
  %%%
  % obliczanie sezonowości
  % (1) średnia odchyleń
  for k = 1 : d
    w(k) = mean(x(q + k: d :length_x - q - 1) - m_hat(1 + k: d :length(m_hat)));
  endfor
 
  % (2) średni okres zmian
  s_hat = w - mean(w);
 
  % (3) sezonowość
  s_hat_long = repmat(s_hat, 1, ilosc_okresow - 1);  
 
  %%%
  % dane bez sezonowości
  d = x(q + 1:length_x - q - 1) - s_hat_long;
 
  %%%
  % reszty 
  % (trend obliczany jest prostą średnią kroczącą)
  floor_stopien_MA_pol = floor(stopien_MA / 2);
  Y = x(q + 1 + floor_stopien_MA_pol:length_x - q - 1 - floor_stopien_MA_pol) 
    - asc_moving_average(d, stopien_MA) 
    - s_hat_long(1 + floor_stopien_MA_pol:length(s_hat_long) 
    - floor_stopien_MA_pol);
 
endfunction

Metoda S3: Różnicowanie z przesunięciem d

Rys.MS3.1: sezonowy szereg czasowy o małym trendzie

Rys.MS3.2: trend \( \hat{m}_t - \hat{m}_{t-d}\)

Rys.MS3.3: reszty modelu (po różnicowaniu) \( Y_t - Y_{t-d}\)

Metodę, którą poznaliśmy nieco wcześniej (różnicowanie do danych stacjonarnych) możemy łatwo dostosować do danych zmiennych sezonowo z okresem zmienności d poprzez wprowadzenie operatora różnicowego z opóźnieniem d, zdefiniowanego jako

\( \nabla_d X_t = X_t - X_{t-d} = (1 - B^d) X_t. \ \)

Uwaga 

nie powinno się mylić operatora \( \nabla_d = (1 - B^d) \ \) z operatorem \( \nabla^d = (1 - B)^d. \ \)

Stosując operator \( \nabla_d \ \) do danych wykazujących sezonowość (i trend) \( \{ X_t \} \ \)

\( X_t = m_t + s_t + Y_t, \ \)

gdzie \( \{ s_t \} \ \) jest komponentem periodyczną funkcję o okresie d dostajemy

\( \nabla_d X_t = m_t - m_{t-d} + Y_t - Y_{t-d}. \ \)

Otrzymujemy w ten sposób nowy szereg czasowy składający się z dwóch komponentów:

• trendu \( m_t - m_{t-d} \ \) i
• szumu (reszt) \( Y_t - Y_{t-d}. \ \)

Aby obliczyć wektor reszt, trend możemy usunąć dowolną metodą zaprezentowaną wyżej, np: używając operatora \( \nabla^k.\ \)

Przykład 

Rozpatrzymy kolejny (sztuczny) szereg czasowy wykazujący zarówno sezonowość, jak i trend (Rys.MS3.1). Trend i tym razem może być dowolnie duży. Aby zilustrować działanie metody S3, usuniemy sezonowość (korzystając z operatora różnicowania z przesunięciem) (Rys.MS3.2) oraz w końcu ostatecznie oszacujemy błąd (podobną metodą - różnicowania) (Rys.MS3.3). Poniżej czytelnik znajdzie propozycję skryptu (m-plik) zawierającego funkcję pozwalającą oszacować wszystkie powyższe wielkości.

function [m, Y] = asc_diffd(x, d, accuracy = 0.001, maxdeg = 5)
 
  if ( isempty (x) || (!isvector(x)) ) 
    error ("asc_diffd: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  if ( d == -1 )
    error("asc_diffd: Jeszcze nie zakodowane...");
  elseif (d == 0)
    error("asc_diffd: d > 0, podałeś d = %d", d);
    return;
  endif  
 
  %%%
  % wektor trendu
  m = x(d : length(x)) - x(1 : length(x)-d+1);
 
  %%%
  % wektor reszt - obliczany przez różnicowanie, do momentu, gdy
  % osiągniemy żądaną dokładność na zerową średnią reszt
  for deg = 1 : maxdeg
    Y = asc_diff(m,deg);
    if mean(Y) < accuracy
      break;
    endif
  endfor 
 
endfunction

Metoda S4: Wygładzanie wykładnicze 3 parametrową metodą Wintersa

Niekiedy dane, które poddajemy dekompozycji posiadają sezonowość multiplikatywną, innym razem addytywną. Właśnie omawiana metoda posiada dwa odrębne mechanizmy dekompozycji sygnału w zależności od typu sezonowości.

Sezonowość multiplikatywna

Zakładamy, że szereg czasowy może być reprezentowany modelem

\( X_t = (b_1 + b_2 t) s_t + Y_t \ \)

gdzie

• \(b_1\) to podstawowy sygnał (komponent stały),
• \(b_2\) to trend,
• \(s_t\) odpowiada za sezonowość (multiplikatywną),
• \(Y-t\) to wektor reszt (szum, wektor błędów).

Jeden sezon ma długość L okresów. Współczynniki sezonowości są zdefiniowane tak, że sumują się do długości sezonu L

\( \sum_{1 \le t \le L} s_t = L. \ \)

Model działa również bardzo dobrze jeżeli pozbędziemy się wcześniej trendu \(b_2\).

Stosowalność 

do szeregów czasowych, w których amplituda sezonowości jest proporcjonalna do średniego poziomu szeregu - do szeregów wykazujących sezonowość multiplikatywną.

Notacja 

• \(R_T \ \) chwilowa wielkość procesu z usuniętą sezonowością pod koniec okresu T,
• \(\hat{R}_t\) oszacowanie szeregu pozbawionego sezonowości,
• \(\hat{G}_t\) oszacowanie trendu,
• \(\hat{S}_t\) oszacowanie komponentu okresowego (sezonowości).

Procedura 

Krok 1. wygładzanie ogólne

\( \hat{R}_t = \alpha \frac{X_t}{\hat{S}_{t-L}} + (1 - \alpha) (\hat{R}_{t-1} + \hat{G}_{t-1}), \)

gdzie \( \alpha \in (0,1) \) to stała wygładzająca. Dzieląc \(X_t\) przez czynnik sezonowy obliczony dla poprzedniego okresu zmian (\(\hat{S}_{t-L}\)) powoduje usunięcie sezonowości z danych - pozostawiamy tylko trend i czynnik stały.

Krok 2. wygładzanie trendu

\( \hat{G}_t = \beta (\hat{S}_{t} - \hat{S}_{t-1}) + (1 - \beta) \hat{G}_{t-1}, \)

gdzie \( \beta \in (0,1) \) to druga stała wygładzająca. Proste wygładzanie trendu z danych pozbawionych sezonowości.

Krok 3. wygładzanie komponentu okresowego

\( \hat{S}_t = \gamma \frac{X_t}{\hat{S}_{t}} + (1 - \gamma) \hat{S}_{t-L}, \)

gdzie \( \gamma \in (0,1) \) to druga stała wygładzająca.

Wartości początkowe dla algorytmu 

aby algorytm działał dobrze wymagana jest znajomość pierwszych dwóch okresów danych (2L). Jeżeli znamy m okresów, oraz niech \(\hat{X}_j\), j = 1,2,...,mL oznacza średnią z j-tego sezonu, to

1. trend

\( \hat{G}_1 = \frac{\hat{X}_m - \hat{X}_1}{(m-1)L}, \)

2. zmienna pozbawiona sezonowości

\( \hat{R}_1 = \hat{X}_1 - \frac{L}{2} \hat{G}_1, \)

3. komponent sezonowości

obliczany dla każdego kroku czasowego \(t = 1, 2, \dots, mL\)

\( \hat{S}_t = \frac{\hat{X}_t}{\hat{X}_i - [(L+1)/2 - j]\hat{G}_1}, \)

gdzie \(\hat{X}_i\) oznacza średnią w okresie odpowiadającym t, a j to pozycja okresu t w sezonie. Powyższe równanie daje m oszacowań komponentu sezonowego na każdy okres

\( \hat{S}_t = \frac{1}{m} \sum_{k=0}^{m-1} \hat{S}_{t+kL}, ~~~t = 1,2,\dots,L,\)

\( \hat{S}_t(0) = \hat{S}_t \frac{L}{\sum_{t=1}^{L} \hat{S}_{t}}, ~~~t = 1,2,\dots,L.\)

Sezonowość addytywna

Zakładamy, że szereg czasowy może być reprezentowany modelem

\( X_t = b_1 + b_2 t + s_t + Y_t \ \)

gdzie

• \(b_1\) to podstawowy sygnał (komponent stały),
• \(b_2\) to trend,
• \(s_t\) odpowiada za sezonowość (addytywną),
• \(Y-t\) to wektor reszt (szum, wektor błędów).

Jeden sezon ma długość L okresów. Współczynniki sezonowości są zdefiniowane tak, że sumują się do długości sezonu L

\( \sum_{1 \le t \le L} s_t = 0. \ \)

Model działa również bardzo dobrze jeżeli pozbędziemy się wcześniej trendu \(b_2\).

Stosowalność 

do szeregów czasowych, w których amplituda sezonowości jest niezależna od średniego poziomu szeregu - do szeregów wykazujących sezonowość addytywną.

Notacja 

• \(R_T \ \) chwilowa wielkość procesu z usuniętą sezonowością pod koniec okresu T,
• \(\hat{R}_t\) oszacowanie szeregu pozbawionego sezonowości (pod koniec okresu t),
• \(\hat{G}_t\) oszacowanie trendu (pod koniec okresu t),
• \(\hat{S}_t\) oszacowanie komponentu okresowego (sezonowości, pod koniec okresu t).

Procedura 

Krok 1. wygładzanie ogólne

\( \hat{R}_t = \alpha (X_t - \hat{S}_{t-L}) + (1 - \alpha) (\hat{R}_{t-1} + \hat{G}_{t-1}), \)

gdzie \( \alpha \in (0,1) \) to stała wygładzająca.

Krok 2. wygładzanie trendu

\( \hat{G}_t = \beta (\hat{S}_{t} - \hat{S}_{t-1}) + (1 - \beta) \hat{G}_{t-1}, \)

gdzie \( \beta \in (0,1) \) to druga stała wygładzająca. Proste wygładzanie trendu z danych pozbawionych sezonowości (tak jak poprzednio).

Krok 3. wygładzanie komponentu okresowego

\( \hat{S}_t = \gamma (X_t-\hat{S}_{t}) + (1 - \gamma) \hat{S}_{t-L}, \)

gdzie \( \gamma \in (0,1) \) to druga stała wygładzająca.

Wartości początkowe dla algorytmu 

jak wyżej.

Uwaga 1

wszystkie powyższe metody możemy stosować wymiennie i łącznie, nie musimy ograniczać się do jednej wybranej na samym początku naszej analizy.

Wszystkie prezentowane kody źródłowe do programów w języku Matlab / Octave jeżeli nie 
sprecyzowano inaczej dostępne są na licencji LGPL v3: ASC