Liniowy model prawdopodobieństwa vs. Regresja logistyczna

Powiedzmy, że zastanawiamy się czy zainwestować w jakieś akcje. Oczywiście standardowym narzędziem do tego jest wycena zdyskontowanych przepływów pieniężnych lub dywidend. Jeśli model wycenił spółkę na poziomie niższym niż bieżąca kapitalizacja, to dostajemy sygnał kupna. W przeciwnym wypadku dostajemy sygnał sprzedaży. Niestety wycena to problem skomplikowany i oparty na wielu założeniach. Nawet wspomagana przez analizę fundamentalną, w której bierzemy pod uwagę nie tylko czynniki wewnętrzne, ale i otoczenie firmy, może zawieść. Łatwiej zastosować samą AF opartą na wskaźnikach lub analizie dyskryminacyjnej. Zagadnienie tej ostatniej już poruszyłem (tutaj , tutaj , dla banków - tutaj) i należy stwierdzić, że jest wyjątkowo użyteczna, będąc jednocześnie techniką prostą i precyzyjną, bo mówi konkretnie czy firma jest w dobrej czy złej kondycji (w przeciwieństwie do analizy wskaźnikowej, która jest niejednoznaczna). Jej wadą jest jednak na stałe ustalenie wielkości parametrów (na podstawie testów na innych próbach z przeszłości). Podobną rolę do analizy dyskryminacyjnej odgrywają liniowy model prawdopodobieństwa oraz regresja logistyczna (model logitowy; logit), ale dostarczają one jeszcze dokładniejszej informacji: jakie jest prawdopodobieństwo, że na przykład firma jest w dobrej albo złej kondycji; jakie jest prawdopodobieństwo bankructwa. A zatem dostajemy nie tyle sygnał kupna lub sprzedaży, co prawdopodobieństwo sygnału (jakkolwiek to dziwnie brzmi). Należy dodać, że logit jest to metoda dużo bardziej skomplikowana i trudniejsza do zrozumienia i żeby ją zastosować konieczne jest użycie programu komputerowego, najlepiej do ekonometrii. Ja używam gretla.

Skoro chcemy uzyskać szansę tylko na to, że będzie dobrze lub źle, to znaczy, że zmienna Y ma tylko dwie wartości: 0 (będzie źle) lub 1 (będzie dobrze). Zacznijmy więc od najprostszego przykładu rzutu monetą, która nie musi być sprawiedliwa. Powiedzmy orzeł - 1, reszka - 0. Chcemy orła. Rzucamy raz monetą - udało się. Rzucamy drugi raz - znów orzeł. Za trzecim razem jest reszka. Rzucamy tak wielokrotnie i zapisujemy wyniki. Teraz chcemy sprawdzić prawdopodobieństwo wyrzucenia orła. Oczywiście nie potrzebujemy do tego żadnych cudów, logitów. Zwyczajnie obliczamy częstość wyrzuconych orłów i dzielimy przez liczbę wszystkich rzutów. Mamy więc to czego szukaliśmy. Powiedzmy, że zastępujemy monetę samochodem jadącym po trudnej drodze - 1 to brak wypadku, 0 to wypadek. Chcielibyśmy sprawdzić szansę, że wypadku nie będzie. Porównując auto z monetą napotykamy problem - nie możemy go testować wielokrotnie, ale tylko raz. Co zrobić? Moglibyśmy wprawdzie sprawdzić statystyki wypadków w kraju w ciągu roku i podzielić na liczbę jadących samochodów w tym samym czasie. Pomijając kwestię aktywnych samochodów, to taka metoda obliczania prawdopodobieństwa przypomina przysłowie o średniej liczbie nóg człowieka i psa. Przecież podejrzewamy, że aby doszło do wypadku, muszą zajść jakieś okoliczności, chcemy ustalić jakieś zmienne, które do tego prowadzą. Bez dokładniejszych informacji, obliczanie średniej częstości nie ma dla nas sensu. Chcemy większej precyzji, stąd pierwszy pomysł odpada. Musimy odkryć związek między jakimiś zmiennymi a zdarzeniem się wypadku, a następnie znaleźć jakiś sposób na obliczenie prawdopodobieństwa wypadku. W przypadku logitu będzie w zasadzie na odwrót: najpierw znajdziemy prawdopodobieństwo, a ono doprowadzi nas do związku między zmiennymi.

Liniowy model prawdopodobieństwa 

Skoro szukamy zmiennej zależnej i niezależnej, to powróćmy na chwilę do monety. Na wynik monety wpływała tylko jej budowa, a nie to, jak nią rzucamy. Budowa mogła powodować większe szanse wyrzucenia orła od reszki lub na odwrót. Co więcej, gdyby moneta była idealnie sprawiedliwa, to nawet budowa nie wpływałaby na zdarzenie oprócz tego, że pozwala jedynie na 0 lub 1. Innymi słowy w tym przypadku Y = X, gdzie Y to wyrzucenie 0 lub 1, a X to zmienna wpływająca na Y i przyjmuje to samo co Y. Jeżeli podejrzewamy, że moneta jest niesprawiedliwa, to to równanie musi się zmienić, np. na Y = a + b*X. Teraz moneta została jakby rozbita na dwie części: same zdarzenia Y oraz zdarzenia, które X wpływają na Y, ale nie są już tożsame z Y.
 
