poniedziałek, 13 lutego 2017

Optymalizacja wstępnych parametrów modelu log-periodycznego

Jeżeli mamy trudności z uzyskaniem wartości krytycznej w modelu log-periodycznym, wtedy możemy spróbować użyć pewnego przekształcenia, które pokazali Filiminov i Sornette [1]. W pierwotnym modelu występowały 3 parametry liniowe: A, B, C oraz 4 parametry nieliniowe a, w, tc, d. Można model przekształcić do takiej postaci, aby uzyskać 4 parametry liniowe i 3 nieliniowe. Nieliniowość powoduje, że obliczenia są bardzo skomplikowane i program komputerowy może nie dać rady rozwiązać problemu numerycznego. Dodatkowo nieliniowość tych parametrów wywołuje quasi-periodyczność, przez co funkcja może mieć wiele punktów optymalnych.

Przypomnę, że pierwotna wersja była następująca:

(1)

Przekształcenie modelu (1) da postać:

(2)

Wzór (2) można zapisać także jako:

(3)

Taka postać jest nawet łatwa do zapamiętania, bo jest to suma funkcji liniowej (A+Bx) oraz funkcji cosinus i sinus ((C1cos + C2sin)x).
Kod w Gretlu na nls (NMNK) dla modelu (3) - nie uwzględniając jeszcze parametrów wstępnych - byłby następujący:
 
...................................
l_Zamkniecie = A+B*(tc-time)^a+C1*(tc-time)^a*cos(w*ln(tc-time))+C2*(tc-time)^a*sin(w*ln(tc-time))

deriv A = 1
deriv B = (tc - time)^a
deriv C1 = cos(w*ln(tc - time))*(tc - time)^a
deriv C2 = sin(w*ln(tc - time))*(tc - time)^a

deriv a = B*ln(tc - time)*(tc - time)^a + C1*ln(tc - time)*cos(w*ln(tc - time))*(tc - time)^a + C2*ln(tc - time)*sin(w*ln(tc - time))*(tc - time)^a

deriv w = C2*ln(tc - time)*cos(w*ln(tc - time))*(tc - time)^a - C1*ln(tc - time)*sin(w*ln(tc - time))*(tc - time)^a

deriv tc = B*a*(tc - time)^(a - 1) + C1*a*cos(w*ln(tc - time))*(tc - time)^(a - 1) + C2*a*sin(w*ln(tc - time))*(tc - time)^(a - 1) + (C2*w*cos(w*ln(tc - time))*(tc - time)^a)/(tc - time) - (C1*w*sin(w*ln(tc - time))*(tc - time)^a)/(tc - time)
..................................


Jednak ciągle pozostaje pytanie jakie przyjąć wstępne parametry. Poprzednio zawsze przyjmowałem A = 1, B = -1, C = 0.5, a = 0.5, a także często w = 9.06, d = -1, tc zaczynałem od tc = T+1, gdzie T jest to liczba obserwacji. Trzeba znaleźć dodatkowo C1 i C2. Zgodnie ze wzorem podanym w [1] wiemy, że zachodzi:


Natomiast dotychczas stosowałem wstępne C = 0,5, a d = -1. Wobec tego wstępne C1 i C2 może wynieść odpowiednio 0,5 oraz -0.0087.

W przytoczonej publikacji [1] pokazany jest także dość skomplikowany algorytm na utworzenie wstępnych optymalnych parametrów liniowych. Należy rozwiązać następujące równanie macierzowe:


Aby nie powielać grafiki, wkleiłem bezpośredni obraz z tej publikacji. Zmienna m to nasza zmienna a, zmienna t(i) to t, indeks i to u nas też t, N to T. Wstępne liniowe parametry A, B, C1 i C2 są de facto funkcją pozostałych nieliniowych parametrów.


Literatura:
[1] Sornette, D., Filimonov, V. - A Stable and Robust Calibration Scheme of the Log-Periodic Power Law Model, 2013.

niedziela, 29 stycznia 2017

Czy model logperiodyczny może służyć do wyceny ryzykowanych akcji?

Dla wielu spółek już samo dopasowanie kursu do modelu log-periodycznego jest trudne albo niemożliwe. Ostatnio przyszedł mi do głowy taki pomysł. Jeżeli teoria Sornette'a-Johansena [1, 2, 3] jest prawdziwa, to są okresy, w których na giełdzie jest więcej spekulantów niż inwestorów (fundamentalnych). Dla jakich spółek takie okresy mogą występować najczęściej? Dość logiczne, że dla takich, które są wysoko spekulacyjne, tj. ich zyski są trudne do przewidzenia. Weźmy np. rynek gier komputerowych. Spółka latami pracuje nad jakimś projektem i w ciągu tego okresu ponosi tylko straty lub niewielkie zyski. Bardzo często przesuwana jest premiera samej gry i w zasadzie można od początku wziąć półroczną poprawkę na wstępną datę, bo ludzie (nawet specjaliści) najczęściej nie doszacowują potrzebnego czasu na zakończenie takiego dzieła. Kiedy już gra ujrzy światło dzienne, nie wiadomo czy odniesie sukces, a jeśli wiadomo, to nie jest pewne ile zarobi. Ale jeżeli odniesie sukces, to zysk może być tak duży, że podwoi czy potroi cały kapitał własny. W takiej sytuacji nie ma innej możliwości - kurs musi mocno pójść w górę. Jeżeli jednak spółka odniesie sukces, to dopiero zaczyna się kłopot. Inwestycje zwyczajnej spółki przyniosłyby nowy poziom cyklicznych przychodów, ale tutaj zysk z inwestycji jest jednorazowy. Z tego powodu wielu inwestorów unika takich spółek jak ognia: punktem zaczepienia dla nich zawsze jest kapitał własny, a w tym przypadku punkt taki nie istnieje. Jest to spowodowane wysoką niepewnością (wariancją) zysków monopolistycznych. W artykule Kiedy większa niepewność zwiększa wartość akcji? wyjaśniłem, że dla spółek o wysokim ROE w stosunku do kosztu kapitału własnego, P/BV naturalnie będzie wyższy: ta różnica zawiera zarówno średnią stopę zmian zysków monopolistycznych jak i jej wariancję, pomijając już inne momenty centralne.

