niedziela, 11 stycznia 2015

Jak sprawdzić stabilność modelu ekonometrycznego w czasie?

Użycie klasycznego modelu ekonometrycznego typu AR-MA (Autoregressive–moving-average model) do progonozy zmian giełdowych może mieć sens tylko wtedy, gdy model spełni co najmniej 5 warunków. Po pierwsze prognozowane zmienne muszą być stacjonarne. Jeżeli warunek ten od razu nie jest spełniony, wtedy stosuje się ogólniejszy od ARMA model ARFIMA (Autoregressive fractionally integrated moving average). Model ten odnajduje "stopień niestacjonarności" (szeregu czasowego), który wykorzystuje do wygererowania zmiennej stacjonarnej, dla której używa już ARMA. Jeżeli pamięć długoterminowa nie występuje, wtedy ARFIMA sprowadza się do ARIMA (natomiast gdy proces jest stacjonarny, wtedy ARIMA sprowadza się do ARMA). ARFIMA może być przydatny w przypadku badania stóp zwrotu, ponieważ często okazują się one posiadać pamięć długoterminową, wywołując efekt wydłużonych ruchów w stosunku do ruchu Browna, a ten z kolei wywołuje coś w rodzaju 'ułamkowej stacjonarności'. Ten efekt jest usuwany przez ARFIMA. Jednakże w poprzednim artykule pokazałem, że efekt ten może być w rzeczywistości bardzo słaby - np. dla MWIG40 występuje istotna, ale bardzo słaba pamięć długoterminowa - stąd użycie standardowego modelu ARIMA nie będzie dla takich ruchów szczególnym błędem. W każdym razie dążymy do uzyskania stacjonarności.

Po drugie parametry modelu, jak i ogólny model muszą być istotne statystycznie.
Po trzecie model musi być stabilny - parametry modelu nie mogą się zmieniać w czasie.

Po czwarte muszą być spełnione odpowiednie warunki rozkładu prawdopodobieństwa składnika losowego. Po piąte trzeba sprawdzić liniowość / nieliniowość modelu. Ten ostatni warunek pomijam, natomiast czwarty zostanie wspomniany.

Należy rozróżniać stacjonarność i stabilność.

Stacjonarność rozumiana jest jako niezmienność rozkładu prawdopodobieństwa w czasie, a co najmniej niezmienność średniej danego procesu w czasie.

Z kolei stabilność oznacza stałość parametrów regresji. Np. jeśli mamy model AR(1):

x(t) = b*x(t-1) + e,
gdzie x(t) to np. stopa zwrotu w okresie t, zaś e to składnik losowy,

to warunkowa wartość oczekiwana zmiennej x(t) jest równa b*x(t-1). Jednakże wyznaczenie parametru b może być sztuczne, jeśli kształt regresji ulega zmianom w czasie. Dlatego musimy mieć pewność, że parametr b będzie stabilny. Powstaje pytanie jak sprawdzić czy parametry modelu pozostają stabilne w czasie? Najbardziej znane są 3 testy na stabilność modelu:
1. Test Chowa
2. Test CUSUM
3. QLR

Ad1) Test Chowa
W 1960 r. ekonomista G. Chow opracował test sprawdzający czy parametry regresji z dwóch różnych próbek są takie same. Najczęściej w ekonomii podejście to służy do testowania czy występuje punkt zwrotny w danym procesie. Cały okres dzielony jest wtedy na dwie podpróby, a moment czasowy rozdzielający je nazywamy punktem zwrotnym. To podejście jest stosowane w Gretlu.

Spróbujemy zrobić przykłady. Zanim jednak zaczniemy, dopowiedzmy, że Gretl automatycznie dzieli daną próbę na połowę, aby porównywać dwa różne okresy. Jest to podejście w miarę poprawne, jednak z praktycznego punktu widzenia lepszym rozwiązaniem wydaje się najpierw na oko oszacować czy w analizowanym okresie mają miejsce istotne zmiany kierunku. Np. możemy starać podzielić na okres hossy i bessy.
W artykule Czy oczekiwane stopy zwrotu w ogóle się zmieniają? analizowałem mWIG40 dla różnych częstości, aby sprawdzić hipotezę stacjonarności. Hipoteza ta dla wszystkich częstości została potwierdzona przez testy ADF i KPSS. Te same dane wykorzystamy teraz do sprawdzenia stabilności modelu. Zaczniemy od miesięcznych stóp zwrotu mWIG40 od początku notowań styczeń 1998-pażdziernik 2014; 198 obserwacji. Poniżej jest wykres notowań mWIG40:


Stopy zwrotu:

Ustalimy hipotetyczny punkt zwrotny w czerwcu 2007 r., gdzie miesięcznie zanotowany został szczyt historyczny. Aby dokonać testu, musimy najpierw wejść w Gretlu w Model -> Klasyczna Metoda Najmniejszych Kwadratów. Tutaj stworzymy klasyczną autoregresję. Pytanie jednak brzmi jakiego rzędu? AR(1), AR(2), a może więcej? Z pomocą przychodzi metodologia Boxa-Jenkinsa, zgodnie z którą teoretyczny model AR(p) ma funkcję autokorelacji cząstkowych (PACF) równą zeru dla opóźnień większych od p (tzn. współczynniki autokorelacji cząstkowej dla opóźnień większych niż p są nieistotne) [1].

Dotychczas nigdy nie wspomniałem o PACF, a jedynie o ACF. Czytelnik, który czytał Czy oczekiwane stopy zwrotu w ogóle się zmieniają? , zapamiętał, że zwróciłem uwagę na to, iż ACF nie pokazuje tak naprawdę autokorelacji, jeśli proces jest niestacjonarny. Jest on właśnie przydatny do naocznego oszacowania czy mamy do czynienia z procesem stacjonarnym poprzez analizę współczynników autokorelacji kolejnych rzędów. Jeżeli są one istotnie dodatnie i spadają z każdym opóźnieniem powoli, to znaczy, że prawdopodobnie mamy proces niestacjonarny (aczkolwiek w artykule Ułamkowa stacjonarność? Co to za twór? wykazałem, że metoda ta może przynieść błędne rezultaty). Z czego to wynika? ACF zawiera wewnętrzne zależności pomiędzy wartościami pośrednimi, np. autokorelacja 5 rzędu, tj. ACF(5) zawierać będzie nie tylko autokorelację pomiędzy obserwacją piątą a bieżącą, ale także autokorelacje pomiędzy obserwacjami pośrednimi, czyli np. 1, 2, 3 i 4 rzędu. Dopiero tzw. cząstkowa autokorelacja, tj. PACF, usuwa te pośrednie zależności., tak że dostaniemy poprawniejszą formę autokorelacji danego rzędu. Przy tym jednak pamiętać trzeba, że PACF jest równa ACF dla 1 rzędu, więc dla procesów niestacjonarnych również ona ma ograniczoną stosowalność. Box i Jenkins doszli do wniosku, że o ile ACF może służyć do oceny stacjonarności procesu (a przez to także determinować właściwy rząd średniej kroczącej w modelu MA), o tyle PACF może służyć za czynnik determinujący rząd w modelu AR.
Oto porównanie ACF(p) i PACF(p) dla miesięcznych zamknięć (czyli nie stóp zwrotu!) mWIG40:



Nie trzeba nawet sprawdzać autokorelacji stóp zwrotu, aby ocenić, że do modelowania AR(p) wystarczy p = 1.

W Gretlu wybierając KMNK, musimy wybrać zmienną zależną, którą jest tu stopa zwrotu. Regresorem będzie stopa zwrotu z opóźnieniem 1. Nie wstawiamy żadnej zmiennej jako regresor, bo poniżej jest opcja "opóźnienia", gdzie wybieramy opóźnienia dla zmiennej zależnej od 1 do 1. Dodatkowo zaznaczamy opcję "Odporne błędy standardowe". Głównym zadaniem regresji odpornej jest umożliwienie uzyskania oszacowań niepodatnych na obserwacje odstające [2]. Jak wiadomo KMNK, która wymaga normalności i homoskedastyczności składnika losowego, nie jest prawidłową metodą dla giełdowych zmian ze względu na występowanie w ich rozkładzie grubych ogonów (rzadkie zdarzenia, które są częstsze niż w rozkładzie normalnym), heteroskedastyczności i efektów ARCH. Pomiędzy tymi 3 cechami nie ma ścisłej zależności, istnieje jednak pewien związek. Modele ARCH można tak skonstruować, by pojawiły się w nich grube ogony [3], [4]. Niestety nie ma możliwości zrobienia testu stabilności przy modelowaniu GARCH, który jest oddzielną funkcją w Gretlu. Dlatego używamy regresji odpornej. Statystyka odpornościowa używa zamiast średniej np. mediany lub obciętej średniej.

