Z Skrypty dla studentów Ekonofizyki UPGOW
Analiza Szeregów Czasowych <<< Techniki analizy szeregów czasowych
Spis treści |
Matlab / GNU Octave
Matlab oraz jego otwarty odpowiednik GNU Octave stanowią środowisko programistyczne oraz język interpretowany w jednym. Nazwa Matlab pochodzi od MATrix LABoratory i w rzeczywistości największa siła owego narzędzia tkwi w obliczeniach bazowanych na macierzach (wektoryzacji obliczeń). Po więcej informacji ogólnych o obu system odsyłamy czytelnika do wikipedii. Poniżej prezentujemy tez kilka bardzo użytecznych linków mogących służyć do podstawowego i nieco bardziej zaawansowanego opanowania Matlab / GNU Octave.
Odnośniki
- na tej platformie
- internet
Wizualizacja danych
Wizualizacja danych to podstawowe i bardzo potężne narzędzie analizy szeregów czasowych. Ciąg liczb nigdy nikomu nic nie powie jeżeli nie zostanie zaprezentowany na wykresie. Ten krótki podrozdział omówi podstawowe formatowanie plików wyjściowych. Wstęp do rysowania wykresów zawarty jest w tym rozdziale o grafice w podręczniku Programowanie. Tutaj nieco rozszerzymy tą wiedzę.
Zapoznaj się z komendami: figure, close, hold, clf.
Podstawową instrukcją do rysowania wykresów jest komenda plot. Przykładowe komendy
x = 0:2*pi/100:2*pi; y = power(sin(x),2); plot(x,y)
spowodują narysowanie linii składającej się ze 100 punktów a przedstawiającą kwadrat funkcji trygonometrycznej sin na przedziale \([0, 2 \Pi] \) w aktywnym oknie graficznym. Następne komendy posłużą do podstawowego opisu osi wykresu i tytułu. Proszę zauważyć, że Octave rozumie formatowanie znane z latexa, takie jak podniesienie (sub: ^), indeksowanie (sup: _) czy symbole matematyczne (\(\sqrt{\cdot}, \gamma, \dots\)).
xlabel('x', 'FontName', 'Times', 'FontSize', 18) ylabel('sin(x)^2', 'FontName', 'Times', 'FontSize', 18) title('Funkcje trygonometryczne', 'FontName', 'Times', 'FontSize', 22)
Ustawiliśmy też czcionkę, jaką chcemy wykorzystać i jej rozmiar (dobrze jest jeżeli czcionka istnieje w systemie ;) ). To teraz lekko podrasujemy wygląd naszego wykresu
plot(x, y, '--rs;sin^2(x);', 'LineWidth', 2, 'MarkerSize', 10) grid on legend('Location', 'North')
Po kolei:
- '--' to linia przerywana
- 'r' to kolor wykresu (czerwony)
- 's' to znaczniki (kwadraty) pojawiające się w miejscu punktów
- ';sin^2(x);' to legenda wykresu
- 'LineWidth' i następujący po nim numer ustawiają szerokość linii
- 'MarkerSize' i następujący po nim numer ustawiają wielkość znaczników
- grid on - włącza siatkę (off wyłącza)
- legend('Location', 'North') - ustawia legendę na środku u góry wykresu
Czasami potrzebne jest formatowanie napisów tekstowych w np: legendach, po ty by przykładowo wstawić automatycznie obliczoną średnią czy wartość jakiejś zmiennej do ciągu znaków. Można to zrobić wykorzystując np: funkcję sprintf. Działa ona identycznie jak z języku c.
plot(x, y, sprintf('--rs;{/Symbol m} = %.2f;', mean(y)), 'LineWidth', 2, 'MarkerSize', 10)
Przy okazji pokazaliśmy jak można uzyskać symbole greckich liter w ciągach znaków i w ogólności jak lokalnie zmodyfikować czcionkę (font).
Zapisywanie grafiki w pliku
Funkcja print umożliwia zapisywanie do wielu formatów graficznych zarówno rastrowych (png, jpg, gif, ...) jak i wektorowych (cdr, svg, eps, ....). Przykładowo, poniższa komenda zapisze nam wynik naszej pracy do pliku png, z rozmiarem 640x480,
print('wykres.png', '-dpng', '-S640,480')
- Najważniejsze formaty
- -dformat
- ps, ps2, psc, psc2 - Postscript (level 1 i 2, mono i color)
- eps, eps2, epsc, epsc2 - Encapsulated postscript (level 1 i 2, mono i color)
- tex, epslatex, epslatexstandalone, pstex, pslatex - generuje plikiu LaTeX (TeX) i odpowiednie pliki EPS (plik wyprodukowany przez epslatexstandalone może być przetwarzany bezpośrednio przez LaTex).
- ill, aifm - Adobe Illustrator
- cdr, corel - CorelDraw
- dxf - AutoCAD
- png - Portable Network Graphics
- jpg, jpeg - JPEG
- gif - GIF
- svg - Scalable Vector Graphics
- pdf - Portable Document Format
- inne przełączniki
- -Sxsize,ysize - ustawia rozmiar obrazka w pikeslach (PNG i SVG)
- -Ffontname:size - ustawia czcionkę i (opcjonalnie) jej rozmiar
Obiekty graficzne
Odpowiednie sformatowanie obiektu graficznego pozwoli nam stworzyć czytelny rysunek. Jeżeli np: zostawimy wszystkie ustawienia jako podstawowe Octave wydrukuje nam plik wielkości Letter (A4), nie dostosowując automatycznie czcionek do wielkości rysunku, rozdzielczości, aspektu (wysokości do szerokości) czy wreszcie ułożenia marginesów. Wszystkie te wielkości możemy dostosować samodzielnie. Co więcej powinniśmy to robić!
Do podstawowej komendy obsługującej format wykresu należy uchwyt gcf. Komenda
h = gcf();
zwraca uchwyt do obecnie budowanego obiektu typu figure. Jeżeli takowy nie istnieje to go tworzy. Możemy teraz na tym obiekcie dokonywać zmian (set) lub uzyskiwać o nim informacje (get). Kod
h = gcf(); set (h, 'papertype', 'a4'); set (h, 'paperunits', 'centimeters'); set (h, 'papersize', [4 3]); set (h, 'paperposition', [0,0,[4 3]]); set (h, 'defaultaxesposition', [0.15, 0.15, 0.75, 0.75]);
definiuje kilka podstawowych własności graficznego okna. Po dokładny opis możliwych modyfikacji odsyłam do literatury.
Podobną funkcją jest funkcja gca zwracająca uchwyt do obsługi osi na aktualnie modyfikowanym obiekcie. Np: ustawienie
ax = gca(); set (ax, 'xlim', [0, 2*pi]);
spowoduje ustawienie zakresu zmienności osi x na [-10:10] niezależnie od danych na wykresie. Natomiast ustawienie
a = gca; set(a,'XTick', 0 : pi / 2 : 2 * pi) set(a, 'XTickLabel', {'0', '{/Symbol p}/2', '{/Symbol p}', '3{/Symbol p}/2', '2{/Symbol p}'} )
spowoduje wygenerowanie poprawnego opisu osi dla naszego przykładu.
Wszystkie powyższe rysunki są wielkości 640x480px i pomniejszone poprzez skalowanie HTML w tym dokumencie. Łatwo zauważyć nieostre linie i czcionki. Aby uniknąć takiej sytuacji można generować (poprze set) rysunki o dokładnie takiej wielkości jaką chcemy użyć na stronie czy w innej dowolnej publikacji. Jeżeli chcemy drukować to większe rysunki (z zachowaniem proporcji!!!) będą bardziej odpowiednie (np: 1200x800px, 3600x2400px). Poniżej znajdziecie odpowiednio sformatowany rysunek, tak, że wszystkie linie są ostre i rysunek dalej pozostaje czytelny (czcionki są odpowiednio duże i czytelne). Poniższy rysunek ma 320x240 px.
Najczęściej używane własności
Style linii
- linestyle
- styl linii może być jedną z
- "-" linia ciągła,
- "--" linia przerywana,
- ":" linia wykropkowana,
- "-." kropka - kreska,
- linewidth
- określa szerokość linii, podstawowa wartość to 1, 2 oznacza dwukrotnie grubszą.
plot(..., 'Linewidth', 2, ...)
Znaczniki
- marker
- może przyjmować wartości
- '+' - krzyż
- '*' - gwiazdka
- 'o' - otwarte kółko
- 'x' - iks
- '^' - trójkąt
- 's' - pełny kwadrat
- 'p' - pusty kwadrat
- 'd' - diament
- 'h' - diament pusty
- MarkerSize
- definiuje wielkość znaczników w punktach, np:
plot(..., 'MarkerSize', 14, ...)
Kolory
Można specyfikowac jako tryplety RGB, lub użyć jednego z predefiniowanych (numerów, skrótów lub pełnych nazw):
- 'r', 'red' - czerwony
- 'g', 'green' - zielony
- 'b', 'blue' - niebieski
- 'm', 'magenta' - magenta
- 'c', 'cyan' - cyan
- 'w', 'white' - biały
- 'b', 'black' - czarny
- 'y', 'yellow' - żółty
Pakiet ASC
Pakiet Analiza Szeregów Czasowych (ASC) składa się obecnie z 14 plików - funkcji w języku MATLAB / GNU Octave. Podrozdział ten będzie prawdopodobnie najszybciej i najbardziej intensywnie rozwijanym rozdziałem tego skryptu.
Opis funkcji
- asc_armafit
Funkcja [phi, theta, v] = asc_armafit(x, p, q) dopasowywuje proces ARMA(p,q) do danego szeregu czasowego x przyjmuje jako agrumenty x - stacjonarny szereg czasowy (niepusty wektor danych) p - rząd wielomianu AR, p >= 0 q - rząd wielomianu MA, q >= 0 Ten program wymaga dowolnego programu obliczającego wstępne wartości 'theta' np: algorytmu innowacyjnego zaimplementowanego w asc_innovation(x,q) Zwraca phi - wektor parametrów AR(p) theta - wektor parametrów MA(q) v - wariancję procesu ARMA(p,q)
użycie
octave:44> x = arma_rnd([0.9,-0.3], [0.63], 1, 50000, 10000); octave:45> [phi, theta, v] = asc_armafit(x, 2, 1) armafit(x,p,q) p = 2 q = 1 -------------- phi = 0.17814 0.73459 theta = 0.41018 v = 3.7739
- asc_diff
Funkcja [Y] = asc_diff(x, k = 2) oblicza wektor reszt szeregu czasowego stosując metodę różnicowania do danych stacjonarnych Przyjmuje jako argumenty x - szereg czasowy (niepusty wektor danych) k - rząd różnicowania, jeżeli nie podany to k = 2 Zwraca Y - wektor reszt
- asc_diffd
Funkcja [m, Y] = asc_diffd(x [, d = -1, accuracy = 0.001, maxdeg = 5]) oblicza wektor reszt szeregu czasowego wykazujacego sezonowość stosując metodę różnicowania z przesunięciem d do danych stacjonarnych Przyjmuje jako agrumenty x - szereg czasowy (niepusty wektor danych) d - opcjonalnie: okresowość danych, jezeli nie podane program sam stara się oszacować okresowość !!! jeszcze nie działa !!! accuracy - opcjonalnie; dokładność średniej wektora reszt, decyduje o stopniu różnicowania szeregu z usuniętą sezonowością [= 0.001] maxdeg - opcjonalnie; maksymalna wartość stopnia różnicowania szeregu z usuniętą sezonowością [= 5] Zwraca m - wektor danych po usunięciu sezonowości Y - wektor reszt
- asc_durbin_levinson
Funkcja [phi, v] = asc_durbin_levinson(x, p) dopasowywuje proces AR(p) do danego szeregu czasowego metodą Durbina - Levinsona Przyjmuje jako argumenty x - stacjonarny szereg czasowy (niepusty wektor danych) p - rząd wielomianu AR, p >= 0 Zwraca phi - wektor parametrów procesu AR(p) v - wariancję procesu AR(p)
- asc_exponential_smoothing
Funkcja [m] = asc_exponential_smoothing(x, a = 0.5) wygładza szereg czasowy Przyjmuje jako agrumenty x - szereg czasowy (niepusty wektor danych) a - parametr wygładzania wykładniczego 0 < a < 1 [= 0.5] Zwraca m - wygładzony wektor procesu {x}
- asc_exponential_smoothing_Holt
Funkcja [F, S] = asc_exponential_smoothing_Holt(x, a = 0.5, b = 0.5) wygładza szereg czasowy przyjmuje jako agrumenty x - szereg czasowy (niepusty wektor danych) a - parametr wygładzania wykładniczego 0 < a < 1 [= 0.5] b - parametr wygładzania wykładniczego 0 < b < 1 [= 0.5] initial - wybór parametrów początkowych dla F i S vec: F(1) = x(1), S(1) = x(2) - x(1) lin : F(1) - współczynnik wolny liniowego dopasowania do x S(1) - współczynnik kierunkowy liniowego dopasowania do x Zwraca F - wygładzona wartość zmiennej {x} S - wygładzona wartość trendu
- asc_fakc
Funkcja [z] = asc_fakc(x, stopien = 1) oblicza funkcję autokorelacji częściowej próby procesu stacjonarnego X, na podstawie realizacji {x} Przyjmuje jako argumenty x - szereg czasowy (niepusty wektor danych) stopien - stopień korelacji [= 1] Zwraca z - wektor autokorelacji częściowej próby, \alpha(h)
- asc_innovation
Funkcja [theta,v] = asc_innovation(x, q) funkcja dopasowywuje proces MA(q) do procesu stacjonarnego X, na podstawie realizacji {x} korzystając z algorytmu innowacyjnego Przyjmuje jako argumenty x - stacjnarny szereg czasowy (niepusty wektor danych) q - rząd wielomianu MA, q > 0 Zwraca theta - wektor współczynników modelu MA(q) v - wariancję procesu MA(q)
- asc_kor
Funkcja [z] = asc_kor(x, stopien = 0) oblicza funkcję autokorelacji próby procesu stacjonarnego X, na podstawie realizacji {x} przyjmuje jako argumenty x - szereg czasowy (niepusty wektor danych) stopien - stopień korelacji [= 0] Zwraca z - wektor autokorelacji \rho(h)
- asc_kow
Funkcja [z] = asc_kow(x, stopien = 0) oblicza funkcję autokowariancji próby procesu stacjonarnego X, na podstawie realizacji {x} przyjmuje jako argumenty x - szereg czasowy (niepusty wektor danych) stopien - stopień korelacji [= 0] Zwraca z - wektor autokowariancji \gamma(h)
- asc_moving_average
Funkcja [z] = asc_moving_average(x, c = 7) oblicza c-punktową średnią kroczącą centralną Przyjmuje jako argumenty x - szereg czasowy (niepusty wektor danych) c - ilość punktów do obliczenia średniej c < length(x) [= 7] c musi być nieparzyste Zwraca z - krzywą dopasowaną do danych {x} uwaga: length(z) = length(x) - c + 1
- asc_seasonal_moving_average
Funkcja [m_hat, s_hat, d, Y] = asc_seasonal_moving_average(x [, d = -1, stopien_MA = 17]) szacuje waktory trendu, sezonowości oraz reszt danych wykazujących sezonowość metodą średniej kroczącej x - wektor danych (szereg czasowy) d - opcjonalnie: okresowość danych, jezeli nie podane program sam stara się oszacować okresowość [= -1] stopien_MA - opcjonalnie, ilość punktów dla średniej kroczącej [= 17] zwraca m_hat - wstępny wektor trendu s_hat - średnią okresowość d - wektor trendu Y - wektor reszt
- asc_small_trend
Funkcja asc_small_trend(x, d) szacuje wektory trendu, sezonowości oraz reszt danych wykazujących sezonowość metodą małego trendu x - wektor danych (szereg czasowy) d - okresowość danych Zwraca m_hat - wektor średniej z wartości szeregu {x} w danym okresie length(m_hat) = length(x) / d s_hat - średni wektor sezonowości (jeden okres) Y - wektor reszt
- asc_spencer
Funkcja: [z] = asc_spencer(x) zwraca 15-punktową średnią kroczacą w/g Spencera z wagami [74, 67, 46, 21, 3, -5, -6, -3] / 320; x - 1D wektor danych wejściowych [length(X) >= 18]
Źródła
Źródła pakietu ASC dostępne są na licencji LGPL v3 poprzez platformę GIT z poniższego adresu
Pakiet TSA
W internecie istnieje wiele pakietów oraz pojedynczych funkcji w jakiś sposób realizujących to co czytelnik znajdzie w poprzednim paragrafie. Wiele z nich jest niepełnych, wiele projektów zostało zawieszonych kilka lat temu. Wreszcie istnieją programy narzędziowe w wielu innych niż MATLAB językach realizujące dokładnie takie same zagadnienia. Pozostając jednak przy MATLAB / GNU Octave chyba obecnie najlepszym i najbardziej rozbudowanym pakietem jest pakiet Pana Aloisa Schlögla pod nazwą Time Series Analysis (TSA). Wiele funkcji znajdujących się w tym pakiecie będzie się powtarzać z pakietem ASC tworzonym na potrzeby tego wykładu ale znajdą tez Państwo wiele nowych.
W miarę potrzeb w tym miejscu pojawią się tłumaczenia i przykłady użycia bardziej pomocnych funkcji pochodzących z prezentowanego pakietu.
- Miejsce skąd można ściągnąć najnowsze źródła pakietu
Pobieranie TSA - Opis plików (funkcji) wchodzących w skład pakietu
Pomoc on-line - Użytkownicy Ubuntu mogą zainstalować pakiet TSA poleceniem
sudo apt-get install octave-tsa
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