Jeżeli akcje są tak wysoko spekulacyjne, to fundamentaliści ich nie tykają, czyli wycena również być może powinna opierać się na racjonalnym modelu spekulacyjnym. Zamiast klasycznie wyceniać, powinniśmy stworzyć model log-periodyczny. Natomiast dla spółek bezpiecznych, takich właśnie jak mBank, który cyklicznie płaci dywidendy, a jego zyski są w miarę do przewidzenia i beta nie jest zbyt wysoka, model ten słabo się sprawdza, co pokazałem w poprzednim artykule. Dlatego tym razem sprawdzę NMNK w Gretlu 3 spółki: CD Projekt (CDR), CI GAMES (CIG) i 11Bit (11B).


1. CDR (T = 2021). Okres 2009-2017 do 26.01.2017. Użyłem kodu:
.......................
scalar A = 1
scalar B = -1
scalar C = 0.5
scalar a = 0.5
scalar w = 15
scalar d = -1
scalar tc = 2100

l_Zamkniecie = A+B*(tc-time)^a*(1+C*cos(w*ln(tc-time)+d))
deriv A = 1
deriv B = (tc-time)^a*(1+C*cos(w*ln(tc-time)+d))
deriv C = B*(tc-time)^a*cos(w*ln(tc-time)+d)
deriv a = (tc-time)^a*(B+B*C*cos(w*ln(tc-time)+d))*ln(tc-time)
deriv w = -B*C*(tc-time)^a*sin(w*ln(tc-time)+d)*ln(tc-time)
deriv d = -B*C*(tc-time)^a*sin(w*ln(tc-time)+d)
deriv tc = -B*C*w*(tc-time)^(a-1)*sin(w*ln(tc-time)+d)+B*a*(tc-time)^(a-1)*(1+C*cos(w*ln(tc-time)+d))
.......................

gdzie l_Zamkniecie to logarytm z kursu zamknięcia, a wszystkie parametry są takie jak opisałem poprzednio.



 Dostajemy regresję w postaci niebieskiej linii:



Chociaż model nie wskazuje, aby w najbliższym czasie doszło do krachu (tc > 2700), to akcje wydają się przewartościowane, znajdując się powyżej krzywej regresji. Zauważmy, że obecnie C/WK = 8,23. Aby wycena była zgodna z modelem C/WK musi wynieść 5,56.


2a. CIG (T = 1770). Okres 2010-01-04 - 2017-01-27. To trudniejszy przypadek, bo nie mogłem znaleźć odpowiednich wstępnych parametrów. Pytanie brzmi czy moglibyśmy się dowiedzieć, w jakim oknie czasu użyć modelu. Aby na to odpowiedzieć, można użyć tc < T, a więc znaleźć osobliwość w przeszłości. To też nie jest prosta sprawa. Nie mogłem użyć parametru w w typowym zakresie 5-20. Użyłem więc kodu:
.................
scalar A = 1
scalar B = -1
scalar C = 0.5
scalar a = 0.5
scalar w = 35
scalar d = -1
scalar tc = 1780

l_Zamkniecie = A+B*(tc-time)^a*(1+C*cos(w*ln(tc-time)+d))
deriv A = 1
deriv B = (tc-time)^a*(1+C*cos(w*ln(tc-time)+d))
deriv C = B*(tc-time)^a*cos(w*ln(tc-time)+d)
deriv a = (tc-time)^a*(B+B*C*cos(w*ln(tc-time)+d))*ln(tc-time)
deriv w = -B*C*(tc-time)^a*sin(w*ln(tc-time)+d)*ln(tc-time)
deriv d = -B*C*(tc-time)^a*sin(w*ln(tc-time)+d)
deriv tc = -B*C*w*(tc-time)^(a-1)*sin(w*ln(tc-time)+d)+B*a*(tc-time)^(a-1)*(1+C*cos(w*ln(tc-time)+d))
....................

Dostajemy następujący wykres:



Wcale to nie znaczy, że wyznaczony punkt osobliwości jest wielkością absolutną. Np. gdy ustawiłem tc = 1695, otrzymałem:



Jednak tylko dla dość wąskiego zakresu mogłem dostać jakieś wartości, więc można przyjąć, że jest to obszar rozpoczęcia nowego trendu. Dlatego przyjmę nowe okno czasu:

2b. CIG (T = 624). Okres: 2014-08-01 - 2017-01-27. Mimo przyjęcia nowego zakresu, wyznaczenie wstępnych parametrów nadal sprawia trudności. W końcu przyjąłem:

scalar A = 1
scalar B = -1
scalar C = 0.5
scalar a = 0.5
scalar w = 20
scalar d = -1
scalar tc = 2500

