Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 19:26, 18 gru 2010

Sat 13 Feb 2010 14:17:49

Analiza 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.\)

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.

Rysunek 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ą).

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

Ćwiczenie 4.1

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

close 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(Y-polyval(a,X));
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 \ \).

function [z] = asc_movingAverage(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

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.

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.

Przykład - 15 punktowa średnia krocząca Spencera (eliminującą wysokie 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)

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)
 
  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

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}, ~~~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. \ \)