Analiza Szeregów Czasowych/Dekompozycja szeregu czasowego

Z Skrypty dla studentów Ekonofizyki UPGOW

Sat 13 Feb 2010 14:17:49

Analiza Szeregów Czasowych

Indeks Dow Jones Industrial Average w latach 1928-1930

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:

  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".

Tabela 4.1
Matlab / Octave Python 
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.eps','-deps','-FTimes-Roman:24')
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.eps','-deps','-FTimes-Roman:24');
print('cw41reszty.png','-dpng','-FTimes-Roman:24')

TBA

Metoda 2 (wygładzanie krzywych metodą ś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 \gt n -q \ \).

Źródło „http://172.16.3.1/ekonofizyka/index.php/Analiza_Szereg%C3%B3w_Czasowych/Dekompozycja_szeregu_czasowego”