Wyniki:




Tym razem CIG wydaje się prawidłowo wyceniony z punktu widzenia modelu spekulacyjnego. Jest tak mimo wysokiego C/WK = 6,36.


3. 11B (T = 1554). Okres 2010-10-28:2017-01-26. Jest to dobry przykład ruchu Levy'ego, dla którego wariancja jest nieskończona. Tym razem ustalenie wstępnych parametrów było bardzo łatwe, bo wystarczyło przyjąć standardowe (np. tc = T+1):

.....................
scalar A = 1
scalar B = -1
scalar C = 0.5
scalar a = 0.5
scalar w = 9.06
scalar d = -1
scalar tc = 1555
......................

Reszta bez zmian.

Wyniki:



Sytuacja jest ciekawa z dwóch powodów. Po pierwsze Spółka zgodnie z tym modelem jest niedowartościowana mimo wysokiego C/WK = 8,45. Model wskazywałby na wzrost C/WK do 12,95. Po drugie prognozowany tc = 1602,58, a więc zostało niecałe 50 dni do szacowanego załamania. Oczywiście ostatnie testy na WIG (dla początku bessy 2007) dobrze ilustrują jak istotna jest aktualizacja modelingu. Tamte testy sugerują, że trzeba myśleć o wyjściu z rynku jeśli tp = tc - Z, gdzie tp to data prognozowania, tj. tworzenia prognozy, a Z to liczba dni, która waha się między 5 a 7.
To czy hipoteza "spekulacyjnej wyceny" się sprawdzi, przekonamy się w najbliższych tygodniach.


Literatura:

[1] Sornette, D., Johansen, A. - Critical Market Crashes, 1999,
[2] Sornette, D., Johansen A., Ledoit, O. - Predicting Financial Crashes Using Discrete Scale Invariance, Feb 2008,
[3] Sornette, D., Johansen A. - Significane of log-periodic precursors to financial crashes, Feb. 1, 2008.

poniedziałek, 16 stycznia 2017

Zmagania z modelem log-periodycznym

Do teorii log-periodycznej od dawna byłem sceptycznie nastawiony, a swoją krytykę wyraziłem kiedyś w tym artykule. Z drugiej strony jestem zwolennikiem teorii chaosu (choć może nie deterministycznego), a więc także nieokresowych i nietrywialnych cykli. Ekonofizyka dostarcza także uzasadnienie powstania tzw. log-periodyczności. Nie jest wykluczone, że dla niektórych okresów lub niektórych walorów, teoria ta może się sprawdzić i tym samym wspomóc inwestora czy spekulanta w podjęciu trudnych decyzji. Co więcej, Sornette i Johansen [1] rozważają najpierw w swoim modelu "ograniczoną racjonalność" traderów (która wynika nie tylko z ułomności ludzkiego umysłu, ale także z natłoku informacji), a następnie [2, 3, 4] tworzą model w ramach racjonalnych oczekiwań, zgodnie z którym moment krachu można prognozować, ale tylko z pewnym prawdopodobieństwem. Ogólnie autorzy twierdzą, że gracze sami "się nakręcają" i kreują dodatnią warunkową oczekiwaną stopę zwrotu, która staje się tym wyższa im bliżej jest krachu. Ale staje się ona wyższa nie dlatego, że gracze "nakręcają się" w sposób nieracjonalny, ale dlatego, że gracze wymagają od rynku wyższej oczekiwanej stopy zwrotu właśnie za to, że się nakręcają. A nakręcają się w sposób racjonalny, dlatego że zdają sobie sprawę, że - w przypadku trendu rosnącego - kupując w takiej fazie, więcej ryzykują. Im bliżej krachu, tym więcej ryzykują, a więc wymagają większego zwrotu. Rozdzielając oczekiwaną stopę zwrotu na prawdopodobieństwo i wielkość ruchu, można powiedzieć, że prawdopodobieństwo dalszego (czyli warunkowego) wzrostu spada, natomiast zakres ruchu rośnie. Natomiast niewarunkowa stopa zwrotu pozostaje równa zero. Aby dobrze zrozumieć różnicę między warunkową a niewarunkową stopą zwrotu, dobrze najpierw przeczytać wpis Smarujący estymator.

Autorzy podkreślają, że nie wystarczy tylko raz zastosować model, który wyprognozuje moment krytyczny i czekać do tego okresu. Jeżeli model wskazuje moment krytyczny za rok, to w ciągu tego roku może się struktura zmienić, więc należy aktualizować model. To jest dodatkowy stochastyczny element ich teorii.

Jest wiele różnych wersji modelu log-periodycznego, ale właściwie wszędzie pojawia się wzór, który w uproszczeniu wyprowadzę idąc tokiem rozumowania Kutnera [5]. Ogólnie biorąc musimy połączyć dwie koncepcje: jedną ze statystyki - rozkład Levy'ego-Pareto i jedną z fizyki - temperatury krytycznej. Podzielę to na dwie części.

1) Rozkład Pareto-Levy'ego jest związany z grubymi ogonami, które dobrze są nam znane, więc nie będę się o nim rozpisywał. Jeśli cenę akcji P(t) da się opisać rozkładem proporcjonalnym do:

(1)