Po wybraniu odpowiednich opcji otrzymałem nieistotną stałą i istotny parametr autoregresji 0,164.


1 gwiazdka oznacza istotność na poziomie 10%, dwie gwiazdki oznaczają istotność na poziomie 5%, 3 gwiazdki na poziomie 1%. W tym przypadku mamy 5% istotności, czyli pewności, że hipoteza o nieistotności parametru jest prawdziwa - jest to na tyle mało, że ją odrzucamy.
Oprócz sprawdzania istotności poszczególnych parametrów, należy też ocenić całościową istotność modelu, czyli de facto współczynnika determinacji R^2. Taki test wykonuje Statystyka F (test F-Snedecora). Widać, że wartość p dla testu F = 0,018692, zatem jest to znów poniżej 5%. Jednak to co jest uderzające, to niezwykle niska wartość R^2, która jest mniejsza od 0,03. Praktyczna zmienność stóp zwrotu w porównaniu do modelowanej stopy zwrotu staje się więc zbyt duża do praktycznego wykorzystania. Jednakże naszym celem nie jest ocena faktycznej przydatności AR, a jedynie zbadanie stabilności modelu. 
Wchodzimy Testy -> Test zmian strukturalnych Chowa. Oto wyniki:


Hipoteza zerowa to brak zmian strukturalnych, hipoteza alternatywna: zmiany są. Co się okazuje? Wartość p = 0,0026 < 0,05 , czyli odrzucamy hipotezę zerową. Faktycznie, model uległ zmianie. Na oko nie bylibyśmy w stanie tego ocenić patrząc tylko na stopy zwrotu.


Ad 2) Test CUSUM
W wielu sytuacjach niedogodnością testu Chowa jest to, że trzeba wskazywać w nim hipotetyczny punkt zwrotny. CUSUM (Cumulated Sum of Residuals) zwany czasami (nie do końca poprawnie) testem Harveya-Colliera) weryfikuje hipotezę o stabilności modelu bez potrzeby wskazywania punktu zwrotnego. Został on zaproponowany przez Browna, Durbina i Evansa w 1975 r. Wchodzimy w Testy -> test stabilności CUSUM. Dostajemy wyniki w formie tekstowej i graficznej.



Jeżeli wykres nie przekracza wartości progowych, czyli niebieskich linii to znaczy, że model można uznać za stabilny. Jak widać tak się nie dzieje, a więc model pozostaje stabilny. Dodatkowo CUSUM wskazał poziom p = 0,99 > 0,05, czyli brak istotności. Parametr ten jednak jest niezbyt przydatny, co się okaże później.


Zatem CUSUM dał przeciwstawne wyniki w stosunku do Chowa.

Trzeba jednak zaznaczyć, że CUSUM opiera się na założeniu, że zmiany strukturalne modelu mają charakter deterministyczny. Odmianą testu CUSUM jest CUSUMQ (Cumulated Sum of Squares Residuals), który może być bardziej niż CUSUM przydatny w sytuacji gdy zmiany struktury modelu mają charakter całkowicie przypadkowy (ze względu na to, że bierze kwadraty różnic to jest bardziej wrażliwy na progi istotności). Wybieramy Testy -> test stabilności CUSUMQ. Dostajemy wyniki w graficznej formie:







W tym przypadku również patrzymy czy wykres przecina niebieską linię. Trzeba przyznać, że wyniki są zupełnie nieoczekiwane. Test sugeruje, że parametry lekko zmieniły się w 2000 r. - to jest w czasach bessy oraz w 2009-2011 r., czyli w czasach niedawnej hossy. Rok 2007 i 2008 były dla testu całkowicie "normalne" pomimo gwałtownych wzrostów i spadków.

Ad 3) Test QLR
Jeśli nie wiemy kiedy nastąpi załamanie strukturalne modelu, ale spodziewamy się, że kiedyś nastąpi, wtedy możemy użyć testu QLR (Quandt Likelihood Ratio). Test ten powstał jako pierwszy ze wszystkich 3 testów, w pracy Quandta w 1958 i 1960 r. W przeciwieństwie do CUSUM zakłada on, że zmiany strukturalne mają charakter stochastyczny [5]. Wybieramy Testy -> test stabilności QLR. Otrzymujemy wyniki:

Bardzo niska wartość p wskazuje na załamanie strukturalne, które wystąpić miało w lipcu 2000 r. Wprawdzie test QLR działa prawidłowo wtedy, gdy daty punktów krytycznych są stosunkowo odległe od punktów daty początku i końca, lecz w tym przypadku rok 2000 nie jest zupełnie na początku, a ponadto data pokrywa się z tym co przedstawia test CUSUMQ.


Przeprowadzimy teraz krótką analizę AR(1) tygodniowych stóp zwrotu mWIG40 w okresie styczeń 2006:październik 2014. Tym razem Przetestujemy tylko CUSUM. Ogólna charakterystyka podana jest poniżej


Parametr autoregresji jest bardzo istotny (na poziomie 1%) i równa się 0,22 (tyle samo wynosi autokorelacja).
Pomimo że wartość p = 0,32 dla CUSUM, a więc wskazuje na nieistotność, to tak naprawdę nie decyduje o tym, czy model miał czy nie miał załamania. Należy spojrzeć na wykres czy nastąpiło przecięcie linii.






Tym razem wg CUSUM struktura modelu zmieniła się dosłownie od marca 2009 r., czyli od początku hossy. Ciekawe, że tym razem zmiana została uznana za deterministyczną, czyli tak jakby pewna zewnętrzna cecha wpłynęła na model. Dodatkowo spójrzmy na CUSUMQ


Test wskazuje na zmianę parametrów już w 2008 r. i niestabilność ta trwała aż do 2013 r. Dopiero od 2013 r. model się "ustabilizował". W przeciwieństwie do CUSUM, CUSUMQ dostarcza informacji o przypadkowej niestabilności modelu w okresach dla miejsc poza niebieską linią.

Powyższa analiza tygodniowych stóp zwrotu manifestuje "wewnętrzne sprzeczności" i trudności w ocenie ekonometrycznej. Stacjonarność danych, wysoka autokorelacja i istotność statystyczna parametrów stoją w sprzeczności do testów Chowa, QLR, CUSUM, CUSUMQ, które często odrzucają stabilność modelu.


Literatura:

[1] E. M. Syczewska, Wprowadzenie do modeli ARMA/ARIMA, W-wa 2011,
[2] D. Korniluk, Metody regresji odpornej, Prezentacja w ramach spotkan zespołu przygotowujacego sie do Econometric Game 2014, 2014
[3] B. O. Bradley and Murad S. Taqqu, Financial Risk and Heavy Tails, 2001
[4] D. Politis, A Heavy-Tailed Distribution for ARCH Residuals, 2004
[5] Discussion Of The Paper By Dr Brown, Professor Durbin And Mr Evans (jako dodatek do artykułu Browna, Durbina i Evansa "Techniques for Testing the Constancy of Regression Relationships over Time", 1975).

wtorek, 23 grudnia 2014

Ułamkowa stacjonarność? Co to za twór?

Ostatnim razem użyłem terminu "ułamkowa stacjonarność" przy testowaniu pamięci długoterminowej różnych procesów. Specjalista mógłby się przyczepić, bo w literaturze przedmiotu pojęcie to nie jest znane, wydaje się nawet dziwaczne. Stacjonarność procesu rozumie się po prostu jako stałość własności (tj. rozkładu prawdopodobieństwa) wraz z przesunięciem w czasie. Wszystko co nie jest stacjonarne, jest niestacjonarne - i ta definicja przez przeciwieństwo zdaje się zaprzeczać istnieniu czegoś takiego jak ułamkowa stacjonarność.

