Z Skrypty dla studentów Ekonofizyki UPGOW
(→Ćwiczenie 4.1) |
(→Ćwiczenie 4.1) |
||
Linia 71: | Linia 71: | ||
close all | close all | ||
y = csvread('PopulacjaUSA1790-1980.csv'); | y = csvread('PopulacjaUSA1790-1980.csv'); | ||
- | X = y(:, | + | X = y(:,1); |
- | Y = y(:, | + | Y = y(:,2)/10^6; |
plot(X,Y,"sr; Dane;", | plot(X,Y,"sr; Dane;", |
Wersja z 21:32, 13 lut 2010
Sat 13 Feb 2010 14:17:49
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=1}^n (x_t - m_t)^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ą).
Ć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 |