to mamy do czynienia z rozkładem Pareto dla t > 0. Zwróćmy uwagę, że z jednej strony wzór (1) opisuje prawdopodobieństwo (zmiany) ceny, ale z drugiej wartość oczekiwana będzie proporcjonalna do tego prawdopodobieństwa. Zgodnie ze wzorem (1) dla każdego okresu t mamy inne prawdopodobieństwo. Gdybyśmy przemnożyli obie strony przez cenę P(t), to dostalibyśmy wartość oczekiwaną ceny dla okresu t. Stąd w miejsce prawdopodobieństwa możemy wstawić po prostu oczekiwaną cenę.

Znaną własnością rozkładu Pareta jest równanie skalowania:

(2)

gdzie L^a to czynnik skalujący funkcję P(t).

Najogólniejszym rozwiązaniem osobliwym (tzn. rozwiązaniem, które nie posiada zapisu zgodnego z całką ogólną - gdyż mamy tu do czynienia z równaniem różniczkowym zwyczajnym) równania skalowania (2) jest

(3)
Funkcję F(u) można rozwinąć w szereg Fouriera:

(4)


2) W złożonych systemach istnieją różne fazy, w których pewne parametry fizyczne / ekonomiczne pozostają względnie stałe. W fizyce takimi parametrami są temperatura lub ciśnienie. Np. faza topnienia następuje w temperaturze poniżej 0 stopni Celsjusza, faza zamarzania - od 0 stopni wzwyż, a faza wrzenia - 100 stopni. Hołyst, Poniewierski i Ciach definiują fazę jako "makroskopowo jednorodny stan układu, odpowiadający danym parametrom termodynamicznym." [6]. Dalej: "Proces przekształcania się jednej fazy w inną nazywa się przejściem fazowym lub przemianą fazową". Różne fazy mogą ze sobą współistnieć (stąd często w przejściach między jesienią a zimą oraz zimą i wiosną utrzymuje się temperatura 0 stopni), ale w przypadku ciągłych przejść fazowych, jedna faza przechodzi bezpośrednią w inną. "Ciągłe przejścia fazowe mają na ogół charakter przemiany typu porządek-nieporządek. Przemianą tego typu jest np. przejście paramagnetyk-ferromagnetyk, występujące w materiałach magnetycznych, takich jak żelazo. Uporządkowanie dotyczy tutaj mikroskopowych momentów magnetycznych umiejscowionych w węzłach sieci krystalicznej. W fazie paramagnetycznej orientacje momentów magnetycznych są chaotyczne, więc jeśli nie ma zewnętrznego pola magnetycznego, to magnetyzacja układu znika. Jest tak wówczas, gdy temperatura układu jest wyższa od pewnej charakterystycznej temperatury Tc, zwanej temperaturą Curie. Natomiast w temperaturze niższej od Tc, mikroskopowe momenty magnetyczne porządkują się spontanicznie wzdłuż pewnego wspólnego kierunku, i magnetyzacja przyjmuje wartość różną od zera. Istotne jest to, że magnetyzacja pojawia się spontanicznie, bez udziału zewnętrznego pola magnetycznego, jako efekt oddziaływania pomiędzy mikroskopowymi momentami magnetycznymi." Tak więc temperatura Curie to temperatura krytyczna dla przemian magnetycznych, która charakteryzuje jednocześnie obszar krytyczny.

Kutner pisze, że w obszarze krytycznym większość wielkości fizycznych zmienia się w zależności od temperatury T według prawa potęgowego:

(5)

gdzie Tc to temperatura krytyczna.

Aby zastosować pojęcie temperatury krytycznej do rynku, zastąpimy temperaturę T czasem t. Jest to chyba najważniejszy moment w zrozumieniu zastosowania fizyki do giełdy: temperaturę określamy jako wielkość proporcjonalną do czasu.

Punkt 1, który jest czysto matematyczny, wydaje się zupełnie logiczny, ale zawiera założenie, że szukamy tylko rozwiązania osobliwego (a nie ogólnego). To założenie jest potrzebne, aby spełniona została interpretacja fizyczna z pktu 2, mianowicie istnienia temperatury krytycznej, a więc dla ekonofizyki - czasu krytycznego. Natomiast czas krytyczny jest potrzebny po to, aby spełniona została koncepcja "racjonalnej (lub ograniczenie racjonalnej) bańki spekulacyjnej".

Technicznie rzecz biorąc zmienna (t - tc) stanie się ujemna, bo wielkość tc z założenia ma być większa od każdego t w którym operujemy: model ma prognozować przyszłość, a nie przeszłość. Dlatego od razu użyjemy zmiennej -(t-tc) = (tc - t). Następnie, nie interesuje nas za bardzo wykres z argumentami tc-t, bo przecież chcemy widzieć normalny przebieg czasu na wykresie. Punkt krytyczny tc jest parametrem do oszacowania, a więc stałą. Dlatego stworzymy zmienną P(t) zamiast P(tc - t). Stąd za P(tc - t)  podstawimy P(t). Jeżeli uprościmy szereg Fouriera w (4) do pierwszego rzędu, dostaniemy w końcu model log-periodyczny:

(6)

gdzie A, B, C, d to pewne stałe oraz

(7)






Do oszacowania parametrów zastosuję NMNK w Gretlu, której praktyczne zastosowanie szczegółowo opisałem w poprzednim artykule
Nieliniowa metoda najmniejszych kwadratów dla modelu trendu. Algorytm Levenberga–Marquardt może w końcu przydać się do czegoś praktycznego, bo poprzednie modelowanie trendu nieliniowego za pomocą funkcji wykładniczej, jak pokazałem, nie ma większego sensu. W przypadku skomplikowanych nieliniowych trendów z cyklami i krachami sprawa wygląda zupełnie inaczej. Tak więc algorytm jest analogiczny do tamtej funkcji. Algorytm Levenberga–Marquardt został także użyty przez Sornette'a i Sammisa [7], którzy jako jedni z pierwszych używali modelu (6) - ich pierwsze próby dotyczyły przewidywania trzęsień ziemi.