Wróćmy do przykładu z samochodem. Podobnie jak z monetą, musimy spytać co będzie miało wpływ na wypadek. Gdyby Y = X, to znaczyłoby, że na wypadek wpływa jedynie sama budowa samochodu i jego "wewnętrzne zagrożenia" albo sam kierowca. Głównym czynnikiem tworzącym zagrożenie są umiejętności kierowcy i jego stan psychiczny. Dla uproszczenia połączymy obie te cechy w jedną zmienną X: kierowca o wysokich umiejętnościach i w dobrym stanie psychicznym (1) lub nie spełniający tych wymagań (0). Gdyby nie było innych zakłócających czynników, to równanie Y = X mogłoby zostać spełnione: dobry kierowca (X = 1) nigdy nie miałby wypadku (Y = 1), natomiast zły kierowca (X = 0) zawsze miałby wypadek (Y = 0). Ta idealna korelacja rzecz jasna nie jest prawdziwa, ale pozwala dobrze zobrazować początkową sytuację. Znalezienie prawdopodobieństwa sprowadzałoby się do określenia z jakim kierowcą mamy do czynienia. I nagle spostrzegamy, że szansa musiałaby wynosić dokładnie 1 albo 0, bo utożsamiliśmy idealnie zdolności kierowcy z zagrożeniem wypadku. Skąd ta różnica w porównaniu z monetą? Rzucanie sprawiedliwej monety (bo o takiej teraz mówimy) to czysty los, natomiast zdolności i stan kierowcy losowe nie są. Można te zdolności oczywiście traktować jak zmienną losową, jednakże to właśnie fakt, że nie są w pełni losowe, pozwala na lepsze przewidywanie zdarzeń.

Oczywiście taka idealna zależność nie występuje, dlatego równanie Y = X należy także tutaj poprawić. Pytanie jak? Zacznijmy od liniowej zależności Y = a + b*X. Jeżeli X = 0, to Y = a, zaś gdy X = 1, to Y = a + b. To jednak ciągle oznacza doskonałą korelację równą 1: X to zawsze przekształcony Y. Załóżmy, że a = 0. Wtedy zawsze Y = 0 gdy X = 0. Dla b = 1, gdy X = 1, Y = 1. Zmienne zachowują się identycznie. W rzeczywistości relacja X-Y jest statystyczna, więc powinna być modelem regresji. To znaczy, że powinna zawierać dodatkowo składnik losowy. Czyli moglibyśmy zapisać

gdzie ε to składnik losowy. Wtedy dostajemy linię regresji jako warunkową wartość oczekiwaną Y względem X:

Zauważmy, że wartość oczekiwana Y równa się prawdopodobieństwu zajścia zdarzenia Y:

(1)

gdzie p to prawdopodobieństwo, że Y przyjmie wartość 1.

W tym miejscu ważne spostrzeżenie: skoro p to prawdopodobieństwo, że Y = 1, to znaczy, że:

(2)

 
Analogicznie warunkowa wartość oczekiwana Y będzie się równać prawdopodobieństwu warunkowemu zajścia zdarzenia Y pod warunkiem X:

(3)

Znów ważne jest to, że mamy tu do czynienia z prawdopodobieństwem warunkowym p, że Y wyniesie 1:

(4)

To kluczowy moment do zrozumienia dalszego wywodu. Warunkowa wartość oczekiwana Y to pewna średnia pod warunkiem, że X wyniesie (0 lub 1). A skoro średnia, to jest to tylko coś między 0 a 1. Ta średnia będzie jednak równa prawdopodobieństwu, że Y właśnie wyniesie 1 (pod warunkiem, że X równa się coś). To trudny, ale ważny element teorii. 

Jeżeli regresja Y byłaby liniowa, to dostajemy:

Taki model nazywa się liniowym modelem prawdopodobieństwa. Parametry a i b znajdujemy normalnie metodą MNK lub inną.

Przykład 1.
Stworzyłem przykładowy zestaw 50 danych: X - Czy kierowca ma zdolności? 1 - Tak, 0 - Nie. Y - Czy brak wypadku? 1 - Tak, 0 - Nie.

Wchodzimy do gretla. Model ->Klasyczna metoda najmniejszych kwadratów. Zaznaczamy jako zmienną niezależną Czy_kierowca_ma_zdolnosci, a jako zależną - Czy_bez_wypadku. Wynik:

 Z uzyskanego okna wyników wchodzimy w Prognozę. Otrzymujemy:

Model wskazuje teoretyczne prawdopodobieństwa (prognozy), jego błędy ex-ante i przedział ufności. Weźmy pierwszą obserwację: był wypadek (0), a prognoza, że nie będzie wypadku to 30%. Skąd wiem, że o tym mówi prognoza? Właśnie stąd, co wcześniej napisałem, że prognoza to prawdopodobieństwo, że Y = 1 pod warunkiem X. To pod warunkiem X powoduje, że nie powiedziałem wszystkiego. Prognoza, że nie będzie wypadku wynosi 30%, ale tylko dla tej jednej obserwacji. A ta jedna obserwacja spełnia pewne warunki. Tym warunkiem jest zmienna X. To znaczy, że dostajemy szansę 30%, że nie będzie wypadku, ale tylko wtedy, gdy X równa się... i tu musimy wrócić do pierwszego rysunku, gdzie jest lista wartości zmiennej Czy_kierowca_ma_zdolnosci. Pierwsza wartość wynosi zero. Czyli wiemy już, że gdy X = 0 (gdy kierowca nie ma zdolności) to prawdopodobieństwo, że Y = 1 (że nie będzie miał wypadku), wynosi tylko 0,3, co jest równoznaczne z tym, że będzie miał wypadek na 70% (1 - 0,3 = 0,7).

 Dla drugiej obserwacji Y = 1, a p = 0,6. Sprawdzamy pierwszy rysunek, wartość X(2) = 1. Czyli wiemy już, że gdy X = 1 (gdy kierowca ma zdolności) to prawdopodobieństwo, że Y = 1 (że nie będzie miał wypadku), wynosi 0,6.