Tyle że praktyczna rzeczywistość sygnalizuje, że kwestia ta nie jest oczywista. Dla przykładu wygenerowałem najpierw proces ruchu Browna z wykładnikiem Hursta = 0.5, czyli zwykły ruch przypadkowy (2048 danych):



Przekształciłem go w stopy zwrotu, czyli uzyskałem Gaussowski biały szum:



Pierwszy surowy test to spojrzenie na korelogram czyli autokorelacje stóp zwrotu:




Brak istotnych dodatnich autokorelacji dla początkowych danych sugeruje już stacjonarność.

Sprawdziłem wszystkimi metodami w Gretlu stacjonarność stóp zwrotu: Zmienna->Testy pierwiastka jednostkowego. ADF wskazuje p=0, a więc na stacjonarność. KPSS daje statystykę 0,071 < 0,46 dla 5% istotności, co oznacza stacjonarność. Testy na ułamkowy rząd integracji mówią, że nie ma pamięci długoterminowej:

Estymator lokalny Whittle'a (m = 96)
  Estymowany ułamkowy rząd integracji = 0,0654853 (0,051031)
  Statystyka testu: z = 1,28324, z wartością p = 0,1994

Test GPH (Geweke'a, Portera-Hudaka) (m = 96)
  Estymowany ułamkowy rząd integracji = 0,095589 (0,0583993)
  Statystyka testu: t(94) = 1,63682, z wartością p = 0,1050

 m jest to rząd opóźnienia, estymowany rząd integracji to rząd pochodnej ułamkowej - jeżeli jest istotnie większy od 0, to występuje długoterminowa pamięć. Ponieważ p > 0,05 dla obu testów, widać, że brak tej pamięci.

Wnioski: to czego należało się spodziewać.

Następnie wygenerowałem proces Browna z wykładnikiem Hursta = 0.86 (2048 danych):



Jest to w zasadzie ten sam wykres co poprzednio, ale znacznie wygładzony, przypomina to filtrację za pomocą średniej ruchomej (w istocie jest w tym sporo prawdy). Można się więc domyślić, że będzie to miało wpływ na wygląd wykresu stóp zwrotu z tego procesu. Ich wykres jest następujący:



Dopóki wykładnik Hursta jest pomiędzy 0 a 1, stopy zwrotu tego procesu są teoretycznie stacjonarne [1]. Jednakże w rzeczywistości stopy te wykazują się silną autokorelacją i to wielu rzędów:



Powolny spadek ACF wskazuje - o dziwo - na niestacjonarność procesu. Należałoby jednak sprawdzić to rzetelniej za pomocą testów. Okazało się, że ADF odrzucił hipotezę niestacjonarności, gdyż  p = 0, ale już KPSS wskazał, że przy istotności=0.05 mamy do czynienia z procesem niestacjonarnym (Statystyka testu = 0,51 > 0,46). Testy na ułamkowy rząd integracji tym razem mówią, że pamięć długoterminowa występuje:

Estymator lokalny Whittle'a (m = 96)
  Estymowany ułamkowy rząd integracji = 0,426418 (0,051031)
  Statystyka testu: z = 8,35606, z wartością p = 0,0000

Test GPH (Geweke'a, Portera-Hudaka) (m = 96)
  Estymowany ułamkowy rząd integracji = 0,458116 (0,0581905)
  Statystyka testu: t(94) = 7,8727, z wartością p = 0,0000

[*]

Wnioski: Dla badanego procesu z H = 0.86 ADF odrzuca niestacjonarność (hipoteza zerowa to niestacjonarność), natomiast KPSS odrzuca stacjonarność (hipoteza zerowa to stacjonarność). Wyniki są sprzeczne.

Powyższa analiza unaocznia, że długozasięgowe korelacje mają pewien wpływ na ocenę stacjonarności procesu, zbliżając nas do pojęcia 'ułamkowej stacjonarności'.

Ciekawostka: Ding i Pei [2] oraz Wang [3] konstruują właśnie pojęcie ułamkowej stacjonarności, dowodząc, że:
1. każdy proces niestacjonarny można przedstawić jako sumę wielu ułamkowo stacjonarnych procesów;
2. biały szum staje się ułamkowym stacjonarnym procesem losowym po dokonaniu filtrowania fraktalnego (w Gretlu nazywa się to Filtr różnicowy ułamkowy).

Kurs papierów wartościowych jest więc sumą wielu ułamkowo stacjonarnych procesów, coś na zasadzie sumowania wielu fal o różnych skalach czasowych lub o różnych częstościach.


Literatura:

[1] BB Mandelbrot, JW Van Ness - , Fractional Brownian motions, fractional noises and applications, SIAM review, 1968
[2] J. J. Ding and S. C. Pei, Fractional Fourier Transforms and Wigner Distribution Functions for Stationary and Non-stationary Random Process, ICASSP, vol. 3, pp. 21-24, 2006
[3] Yen-Chieh Wang, Fractional Fourier Transform for Random Process Signal Analysis

[*] Dodatek:
Przypomnę, że wykładnik Hursta jest dany wzorem: H = d + 1/a. gdzie d - ułamkowy rząd integracji, a - parametr alpha rozkładu Levy'ego jako miara rozciągliwości rozkładu. Teoretycznie w tym wypadku a = 2, ponieważ wygenerowałem proces Browna, a ten zawsze jest gaussowski niezależnie od H. Ale na wszelki wypadek sprawdziłem metodą Mccullocha parametry rozkładu badanego procesu. Mcculloch pokazał, że a = 2, a więc to czego należało się spodziewać.
Jeśli chodzi o poziom H to oczywiście wynosi on 0.86. Zatem możemy obliczyć ile wynosi teoretyczny poziom d:
d = H - 1/a
d = 0.86 - 0.5 = 0.36

Porównując to z wynikami testu Whittle'a i GPH, widać, że testy zawyżyły poziom d, przy czym Whittle wskazał najbliższą prawdziwą wartość 0,426.
To zawyżenie wynika z faktu, że mamy tu do czynienia nie tylko z autokorelacją długoterminową (nieliniową), ale też autokorelacją krótkoterminową (liniową). Pewności tej nabrałem po przetestowaniu danych w programie Matrixer, który pozwala zbudować model ARFIMA. Najpierw w programie tym usunąłem wszystkie parametry do obliczenia z wyjątkiem d. Matrixer wyliczył d = 0.464, czyli bardzo zbliżoną do GPH w Gretlu. Następnie dodałem do obliczenia parametr AR(1). Matrixer wskazał, że parametr AR(1) = 0.148 oraz d = 0.366. A więc po usunięciu efektu autoregresji niemal idealnie dopasował wartość d do teoretycznej (0.36). Świadczy to oczywiście o tym, że w Gretlu metoda szukania długoterminowych zależności jest obciążona korelacjami krótkoterminowymi.

W związku z powyższym do wyliczeń d dla mWIG40 w ostatnim artykule należy podejść krytycznie. Przypominam, że uzyskana wartość d była naprawdę spora i wyniosła ponad 0.17. Teraz jednak widać, że mogła być to wartość zawyżona. I faktycznie, gdy przeprowadziłem jeszcze raz test na tych samych danych w Matrixerze, to bez AR(1) dostałem d = 0.14, a po wprowadzeniu AR(1) już d równa się ok. 0.08. Obydwa parametry (AR(1) i d) były istotne, ale wartość d jest dużo słabsza.

Teoretycznie ułamkowy ruch Browna (w sensie stóp zwrotu) jest stacjonarny dopóki -0,5 < d < 0,5 (a więc dla 0 < H < 1). Stąd badanie stacjonarności przez pryzmat ułamkowgo rzędu integracji polega po prostu na sprawdzeniu czy d < 0,5. Gretl jednak jak to wykazałem słabo sobie radzi w tej metodzie, a niestety obecnie nie ma możliwości budowania modelu ARFIMA i trzeba korzystać z innych oprogramowań.