Z Skrypty dla studentów Ekonofizyki UPGOW
Analiza Szeregów Czasowych <<< Procesy stochastyczne | Modelowanie szeregów czasowych >>>
Spis treści |
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:
- 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.
- 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.
- 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.\)
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. \)
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".
| Tabela 4.1 | |
|---|---|
| Matlab / Octave | |
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
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 \ \).
| Matlab / Octave |
|---|
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.
| Matlab / Octave |
|---|
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.
| Matlab / Octave |
|---|
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
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]. \ \)
| Matlab / Octave |
|---|
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)
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 \} \ \).
| Matlab / Octave |
|---|
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
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
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.
| Kod MS1 (Matlab / Octave) |
|---|
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
Technika ta ma lekką przewagą nad swoją poprzedniczką, bowiem nie musimy zakładać, że średnia liczona na okresie zmian jest praktycznie stała.
- Przepis
- Zakładamy, że mamy wektor obserwacji \( \{ X_t \}, ~~~t=1,\dots,N \ \).
- Na początku szacujemy trend poprzez zastosowanie specjalnego filtra, tak, by wygładzić część odpowiedzialną za sezonowość i by zmniejszyć szum.
- 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. \ \)
- Jeżeli okres jest liczbą nieparzystą, to stosujemy zwykłą średnią kroczącą.
- 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 \ \).
- 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 \ \).
- 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.
| Kod MS2 (Matlab / Octave) |
|---|
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
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.
| Kod MS3 (Matlab / Octave) |
|---|
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