Żeby pokazać zmienne niezależne i prognozę, zapisuję prognozy (+ w oknie prognozy) jako zmienną Pr_Czy_bez_wypadku_lin1. W oknie głównym zaznaczam tę nową zmienną razem z Czy_kierowca_ma_zdolnosci, i wyświetlam tabelę:

Faktycznie, dostaję to czego oczekiwaliśmy.

Podsumowując ten przykład, jeżeli mamy dobrego kierowcę, to spokojnie przejedziemy z prawdopodobieństwem 0,6, a jeżeli mamy złego kierowcę, to będziemy mieć wypadek z prawdopodobieństwem 0.7. Widzimy więc, że nie ma symetrii między odwrotnymi przypadkami.

Przykład 2.
Spróbujmy model poprawić, dodając drugi czynnik - otoczenie. Nowa zmienna niezależna X2 będzie pytać czy inni kierowcy mają zdolności? Tak - 1, Nie - 0. Dane tak samo własnoręcznie stworzyłem. Tabela wszystkich zmiennych z wartościami jest teraz taka:

Wszystko robimy tak samo jak poprzednio, tylko dodajemy nową zmienną niezależną do modelu.

Model się znacznie poprawił (AIC, BIC, HQC). Spójrzmy na prognozy:

Jak widać prognozy wypadków już dużo bardziej się zmieniają niż w pierwszym modelu. Mimo wszystko możemy bez nowych obserwacji sami wychwycić p. Zapiszmy prognozę w bazie danych (znak + w oknie prognozy), np. o nazwie Pr_Czy_bez_wypadku_lin2. W oknie głównym zaznaczmy zmienne niezależne i nowo dodaną zmienną. Wyświetlamy:

Gdy X1 = 0 i X2 = 0, to p = 0,13. Gdy X1 = 1 i X2 = 1, to p = 0,826. Gdy X1 = 1 i X2 = 0, p = 0,4. Gdy X1 = 0 i X2 = 1, p = 0,554. Te wielkości potem się powtarzają, stąd łatwo zgadniemy jakie będzie p dla nowych obserwacji (nowych kierowców). Jednak przy 3 zmiennych będzie to już na tyle skomplikowane, że nie opłaca się samemu tego robić i lepiej skorzystać z automatycznej prognozy.

Powiedzmy, że dostajemy nową obserwację X, tzn. nowego kierowcę, o którym wiemy, że bywa pijany, czyli X1 = 0, a do tego kierowcy, którzy jadą obok niego łamią prawo, bo stracili prawo jazdy, tj. X2 = 0. Nie wiemy, czy nastąpi wypadek, czyli nie posiadamy wartości Y. Dostajemy realną sytuację, w której przewidujemy przyszłość, tzn. chcemy się dowiedzieć, jaka jest szansa, że nie będzie wypadku. W Dane -> Dodaj obserwacje dodajemy 1 obserwację, zaznaczamy zmienną X1, X2 i Dane - > Edycja wartości. Możemy wpisać 51 obserwację, tj. X1 = 0, X2 = 0. Co ważne, nie dodajemy Y. Teraz ponownie wchodzimy w KMNK, zaznaczamy wszystko tak jak poprzednio i wchodzimy w prognozę. Dostałem:

Otrzymaliśmy to, czego oczekiwaliśmy, tj. p = 0,13. Jak wyżej napisałem, można samemu do tego dojść przy małej liczbie zmiennych, ale przy większej już taka procedura byłaby niezbędna. 
 
W ten sposób uzyskujemy prosty model do szacowania poziomu zagrożenia.

Liniowy model prawdopodobieństwa jest prosty i dość logiczny, ale ma dwie wady. Po pierwsze, w tablicy prognoz, przedział ufności zawiera liczby ujemne, a także powyżej 1. A przecież prawdopodobieństwo nie może być ujemne ani wyższe od 1. W omawianym przypadku problem ten ma charakter czysto teoretyczny, bo nie dostaniemy prognoz poniżej 0 ani powyżej 1. Niemniej mogą się pojawić przykłady, że gdy takie wartości się pojawią.  

Druga wada jest z punktu widzenia analizowanego przykładu poważniejsza. Bo na jakiej właściwie podstawie założyliśmy, że zależność Y-X jest liniowa? I stąd, dlaczego prawdopodobieństwo zdarzenia Y ma być funkcją liniową X? Gdy stosujemy klasyczną regresję liniową, to zazwyczaj patrzymy graficznie czy obserwacje X układają się liniowo względem Y. A jak się dzieje w tym przypadku? Wykres punktowy wygląda następująco:

Oczywiście wiele z 50 danych nakłada się na siebie, bo wszystkie mają tylko 0 lub 1.

Regresja logistyczna 

Czy można rozwiązać nasz problem bez założenia liniowej zależności X-Y? Powiedzmy, że mamy jakąś ogólną statystyczną zależność X od Y:

(5)

 gdzie ε to składnik losowy.

Rozwiązanie tego zadania wiąże się z wykorzystaniem zasady maksymalnej entropii. To ona doprowadza do wzoru na regresję logistyczną, który wyprowadziłem poprzednim razem:

(6)

gdzie:
p(i) - prawdopodobieństwo i-tego zdarzenia dla zmiennej X, która przyjmuje wartości x(i, 1), x(i, 2), ..., x(i, K).
Reszta to stałe.
Ponieważ są dwa zdarzenia Y, to i = 1, 2. Natomiast dla X jest w sumie K*i + 1 zdarzeń.

Mamy tu czynienia z K-wymiarową zmienną X. W naszym przypadku K = 1:

(7)

 W uproszczeniu:

(8)

Model (7) można też przekształcić do postaci:

(9)

Ponieważ szukamy p(Y = 1 | X), to (7), dla którego i = 1 (bo i = 1, 2), zapiszemy:

(10)

Pamiętajmy przy tym, że ciągle jest spełniona tożsamość (4), stąd:

(11)

Sumarycznie moglibyśmy zapisać:

Zmienną z formalnie należy traktować jako stałą wartość. Dla wszystkich zdarzeń, tj. X, p stanie się funkcją prawdopodobieństwa - dystrybuantą rozkładu logistycznego:

(12)

Powstaje nieco zagmatwana sprawa. Zgodnie z (12) zmienna Z ma rozkład logistyczny. Ale jedyną zmienną losową tworzącą Z jest X proporcjonalny do Z, zatem X też musi mieć rozkład logistyczny (tzn. rozkład prawdopodobieństwa zmiennej Z jest rozkładem logistycznym). Co się dzieje jednak z Y? Zmienna ta nie może być tożsama ze zmienną Z, bo by idealnie korelowała z X bez błędów losowych. To znaczy, że brakuje jej składnika losowego w Z.  Pytanie czy ten składnik można po prostu dodać do Z? Okazuje się, że dla jednej zmiennej X można tak zrobić, ale nie gdy występuje więcej zmiennych. Np. jeśli X = 1, a beta(0) = 0, beta(1) = 0,5 to Y = 0,5 + składnik losowy. Pamiętajmy, że Y to zmienna zero-jedynkowa, więc składnik losowy przyjmie wartość 0,5 albo -0,5. Wartość oczekiwana składnika losowego równa się zero i powinien być niezależny od X, co powoduje, że E(Y | X) = 0,5. Jeżeli jednak występują dwie zmienne niezależne, X1 i X2, wtedy jeśli beta(0) = 0, beta(1) = 0,5 oraz beta(2)  = 0,7, to E(Y | X) = p = 0,5 + 0,7 = 1,2, czyli p > 1. Do tego nie może dojść, gdyż nie dopuszcza do tego sama konstrukcja wzoru (12). Prawdopodobieństwo po prostu nie może przekraczać 1. Wobec tego niemożliwe jest zwyczajne dodanie składnika losowego do Z, aby dostać Y. Oznacza to, że Y będzie funkcją nieliniową i poszukiwana funkcja g będzie postaci:

W sumie możemy jedynie napisać:

I tak jak stwierdziliśmy wyżej, w przypadku tylko jednej zmiennej niezależnej X - ale tylko zero-jedynkowej - prawdziwe będzie równanie:

co oznacza, że dla funkcji jednej zmiennej regresja logistyczna pokrywa się z modelem liniowym prawdopodobieństwa. Jaka jest tego przyczyna? Po prostu X ma tylko dwa punkty, a więc Y też musi mieć! Skacze więc z jednego punktu do drugiego i stąd sztucznie zachowuje się jak funkcja liniowa.
Dla zmiennej X nie-zerojedynkowej już tak nie będzie. 

Pozostało zadanie estymacji parametrów. Chociaż doszliśmy do dość zadziwiającego liniowego wyniku dla 1 zmiennej X, który prowadzi do wniosku, że można parametry poszukać MNK, to w ogólnym przypadku korzysta się z metody największej wiarygodności. Poszukuje się w niej w takich parametrów, aby prawdopodobieństwo obserwowanych danych (zmiennych o rozkładzie logistycznym) było największe. Na szczęście to trudne zadanie wykonuje program.


Przykład 3. 
Wykorzystamy te sam zestaw 50 danych co w przykładzie 1: X - Czy kierowca ma zdolności? 1 - Tak, 0 - Nie. Y - Czy brak wypadku? 1 - Tak, 0 - Nie:

Wchodzimy w gretla. Model -> Ograniczona zmienna zależna -> model logitowy -> binarny. Zaznaczamy jako zmienną niezależną Czy_kierowca_ma_zdolnosci, a jako zależną - Czy_bez_wypadku. Dodatkowo zaznaczamy opcję Pokaż wartość p (zamiast Pokaż efekt krańcowy dla średniej). Dostałem taki wynik:

Model okazuje się istotny statystycznie. Poniżej widzimy macierz błędów: Dane przewidywane (Y) kontra empiryczne (X): wiersz drugi, kolumna 2, czyli 1x1, oznacza ile razy model trafnie przewidywał brak wypadku (18). Mówiąc inaczej, ile razy występowało Y = 1 oraz X = 1. Wiersz pierwszy, kolumna pierwsza, czyli 0x0, oznacza ile razy model trafnie przewidywał wypadek (14). Mówiąc inaczej, ile razy występowało Y = 0 oraz X = 0. Pierwszy wiersz, druga kolumna, czyli 0x1 pokazuje ile razy model błędnie przewidywał, że kierowca przejedzie bez wypadku (12). Drugi wiersz, pierwsza kolumna, czyli 1x0 pokazuje ile razy model błędnie przewidywał, że kierowca będzie miał wypadek (6). Jak widać najwięcej błędów model popełnił przewidując brak wypadku, gdy kierowca miał zdolności, natomiast stosunkowo niedużo błędów popełnił szacując wystąpienie wypadku, gdy kierowca nie miał zdolności. Jest to sygnał, że słabe zdolności kierowcy rzeczywiście będą prognozować wypadek, natomiast wysokie zdolności mogą nie być wystarczającym czynnikiem do oceny czy nastąpi wypadek.

Zanim poprawimy model, sprawdźmy w końcu to, do czego od początku dążyliśmy - jakie jest prawdopodobieństwo Y = 1. Z uzyskanego okna wyników wchodzimy w Prognozę. Otrzymałem:

A zatem tak samo jak w liniowym modelu, każda obserwacja generuje prawdopodobieństwo wypadku, p(i, t). Jest to możliwe, bo wykorzystujemy teoretyczny poziom prawdopodobieństwa, dla którego, jak wyżej napisałem, X przyjmuje konkretne wartości dla t-tej obserwacji. Wg modelu, jeżeli kierowca ma zdolności, to szansa na brak wypadku, wynosi tylko 60%. A jeśli nie ma zdolności, to 30%. Czyli dostajemy to samo co w liniowym modelu, co wyżej sami przewidzieliśmy.

Przykład 4.
Spróbujmy model poprawić dodając nową zmienną niezależną X2, tę samą co w przykładzie 2: czy inni kierowcy mają zdolności? Tak - 1, Nie - 0. Tabela wszystkich zmiennych z wartościami jest teraz taka:

Wszystko robimy tak samo jak poprzednio, tylko dodajemy nową zmienną niezależną do modelu. Wynik dla logit:

Wszystkie kryteria, od BIC do AIC się poprawiły, czyli model się poprawił. Widać to zresztą po macierzy błędów: sumarycznie wcześniej było 18 błędów (12+6), teraz 14 (6+8).  Nowa zmienna okazuje się bardzo istotna stat. Teraz prognoza:

Zapiszmy prognozę w bazie danych (znak + w oknie prognozy), np. o nazwie Pr_Czy_bez_wypadku_log2. Wszystkie kolejne kroki możemy powtórzyć z przykładu 2. Czyli znów, chcąc prognozować wypadek, mając nowe informacje o X1 i X2, dopisujemy je do obserwacji i sprawdzamy jaka będzie prognoza Y.

Widzieliśmy, że dla X o 1 zmiennej prognozy w modelu liniowym i logicie były identyczne. Porównajmy teraz prognozy z przykładu 2 i przykładu 4:

Okazuje się, że prognozy są bardzo podobne, choć tym razem nieco się różnią.

Na koniec zaznaczę, że zarówno liniowy model prawdopodobieństwa jak i regresja logistyczna nie muszą być stosowane tylko dla zero-jedynkowej zmiennej X. To Y musi być zawsze binarny. Więcej o tym, z przykładami, napiszę następnym razem.

Od zasady maksymalnej entropii do regresji logistycznej

Wiemy już skąd się bierze idea entropii i zasada maksymalnej entropii (ZME). Jeśli system jest "sprawiedliwy", tak że prawdopodobieństwo każdego zdarzenia jest identyczne, to w takim systemie entropia będzie zawsze najwyższa. Podobnie, jeśli nie mamy żadnych informacji o systemie i z góry zakładamy rozkład jednostajny, to w takim systemie entropia będzie najwyższa. To twierdzenie można odwrócić: maksymalna entropia będzie prowadzić do rozkładu jednostajnego, jeśli nie mamy żadnych informacji o nim (może poza tym, że zdarzenia są niezależne od siebie). Gdy pojawią się nowe informacje, ZME doprowadzi do powstania innych rozkładów.

Powiedzmy, że rzucamy sześcienną kostką do gry, ale nie wiemy czy jest sprawiedliwa, a więc nie mamy pewności czy szansa każdego zdarzenia wynosi 1/6. Mimo to nie posiadamy żadnych innych informacji o tej kostce oprócz tego, że pociąga 6 różnych zdarzeń. Wspomaga nas jeszcze jedna informacja, ale spoza układu, wynikająca z samej teorii: prawdopodobieństwa sumują się do 1. Dostajemy dwa równania:

gdzie H to entropia Shannona, a p(i) - prawdopodobieństwo zdarzenia x(i).

Ten układ można rozwiązać metodą Lagrange'a. ZME zostaje zapisana jako zadanie optymalizacyjne:

Rozwiązaniem tego zadania będzie p = 1/6.

Ponieważ standardowa maksymalizacja oparta na różniczkowaniu wymaga ciągłości funkcji, to należałoby właściwie zamienić sumę na całkę. Dlatego w przypadku dyskretnym najłatwiej użyć nierówności Jensena do dowodu.