Kiedy zacząłem testować (6) szybko okazało się, że wiele zależy od ustawienia początkowych parametrów: a, d , w i tc. Algorytm wprawdzie sam je szacuje, ale potrzebuje danych wejściowych. W wielu sytuacjach model nie zbiega do wybranej funkcji. Stąd warto posłużyć się wynikami uzyskanymi przez innych badaczy. Drożdż, Grummer, Ruf i Speth [8] empirycznie dowodzą, że dla większości znaczących, historycznych finansowych zdarzeń, preferowana wartość parametru L = 2. Na podstawie wzoru (7) dostaniemy więc, że preferowane w równa się ok. 9,06. Następnie, w [9] autorzy stwierdzają, że jeśli 0 < a < 1, to cena w skończonym czasie dojdzie do tc. I jeszcze na koniec znalazłem jedną z najnowszych publikacji Sornette'a i Filimonova [10], którzy po zebraniu wielu danych, skompilowali parametry do następujących ograniczeń:
0.1 ≤ a ≤ 0.9,
6 ≤ w ≤ 13,
|C| < 1
B < 0.

Jeśli chodzi o ostatnią nierówność, to jest prawdopodobne, że B > 0 dla trendu spadkowego - znalazłem w [11]. Piszę prawdopodobnie, bo autor popełnił tam błąd we wzorze na w, więc nie mam pewności czy gdzieś indziej nie ma błędu.

Jeśli chodzi o początkowe tc, to jak sądzę preferowane będzie T+e gdzie T to liczba obserwacji, a e to liczba dni, po której spodziewalibyśmy się załamania (lub ewentualnie zmiany bessy na hossę) od momentu T.

Modele log-periodyczne w trzech ostatnich wspomnianych pracach opierają się na logarytmach ceny, czyli mają postać:

(8)
Tego właśnie modelu użyję.


Przykład 1. WIG: notowania dzienne 1994-13.01.2017 (T = 5704 obserwacji). Przy preferowanych parametrach modelu nie mogłem otrzymać. Dopiero po wpisaniu następujących wstępnych parametrów:
A = 1, B = -1, C = 0.5, a = 0.5, w = 15, d = -1, tc = 5800 uzyskałem model:



 Oszacowane tc = 7405,03 (p-value 1%). Z tego punktu widzenia krach nam nie straszny: 7405 - 5704 = 1701 dni pozostało do ewentualnego krachu. Jest to oczywiście tak odległa perspektywa, że nie ma sensu brać tego na poważnie. 

Zamiast tego cofnijmy się do końca czerwca 2007. Liczba danych = 3316. Ustawione preferowane początkowe parametry wystarczyły: A = 1, B = -1, C = 0.5, a = 0.5, w = 9.06, d = -1, tc = 3317. Tym razem model robi wrażenie:



Otóż oszacowane tc wyniosło 3321,68 (p-wartość 1%). Natomiast faktyczna bessa zaczęła od 3322 dnia (7.07.2007). Trafiło w punkt. W momentu prognozy do oszacowanego dnia załamania pozostało 5-6 dni (3321,68 - 3316).
Wydawałoby się, że potrafimy przewidywać przyszłość. Nie do końca. Gdy ustawimy dane do końca maja (3296 danych), to dostaniemy tc = 3301,19, czyli o 20 dni za wcześnie, ale 5 dni od dnia, w którym dokonalibyśmy prognozy (3301-3296=5). To samo, gdy ustawimy datę końcową koniec kwietnia (3275): dostaniemy tc = 3282,11. Tym razem jednak czas do prognozowanego krachu jest nieco większy (3282-3275=7 dni). Weźmy jeszcze koniec grudnia 2006 (3192 danych). Dostałem tc = 3213,42. Czas do szacowanego krachu jest znacznie większy (3213-3192=21 dni).

Zwracam uwagę, że nie ma znaczenia jaką ustawimy datę początkową tc, zawsze wyjdzie to samo. Odpowiednia wstępna wielkość tc jest potrzebna tylko po to, by algorytm "załapał" ten parametr. W powyższych przykładach najczęściej wpisywałem tc = T+1.

Tak więc decydującym jest ciągła aktualizacja modelu (8), tak jak wcześniej podkreślali to Sornette i Johansen. A mimo to możemy uzyskiwać fałszywe sygnały, o czym przekonamy się, generując prognozę w dniu poprzedzającym początek bessy, czyli 6.07.2007: znów przesunęłaby się prognoza krachu o ok. 5-6 dni (T = 3321 oraz tc wyniósłby 3326,6, czyli 3326,6-3321=5,6). Ta sytuacja będzie codziennie się powtarzać, przez cały lipiec, a nawet w kolejnych miesiącach.

Przykład 2. Mbank: roczne 1994-2015 (T = 22 obserwacje). Wróćmy jeszcze raz do mbanku, który ostatnio modelowałem regresją nieliniową, tyle że zwykłą funkcją wykładniczą. Tym razem możemy znaleźć o wiele dokładniejsze przewidywania przy tych samych danych. Aby dostać model użyłem: A = 1, B = -1, C = 0.5, a = 0.5, w = 16, d = -1, tc =23. NMNK oszacował tc = 24,87, a wykres jest taki:



Jednakże obecnie możemy już dodać rok 2016. Pytanie czy wielkość tc zmieni się? Użyłem tych samych wstępnych parametrów. NMNK oszacował ponownie tc = 24,87, a wykres jest taki:


 Czyli wg modelu mamy spodziewać się załamania za 1-2 lata.

Aby jednak nie napawać się zbytnio optymizmem, pokażę teraz, że model (8) kiepsko radził sobie z dziennymi wahaniami.


Przykład 3. Mbank: dzienne 1994-koniec 2013 (T = 4941 obserwacji): Wpisałem następujące wstępne parametry:
A = 1, B = -1, C = 0.5, a = 0.5, w = 9.06, d = -1, tc = 5000. Uzyskałem model o wykresie:



Ogólnie model na rysunku wydaje się całkiem dobry. Parametr krytyczny tc został oszacowany na 5507,77 i jest istotny statystycznie na poziomie 1% (tak samo jak reszta parametrów). czyli krachu powinniśmy oczekiwać po 5507,77-4941 = 566,77 dniach. Dodatkowo: Błąd standardowy reszt  = 0,252757, Skorygowany R-kwadrat   0,890034. Co więc się stało po tym czasie? Poniżej zamieściłem ten sam wykres powiększony o 567 dni.



A do dziś mamy:



Już dawno po "krachu". Model powinien przewidywać go w momencie gdy tworzyliśmy prognozę, bo od tego momentu trend się zmienił. Ale można byłoby starać się obronić model, bo kurs od momentu wykonania prognozy, czyli końca 2013, przez kolejne pół roku stał w miejscu, więc należałoby aktualizować model i być może wtedy uzyskalibyśmy moment wyjścia? Nic z tego. Przesuwając dane do końca lipca 2014 (5087) dostaniemy tc = 5588,87 i wykres:


Zatem czekalibyśmy, bo przecież różnica 5589-5087 = 502, czyli jeszcze ok. 2 lata dni roboczych.

Powyższe bardzo skrótowe testy pokazały, że model log-periodyczny może rzeczywiście służyć jako wskaźnik "temperatury rynku" czy zwiększonej niepewności dla racjonalnego inwestora, tak że spekulacja zaczyna przeważać nad wyceną fundamentalną. Ale kontrowersyjne jest twierdzenie, że pozwala przewidywać zmianę trendu. Możliwe, że rynek zachowuje się tak jak to opisuje teoria Sornette-Johansena [2, 3, 4]. Fakt pozostaje faktem, że analiza WIG wykazała, że im bliżej był dzień zmiany trendu w 2007, tym model (8) wskazywał na coraz bliższą datę od momentu prognozy, przy czym krańcowym czasem pozostałym do załamania był mniej więcej 1 tydzień.

Na sam koniec stworzę prognozę dla indeksu USA dla częstości rocznej i miesięcznej.

Przykład 4. S&P500 roczne: 1933-2016 (T=84 dane). Wstawiłem standardowe parametry, w tym tc = 85. W odpowiedzi otrzymałem tc = 94,355. Czyli do krachu stulecia jeszcze 10 lat.


Najgorsze rezultaty model osiągnął w drugiej połowie lat 90 XX w. Można to nawet potraktować jako argument za przewartościowaniem akcji w tym okresie.

Przykład 5. S&P500: miesięczne 1933-koniec 2016 (T=1008). Ustawienia: identyczne co poprzednio, ale nie chciało wejść z tc = 1009 i dopiero dla tc = 1020 "załapało". Prognoza tc = 1126,55, a wykres przedstawia się:


W powiększeniu ostatnich okresów widać, że model wchodzi w fazę wzrostu:


Jeżeli prognoza się sprawdzi, to czeka nas krach w USA za 119 miesięcy, czyli za 119/12 = 9,9 lat. Prognoza miesięczna idealnie więc pokrywa się z prognozą roczną.


Literatura:
[1] Sornette, D., Johansen, A. - Modeling the stock market prior to large crashes, 1999,
[2] Sornette, D., Johansen, A. - Critical Market Crashes, 1999,
[3] Sornette, D., Johansen A., Ledoit, O. - Predicting Financial Crashes Using Discrete Scale Invariance, Feb 2008,
[4] Sornette, D., Johansen A. - Significane of log-periodic precursors to financial crashes, Feb. 1, 2008,
[5] Kutner, R., Wprowadzenie Do Ekonofizyki: Niegaussowskie Procesy Stochastyczne Oraz Niedebye'Owska Relaksacja W Realu. Elementy teorii ryzyka rynkowego wraz z elementami teorii zdarzeń ekstremalnych, W-wa czerwiec 2015,
[6] Hołyst, R., Poniewierski A., Ciach A., Termodynamika Dla Chemików, Fizyków I Inżynierów, Październik 2003,
[7] Sornette, D., Sammis, C., - Complex critical exponents from renormalization group theory of earthquakes Implications for earthquake predictions, 1995,
[8] Drożdż S., Grummer, F., Ruf, F., Speth, F. - Log-periodic self-similarity an emerging financial law, 2002,
[9] Fantazzini, D. ,Geraskin, P., - Everything You Always Wanted to Know about Log Periodic Power Laws for Bubble Modelling but Were Afraid to Ask, Feb 2011,
[10] Sornette, D., Filimonov, V. - A Stable and Robust Calibration Scheme of the Log-Periodic Power Law Model, 2013,
[11] Górski, M. - Zastosowanie Teorii Log-Periodyczności W Prognozowaniu Krachów Giełdowych, 1/2014.