Świat rynków finansowych znacznie różni się od rzutu kostką. Niektóre zdarzenia są bardziej, inne mniej prawdopodobne. Najczęstsze zdarzenia układają się wokół dominanty, a jeśli rozkład jest symetryczny, to dominanta równa się wartości oczekiwanej. Załóżmy, że kostka jest niesprawiedliwa, tak że premiuje niektóre zdarzenia. Nie wiemy jednak czy rozkład jest symetryczny czy nie. Dostajemy tylko informację o wartości oczekiwanej (choćby na podstawie wielu rzutów). Pamiętajmy, że wartość oczekiwana nie musi leżeć na środku ani być najbardziej prawdopodobną wartością. Otrzymujemy 3 równania:

gdzie x(i) to wartość i-tego zdarzenia, a x z kreską na górze to wartość oczekiwana.
Powstaje zadanie optymalizacyjne:

Rozwiązaniem tego zadania będzie:

(1)

Jest to rozkład Boltzmanna, znany z fizyki. W zależności od parametru lambda(2) i wartości oczekiwanej, otrzymamy różne wyniki. Np. jeśli wartość oczekiwana = 3,5, to dostaniemy z powrotem rozkład jednostajny:

 
Gdy wartość oczekiwana = 2,5 i lambda(2) = 0,371, to rozkład przybiera postać jak na poniższym rysunku:

Gdy zaś wartość oczekiwana = 4,5 i lambda(2) = -0,371, wtedy rozkład przybiera postać:

Widać, że pomimo znajomości wartości oczekiwanej, w dwóch ostatnich przykładach nie ma żadnej symetrii ani nawet zalążka dzwonu (dominanta to odpowiednio 1 i 6). Mówię o tym, bo musimy sobie uświadomić, że sama ZME przy warunku znajomości jedynie wartości oczekiwanej, nie doprowadzi do rozkładu typowego dla rynków (a tym bardziej do Gaussa).

Natomiast uzyskany rozkład prowadzi do interesującego odkrycia. Wzór (1) przekształcamy:

Uzyskany wynik można uprościć. Zauważmy, że x(6) jest zwykłą stałą, a nie zmienną. Zastąpmy stałe nowymi parametrami:

(2)

Licznik w (2) jest to tzw. szansa (ang. odds), a więc prawdopodobieństwo jest to szansa podzielona przez sumę wszystkich szans (na jedynkę w mianowniku możemy spojrzeć jako na wartość wynikającą z redundancji ostatniej - szóstej wartości, której szansa całkowicie zależy od wszystkich poprzednich).
 
W ten sposób zbliżamy się do koncepcji regresji logistycznej (logistic regression), inaczej logit. Wzór (2) przedstawia ogólniejszą jej formę. Gdybyśmy chcieli zapisać ściślej tę funkcję, powinniśmy uwzględnić możliwość występowania wielu zmiennych, jak to często ma miejsce w regresji liniowej, a także N zamiast 6-ciu zdarzeń. A więc np.  dla dowolnego N dostajemy postać tzw. wielomianowej regresji logistycznej (multinomial logistic regression) czy wielomianowego logitu:

(3)

Jak można byłoby sobie wyobrazić tego typu zdarzenie? Np. rzucamy już nie sześcienną ale N-ścienną kostką, w którą uderzają inne N-ścienne kostki. Te inne kostki są od siebie niezależne. Wszystkich kostek jest K: K-1 obcych i jedna nasza.* Formuła (3) wskaże jakie jest prawdopodobieństwo uzyskania i-tego zdarzenia w takim układzie.

Robi się to trochę skomplikowane, dlatego zajmijmy się klasycznym, dwumianowym logitem. Kiedy on zachodzi? W logicie występują tylko 2 zdarzenia: orzeł-reszka albo prawda-fałsz. Domyślamy się więc, że N = 2. Wtedy (3) przybiera postać klasycznego logitu:

(4)

Ponieważ nie ma tu już sumy do N-1, to j = i.
Powiedzmy, że K = 2. Wtedy logit wygląda całkiem przyjaźnie:

Dostajemy model, w którym są tylko dwie zmienne zero-jedynkowe. Dzięki niemu możemy np. sprawdzić jakie jest prawdopodobieństwo, że firma zbankrutuje (albo na odwrót, że zysk wzrośnie), jeżeli spełnione lub niespełnione zostaną jakieś dwa warunki reprezentowane przez dwie zmienne.

Jeżeli dotychczas mieliśmy opór przed użyciem logitu, to po tym tekście powinien zniknąć. Okazuje się być modelem naturalnym i logicznym. Odrębna sprawa to praktyczne jego zastosowanie i tym zajmę się niebawem.


Literatura:
[1] He, F., Maximum entropy, logistic regression, and species abundance, April 2010
[2] https://sgfin.github.io/2017/03/16/Deriving-probability-distributions-using-the-Principle-of-Maximum-Entropy/
[3] https://en.wikipedia.org/wiki/Multinomial_logistic_regression