piątek, 23 grudnia 2016

Nieliniowa metoda najmniejszych kwadratów dla modelu trendu

W wielu pakietach ekonometrycznych można użyć Nieliniowej Metody Najmniejszych Kwadratów (NMNK), aby uzyskać bezpośrednio poszukiwany parametr. NMNK ma podobnie jak MNK podstawy teoretyczne - o ile słuszność MNK wynika z twierdzenia Gaussa-Markowa, o tyle NMNK jest słuszna z powodu uogólnionych twierdzeń Gaussa-Markowa [1, 2, 3], a także dowiedzionej asymptotycznej zgodności jej estymatora [4].

NMNK nie jest jednolitą metodą. W Gretlu stosowany jest algorytm Levenberga–Marquardt. Algorytm ten jest kombinacją dwóch znanych metod numerycznych: metody gradientu prostego oraz algorytmu Gaussa-Newtona. W metodzie gradientu prostego suma kwadratów błędów jest korygowana w kierunku spadku nachylenia funkcji (gradientu), bo tam znajduje się minimum. W metodzie Gaussa-Newtona wykorzystuje się z kolei przybliżenie funkcji za pomocą wzoru Taylora; kwadraty różnic pomiędzy prawdziwą funkcją a funkcją Taylora powinny być jak najmniejsze [5].

Na chwilę zostawmy to. Wcześniej (tu oraz tu) pokazałem, jak poprawnie retransformować model trendu, aby uzyskać warunkową oczekiwaną stopę zwrotu. Mamy następujący model trendu:

(1)


Powinniśmy zdefiniować składnik losowy. Powiedzmy, że jest on stacjonarny, posiada wartość oczekiwaną = 0 oraz stałą wariancję. Jeżeli zmienna losowa P(T) ma rozkład log-normalny, to jej wartość oczekiwana wynosi:

(2)

Dzięki NMNK możemy oszacować parametry bezpośrednio z modelu (1), a więc bez potrzeby transformacji logarytmicznej. Przykładowo, biorąc te same dane dla mbanku co w wymienionych na początku artykułach (rocznie 1994-2015; 22 obserwacje), możemy stworzyć w Gretlu model (1). Po zaimportowaniu danych, trzeba kliknąć: Model -> Nieliniowa Metoda Najmniejszych Kwadratów. Pojawi się okno i tam wpisujemy:
.................
scalar a = 1
scalar b = 1
Zamkniecie = exp(a+b*time)
deriv a = exp(a+b*time)
deriv b = exp(a+b*time)*time
.................