*Jest to abstrakcyjna gra, więc lepszym przykładem byłaby zamiana kostek na samochody ścigające się ze sobą po trudnej drodze i N=3 zdarzenia: (a) że nasz samochód przejedzie bez szwanku,  (b) że nasz samochód zderzy się z innym samochodem i (c) że nasz samochód będzie miał wypadek. Zmienna x1 oznaczałaby poziom umiejętności kierowcy, również z trzema możliwościami: słabym, średnim i dobrym. Reszta zmiennych, x2, x3, ..., oznaczałaby pozostałe auta, a dokładniej poziomy umiejętności pozostałych kierowców (tak samo: słaby, średni lub dobry). Zmienna x1 dotyczy więc naszych wewnętrznych możliwości, a inne zmienne zewnętrznych możliwości. Ponieważ poziom słaby, średni i dobry będzie silnie korelował odpowiednio z ryzykiem wypadku, zderzenia i braku wydarzeń, to przypomina to kostkę, która bez czynników zewnętrznych generuje kilka prawdopodobnych zdarzeń. Trzeba jednak zauważyć, że oba przykłady subtelnie się różnią. Usuńmy "zmienne zewnętrzne". Dla przykładu z kostką otrzymamy tylko N wylosowanych zdarzeń, natomiast w przykładzie z autami występuje w pewnym sensie 2*N zdarzeń, czyli tutaj 2*3 = 6 zdarzeń: brak wydarzeń, zderzenie, wypadek oraz poziom dobry, średni, słaby naszego kierowcy.  Tak więc w tym ostatnim mamy konkretny podział na zmienną zależną i niezależną; w przykładzie z pojedynczą kostką tak nie było. Gdyby korelacja między zmienną zależną a niezależną wyniosła 1, to moglibyśmy dopiero utożsamić obydwie zmienne ze sobą. I rzeczywiście tak można by zrobić, gdyby nie występowały inne zakłócające czynniki, jak np. poziom koncentracji kierowcy, który staje się de facto nową zmienną modelu (np. z podziałem na słabą, średnią i dobrą koncentrację). Niemniej, gdyby założyć, że koncentracja jest zawsze dobra, to można założyć, że zawsze dobry kierowca w trudnych warunkach nie będzie miał ani wypadku ani zderzenia, średni będzie miał zawsze zderzenie, a nie będzie miał wypadku, a słaby zawsze będzie miał wypadek. W ten sposób zmienna zależna staje się tożsama z niezależną - wtedy dokładnie dostaniemy przykład z kostką. Wówczas możemy dołożyć kolejne zmienne niezależne, co spowoduje złamanie tej doskonałej korelacji między zmienną zależną a niezależną, czyli wracamy do typowego modelu regresji.

Szybkie wyprowadzenie entropii informacyjnej Shannona

Entropia informacyjna jest ważnym pojęciem w statystyce i ekonometrii, ale raczej mało znanym w finansach. Czasami wykorzystuje się ją jako miarę ryzyka inwestycyjnego, lecz zanim jej użyjemy powinniśmy zrozumieć jej istotę. Spójrzmy najpierw na ten prosty wzór na entropię Shannona:

(1)

gdzie p(i) to prawdopodobieństwo zajścia pewnego zdarzenia x(i).

Pierwsze pytanie jakie przychodzi na myśl to: dlaczego w ogóle taki wzór? W polskich, a nawet zagranicznych opracowaniach internetowych nie znajdziemy szybkiej odpowiedzi. Oczywiście warto zajrzeć do oryginalnej pracy Shannona [2], która szczegółowo wyjaśnia użycie takiej a nie innej postaci, jednak dla laika w teorii informacji praca ta jest ciążka do przetrawienia i brakuje mi tam czysto matematycznego wyprowadzenia bez pojęcia bitów. Dlatego zacząłbym od prac, które wyprowadzają entropię wprost z rozkładów prawdopodobieństwa. Np. Hanel et al. [1] zaczynają od rozkładu dwumianowego, w którym występują tylko dwa zdarzenia: wyrzucenie orła albo reszki. Z tym że nie mamy pewności ani założenia, że rzucana moneta jest symetryczna i że oba zdarzenia są równie prawdopodobne, tj. nie muszą wynosić po 50%. Taką monetą rzucamy N razy (i każdy rzut jest niezależny od siebie). Powiedzmy, że k1 wskazuje ile razy wyrzucimy orła, a k2 ile razy wyrzucimy reszkę (czyli N = k1 + k2). Prawdopodobieństwo wyrzucenia k1 razy orła wynosi:

(2) 

gdzie p(1) to prawdopodobieństwo wyrzucenia raz orła, a p(2) to prawdopodobieństwo wyrzucenia raz reszki.

Jak wyżej powiedziano nie znamy teoretycznego p(i), jednak przy dużym N empiryczne prawdopodobieństwo będzie się zbliżać do teoretycznego. Dlatego uznajemy, że p(i) = k(i) / N. Dla dużego N można aproksymować wzór (2), korzystając z formuły Stirlinga. Wtedy (2) można zapisać w postaci:

(2)

gdzie H = -p(1)*ln(p1) - p(2)*ln(p2) = -[p(1)*ln(p1) + p(2)*ln(p2)].

Czyli widzimy, że:

(3)

Dlaczego entropia stanowi miarę nieporządku? Wiadomo, że N! / (k1! / (N - k1)!) to liczba kombinacji (bez powtórzeń) ułożenia k1 orłów w N rzutach. Czyli przy stałym N to właśnie entropia H mierzy liczbę tych kombinacji. Natomiast to, że im więcej kombinacji tym większy nieporządek, jest czymś intuicyjnie zrozumiałym. Oznacza to więcej możliwych niezależnych od siebie zdarzeń (trzeba zwrócić uwagę, że gdy zdarzenia są od siebie zależne, to liczba kombinacji nie będzie już oznaczać nieporządku, ale raczej mierzyć ilość informacji i wtedy taką definicję może przyjąć entropia - zob. ten wpis dotyczący termodynamiki nierównowagowej).