gdzie:
Zamkniecie - to szereg czasowy, w naszym przypadku kurs mbanku (zmienna P(T));
time - zmienna czasowa (wygenerowana w opcji Dodawanie zmiennych -> time - zmienna czasowa.
a, b - współczynniki regresji, które trzeba od początku zadeklarować. Są to skalary w naszej funkcji, które przyjmą jakąś początkową wartość, np. 1. Standardowo deklarujemy więc scalar a = 1, scalar b = 1. Mogą być one dowolne, ale może się zdarzyć, że nie wystarczy iteracji, aby zbiegły do prawdziwej wartości wg algorytmu. Jeśli tak się stanie, zmieniamy metodą prób błędów.

deriv - pochodna funkcji po danym współczynniku. Np. pochodna d(Zamkniecie)/d(a) = exp(a+b*time), dlatego jest zapis:
deriv a = exp(a+b*time).
I tak samo dla b: pochodna d(Zamkniecie)/d(b) = exp(a+b*time)*time. Dlatego jest zapis:
deriv b = exp(a+b*time)*time.

Mimo że a i b są to skalary modelu, to są one w funkcji szacowane, czyli pierwotnie są zmiennymi.

Niezbyt przyjemne. Można się obyć bez pochodnej, ale w helpie Gretl rekomenduje jej obliczanie.

Przed kliknięciem OK, możemy zaznaczyć jeszcze opcję "Odporne błędy standardowe". Wtedy OK.



Jak widać na tablicy zbieżność nastąpiła po 27 iteracjach. Gdybyśmy wpisali scalar a = 10, scalar b = 10, zbieżność nastąpiłaby dopiero po 240 iteracjach.




Uzyskany parametr b = 0,0996, natomiast wolny parametr a możemy zignorować, bo interesują nas stopy zmian. W porównaniu do logarytmicznego modelu, gdzie uzyskaliśmy logarytmiczną stopę ponad 11%, może się wydać dziwne, że teraz jest niższa wartość. Przecież otrzymaliśmy od razu stopę, która powinna być właśnie wyższa od tej z modelu logarytmicznego MNK. Co jest tego przyczyną? Stało się tak, bo nasz model jest niepełny. Założyliśmy błędnie, że wygląda jak na poniższym rysunku:


W tym modelu składnik losowy jest stacjonarny, może więc to być zmienna losowa o rozkładzie normalnym. Ta zmienna po prostu jest argumentem funkcji wykładniczej. Ale w rzeczywistości powinna być funkcją czasu. Gdy za tę zmienną losową wstawimy proces Browna, czyli skumulowaną zmienną losową o rozkładzie normalnym, to otrzymamy geometryczny proces ruchu Browna. Wtedy nasz model wygląda mniej więcej tak:



Ogólnie geometryczny proces ruchu Browna ma postać (zob. np. [6]):

(3)

gdzie B(T) jest zwykłym procesem ruchu Browna. B(0) = 0.

Proces (3) w danym punkcie czasu może przyjąć różne wartości zgodnie z pewnym rozkładem prawdopodobieństwa. Przestaje być więc rozumiany jako proces, a staje się "ruchem" - zmienną stacjonarną. Ten ruch posiada zawsze rozkład log-normalny. Oznacza to, że możemy wyprowadzić analogicznie do (1) wartość oczekiwaną tej zmiennej.


Zauważamy, że składnik losowy w modelu (1) to jest to samo co (odchylenie standardowe razy proces ruchu Browna) w modelu (3), co oznacza, że musi się tak samo przekształcać jak w modelu (2), gdzie występuje wartość oczekiwana procesu. I stąd podstawiając do (3) powyższe wyprowadzenie wariancji, dostaniemy wartość oczekiwaną geometrycznego (procesu) ruchu Browna:

(4)

W tym miejscu rodzi się pytanie czy rzeczywiście znamy wartość początkową P(0). To zależy od tego czy zaczynamy od T = 1 czy T = 0. Musimy zdecydować czy pierwsza wartość już decyduje o funkcji zależnej od czasu czy dopiero kolejna. Jeśli założymy, że P(0) jest znane, to podstawimy pierwszą cenę zamknięcia mbanku z 1994 r., czyli 30, 68. W przeciwnym wypadku P(0) musi zostać oszacowane. Wybór sposobu ma istotne znaczenie dla uzyskanych wyników: dodanie wolnego wyrazu zmienia postać modelu. Ale jeżeli wybierzemy model z wyrazem wolnym, to estymacja NMNK przyniesie ten sam rezultat co dla modelu (2) - czyli powracamy do modelu (1), który przecież uznaliśmy za niepełny. Wstawienie P(0) nie rozwiąże tego problemu, ale możemy sprawdzić do czego doprowadzi.

Wstawiamy pierwszą cenę zamknięcia P(0) = 30,68. Możemy wtedy powtórzyć całą procedurę w Gretlu, używając modelu (4):

.................
scalar b = 1
Zamkniecie = 30.68*exp(b*time)
deriv b = 30.68*exp(b*time)*time
.................

Uzyskane parametry to:


Dostaniemy wykres trendu:



Zatem u = 12,5, czyli efektywna oczekiwana stopa exp(12,5) - 1 = 13,3%.

Mimo iż obliczyliśmy parametr u, to taki estymator jest niezgodny w tym sensie, że nie będzie dążył do prawdziwej wartości parametru u. Żeby była zgodność, składnik losowy musi być stacjonarny (zob. twierdzenia w [4]). Najprawdopodobniej potrzeba użyć uogólnionej nieliniowej metody najmniejszych kwadratów, która pozwoliłaby uwzględnić rosnącą wariancję i autokorelację składnika losowego.

Ostatecznie można ten problem rozwiązać. Łatwo zauważyć, że tylko w punkcie dla T = 1 geometryczny proces ruchu Browna będzie tożsamy z modelem (1), dla którego składnik losowy jest stacjonarny. W tym punkcie obydwie krzywe się przetną. Czyli w tym punkcie model (2) równa się modelowi (4):


Wyrazy wolne możemy do siebie dostosować, aby zachować identyczną interpretację punktu startu. Dlatego przyjmiemy lnP(0) = a. Stąd:

(5)

Ze względu na stacjonarność składnika losowego dla T = 1 estymator u równy b plus połowa wariancji może stać się w końcu zgodny. Wcześniej oszacowany współczynnik b był estymatorem zgodnym (zakładaliśmy początkowo stacjonarność składnika losowego), chociaż występował w modelu niepełnym. Połowa wariancji też jest estymatorem zgodnym. Suma estymatorów zgodnych daje estymator zgodny (wynika to z własności granicy ciągów). Stąd parametr u równy b plus połowa wariancji stanie się estymatorem zgodnym. Parametr b został wcześniej oszacowany, ale wariancję również znamy, bo jest to przecież ten sam składnik, który był obliczany dla MNK. Estymowaliśmy go wcześniej za pomocą błędu standardowego reszt (tu) i wyniósł 0,346. Możemy podstawić:

u = b + 0,5var
u = 0,0996 + 0,5*0,346^2 = 0,1596

Aby uzyskać oczekiwaną (efektywną) stopę zwrotu wstawimy exp(0,1596) - 1 = 17,3%. Dla porównania w Transformacja lognormalnego modelu z nieznanym parametrem wyszła wartość zbliżona, bo 18,7%.


Literatura:
[1] Louton, T., The Gauss-Markov Theorem for Nonlinear Models, Dec. 1982,
[2] Kariya, T., A Nonlinear Version of the Gauss-Markov Theorem, Jun. 1985,
[3] Kariya, T., A Maximal Extension of the Gauss–Markov Theorem and Its Nonlinear Version, 2002,
[4] Wu, C-F, Asymptotic Theory of Nonlinear Least Squares Estimation, May 1981,
[5] Marquardt, D. W., An Algorithm for Least-Squares Estimation of Nonlinear Parameters, Jun. 1963.
[6] https://en.wikipedia.org/wiki/Geometric_Brownian_motion