Teraz kolejne spostrzeżenie pozwoli nam zrozumieć znaczenie entropii. Powiedzmy, że wiemy, że moneta jest sprawiedliwa, czyli otrzymanie orła jest tak samo prawdopodobne jak reszki i musi wynieść 0,5. Podstawmy teraz p(1) = p(2) = 0,5 do wzoru (1):

Następnie przypuśćmy, że moneta nie byłaby sprawiedliwa, tak że p(1) = 0,6, a p(2) = 0,4. Wtedy H równa się:

W trzecim przykładzie niech p(1) = 0,8 i p(2) = 0,2:

W ostatnim przykładzie dla p(1) = 1 i p(2) = 0:

Wniosek nasuwa się sam: dla każdego innego prawdopodobieństwa niż 0,5 entropia będzie zawsze mniejsza. Czyli szansa wyrzucenia orła lub reszki równa 50% zapewnia najwyższą entropię. 

Co w sytuacji, gdy zdarzeń jest więcej niż 2? Dostajemy wtedy rozkład wielomianowy i prawdopodobieństwo wyrzucenia k(i) razy zdarzenia x(i) wynosi:

(4)

Wzór (4) możemy zapisać analogicznie jak (2):

(5)

gdzie H to właśnie wzór nr (1).

Idąc tym samym tropem co poprzednio zauważamy, że dla stałego N entropia H będzie warunkować liczbę kombinacji i tym samym mierzyć poziom nieuporządkowania układu.

Jeżeli układ ma pozostać "sprawiedliwy" dla każdego zdarzenia, to musimy przyjąć rozkład jednostajny, gdzie k1 = k2 = ... = kW = k. Wtedy p(i) = k / N, tzn. każde zdarzenie jest tak samo prawdopodobne. Domyślamy się już, że równa szansa każdego zdarzenia prowadzi do osiągnięcia maksimum entropii.

Powyższe wyprowadzenie stanowi źródło drugiej zasady termodynamiki, zgodnie z którą w układzie izolowanym entropia nie maleje. Gdy dodamy do tego izotropowość i jednorodność przestrzeni i czasu (np. symetryczność monety), wtedy entropia zawsze będzie maksymalna. 

Stąd rodzi się tzw. zasada maksimum entropii - ważna koncepcja w statystyce. O niej jednak więcej opowiem następnym razem.


Literatura:
[1] Hanel, R., Thurner, S., Gell-Mann, M., How multiplicity determines entropy and the derivation of the maximum entropy principle for complex systems, May 13, 2014,
[2] Shannon, C. E., A Mathematical Theory of Communication, July, 1948.

Filtr Butterwortha + ARMA, czyli w górę albo w dół

Od kilku lat jestem negatywnie nastawiony do wzrostów na giełdzie, chociaż ostatnio w związku z zastosowaniem filtru Butterwortha (FB), pisałem, że WIG jest coraz bliżej dołka. Koniec FB sugeruje dołek na indeksie - tym razem wziąłem pod lupę WIG20 (kwartalne stopy od stycznia 2002 do marca 2019 r.):

Parametry: x_c = 27, n = 2. FB możemy teraz prognozować już klasyczną ekonometrią - tj. ARMA. Starujemy od wyrównanego szeregu FB:

Pierwsze różnice FB można uznać za stacjonarne (ADF i KPSS przemawiają za nią). Dodatek armax w gretlu wskazał optymalne parametry ARMA(4, 2) dla BIC i HQC jednocześnie. Możemy teraz użyć wskazanych pierwszych różnic FB do ARMA(4, 2) albo od razu zastosować oryginalny FB do ARIMA(4, 1, 2). Ta ostatnia pozwala graficznie ocenić kierunek FB. I tak otrzymałem prognozę na 4 kolejne kwartały:

Same przyrosty FB, czyli ARMA(4, 2), daje podobnie pozytywny sygnał:

Z powyższego wynikałoby, że mój negatywny sentyment do rynku jest błędny, a giełda chce się wzbijać. Jednakże problem z takimi modelami jak ARIMA jest tego typu, że prognoza silnie zależy od zakresu obserwacji. Pytanie brzmi czy zostałaby ona zachowana, gdyby obserwacje zacząć wcześniej, np. od stycznia 2000 r.? Okazuje się, że wtedy się zmienia. Optymalny model będzie to ARIMA(2, 1, 9), a prognoza:

 

Gdyby zrobić analizę na zmianach FB, czyli ARMA(2, 9), to jeszcze bardziej uwidacznia się ryzyko spadków:

A co, gdyby wziąć jeszcze wcześniejszy okres, np. od stycznia 1995? Wtedy optimum staje się ARIMA(2, 1, 6), a prognoza przedstawia się następująco:

Czyli znowu prognoza trendu rosnącego. Różnica jednak między pierwszą a ostatnią jest taka, że ta ostatnia jest mniej pewna wzrostów: sugeruje bardzo słaby zwyżkowy trend liniowy FB. W takim razie porównajmy to jeszcze z prognozą zmian, tj. użyjmy ARMA(2, 6):

  

Potwierdzona zostaje "prognoza niepewności i wątpliwości". Pamiętać trzeba tylko, że FB wykonywany był na stopach zwrotu, a więc ARMA na zmianach FB wskazuje kierunek zmian stóp zwrotu, a nie indeksu.