piątek, 12 sierpnia 2022

Oczyszczanie zysku netto z niepewnych pozycji

Ostatni wpis skłonił mnie do stworzenia formuły na oczyszczony zysk netto z niepewnych lub uznaniowych elementów, jak zmiany w aktywach rzeczowych. Owszem, mają one jakąś wartość, ale realne korzyści są odłożone w niepewną przyszłość, a w dodatku jeżeli ich wycena została dokonana przez samą firmę, to istnieje ryzyko, że nie odzwierciedla ona faktycznej wartości wewnętrznej tego aktywa. W sprawozdaniu finansowym PKN Orlen można znaleźć metodologię wyceny odpisów aktualizujących wpływającą na wynik finansowy, jednak po pobieżnym przeczytaniu dostrzegłem kilka założeń za nią stojących. Nie chcę ich tu przywoływać, ale żeby sprawdzić ich poprawność, trzeba by dokładnie je prześledzić. W dodatku powstawanie różnicy między sumą zysków z 4 kwartałów a zyskiem za cały rok przez kilka lat rzędu pozwala powątpiewać w słuszność uznawania zmian w aktywach trwałych za coś realnego.

Przypomnę teraz wzór, jaki pokazałem w tym art. (był to wzór nr 5):

(1)

Czyli zmiana stanu gotówki (w ogólnym sensie) = zysk netto + amortyzacja - zmiana stanu zapasów - zmiana należności + zmiana zobowiązań (krótkoterminowych i długoterminowych) - zmiana aktywów trwałych (wraz ze zmianą ich amortyzacji)  - wypłacone dywidendy + emisja akcji.

Przedstawione równanie wskazuje relację między zyskiem księgowym a pieniężnym. Aby zysk netto stał się gotówką, trzeba odjąć lub dodać  wymienione pozycje.

Na wiele z tych pozycji możemy spojrzeć przez pryzmat ryzyka. Po prawej stronie leżą środki pieniężne, których ryzyko otrzymania jest zerowe, bo mamy je na koncie. Ryzyko bliskie zera dotyczy amortyzacji, która odzwierciedla zużywanie się środków trwałych. Ich zakup odbył się w przeszłości i ten koszt został rozbity na wiele okresów do przodu. 

Omówimy teraz po kolei te elementy.

Amortyzacja

W momencie zakupu środków trwałych ponosimy na nie wydatki, ale ponieważ środki te nie są w całości wykorzystane do produkcji, ale stopniowo, to nie stanowią kosztu uzyskania przychodów. Dlatego owe wydatki są rozpisane w formie kosztu poprzez podział na wiele okresów do przodu - stąd powstaje amortyzacja w kolejnych miesiącach czy latach. 

W tym momencie można zauważyć różnicę między spojrzeniem księgowego a inwestora: księgowego interesuje średni zysk firmy ze wszystkich okresów, z kolei inwestora zysk firmy z ostatniego okresu. Inaczej mówiąc księgowego interesują wartości średnie, inwestora - wartości krańcowe. Wynika to z faktu, że księgowy zajmuje się przeszłą oraz teraźniejszą sytuacją finansową, natomiast inwestor przyszłą i korzysta do tego z bieżących danych.

Mimo tej różnicy obie perspektywy się zazębiają. Zakupiony środek trwały pracuje przez lata, ale też fizycznie i ekonomicznie się amortyzuje, co właśnie odzwierciedla księgowa amortyzacja. Jeżeli maszyna zaczyna się np. zacinać, to jest to element tego kosztu. Co więcej, jeśli można by ją zastąpić nowocześniejszą wersją, to także mamy do czynienia z tym kosztem (zob. np. art. 32, p. 2 Ustawy o Rachunkowości) Jeżeli przestaje być ona zdatna do użytku i trzeba ją wymienić, to środek zostaje fizycznie zamortyzowany. Jeżeli firma jest zmuszona do zakupu nowej maszyny, to znowu nie zobaczymy tego wydatku w obniżonym zysku, ale jednak gotówka się zmniejszyła. Zatem amortyzację można traktować dwojako: 

(a) jako koszt fizycznego i ekonomicznego zużycia lub

(b) jako rodzaj kosztu rezerwy albo średni wydatek ważony ryzykiem nowych wydatków.

Powyższe dowodzi, że powszechne twierdzenie, iż firma może wypłacić dywidendę nie tylko z zysku netto, ale też z amortyzacji, jest fałszywe. Chociaż w pewnych okresach będzie prawdziwe, to nie raz się zdarzy, że firma nie będzie mogła wypłacić (dużej) dywidendy, bo albo musiała wydać pieniądze na inwestycje w środki trwałe, albo przewiduje, że za chwilę będzie musiała je wydać. 

Niemniej problemem pozostaje ustalenie okresu amortyzacji i stawek amortyzacyjnych, które mogą rozmijać się z rzeczywistą użytecznością aktywów. To tutaj leży element niepewności, który inwestor powinien rozważyć, tzn. czy dodać do zysku amortyzację, jeśli tak, czy całą, czy tylko część i jaką. 

Możemy przyjąć założenie, że jeśli wartość środków trwałych rośnie o więcej niż zobowiązania finansowe oraz emisja akcji, to nie dodajemy w ogóle amortyzacji, a jeśli o mniej - to dodajemy w zależności od naszego poziomu precyzji - albo całą amortyzację w uproszczeniu, albo taką część, która pozostanie jeszcze po zapłacie za nowe środki trwałe. W tym drugim przypadku dodamy kwotę równą:

amortyzacja - min(amortyzacja, x), 

gdzie x to różnica między wzrostem wartości środków trwałych a wzrostem zobowiązań finansowych i nowych emisji akcji. 

W tym miejscu dość istotna uwaga. We wzorze (1) występuje określenie gross fixed assets oznaczające aktywa trwałe brutto, czyli aktywa trwałe netto powiększone o ich amortyzację. Ta amortyzacja występuje właśnie dlatego, że może finansować owe środki. Ale jest ona odejmowana, wraz ze zmianą środków trwałych, więc całkowita amortyzacja ulega pomniejszeniu. Możemy tę trudność ominąć. Mamy do czynienia ze zmianą tych środków, czyli:

środki netto z okresu 1 + amortyzacja tych środków z okresu 1 - środki netto z okresu 0 - amortyzacja tych środków z okresu 0.

Stąd:

(środki netto z okresu 1 - środki netto z okresu 0) + (amortyzacja środków brutto z okresu 1 - amortyzacja środków brutto z okresu 0).

Inaczej:

zmiana środków trwałych netto + zmiana amortyzacji środków trwałych brutto.

Jeżeli amortyzacja w okresie bieżącym jest zbliżona do amortyzacji z okresu poprzedniego, to zmiana jest zerowa. Przyjmiemy to założenie, które wyeliminuje konieczność korygowania amortyzacji. Innymi słowy będziemy posługiwać się aktywami trwałymi netto, zamiast brutto.

Zmiana stanu zapasów

Często niezrozumiałe jest dla wielu, dlaczego nowe zapasy są ujmowane jako część przychodów, skoro są... no właśnie tylko zapasami - nie zostały sprzedane. Wydaje się to bez sensu. Aby zrozumieć ten mechanizm, trzeba poznać zasady księgowości, jak kontynuacja działalności oraz współmierność kosztów do przychodów. Jeżeli w danym okresie firma kupuje zapasy, zwiększa poziom aktywów, które posłużą do produkcji, ale produkty gotowe będą sprzedane w przyszłości. Z powodu zasady współmierności kosztów i przychodów trzeba więc odjąć te wydatki z kosztów. Czyli koszty się zmniejszają. Ale zmniejszenie kosztów jest równoznaczne ze zwiększeniem przychodów o tyle samo. Stąd zysk netto zawiera sztucznie zmianę stanu zapasów na plus. Ze względu na zasadę kontynuacji działalności, przewidujemy, że w niedalekiej przyszłości (skoro to aktywa krótkoterminowe, to znaczy, że do 1 roku) zapasy przekształcą się w produkty gotowe i faktycznie zostaną sprzedane, generując realny cash.

Jak by nie patrzeć, to zmiana zapasów stanowi czynnik niepewności wyniku. Są branże bardziej pewne i mniej pewne, są firmy, które mają już ugruntowaną pozycję rynkową, ale też takie, które dopiero wchodzą na rynek. W księgach tego nie ujrzymy. Jeżeli mamy pewność, że firma sprzeda towary, to możemy pominąć zmianę zapasów, a jeśli nie mamy pewności, to lepiej odjąć. Analiza ta może być ciekawym zadaniem intelektualnym.

Jeżeli przyjmujemy brak pewności, to odejmujemy zapasy od zysku, czyli zostawiamy je po lewej stronie równania (1). Jeżeli wierzymy, że za chwilę wygenerują zysk, zostawiamy zapasy w zysku netto - wtedy przenosimy na prawą stronę.

Oczywiście możemy przyjąć rozwiązanie pośrednie i wybrać pewne prawdopodobieństwo przejścia zapasów z produkcji w realną sprzedaż. W ogólności wartość oczekiwana sprzedanych nowych zapasów byłaby równa:

(1-p1)*zmiana stanu zapasów , 

gdzie p1 - prawdopodobieństwo przejścia zapasów z produkcji w realną sprzedaż.

Do obliczeń można wspomóc się cyklem zapasów lub jego odwrotnością, tj. rotacją zapasów. O tym więcej przy omawianiu przykładów.

Zmiana należności

Należności w przeciwieństwie do zapasów są pozycją w pewną w tym sensie, że nasz kontrahent musi nam zapłacić za już sprzedane towary lub usługi. Oczywiście zawsze jest ryzyko, że nie zapłaci, ale dużo mniejsze niż to, że nie sprzedamy towarów, jak to było z zapasami. 

Dlatego należności zawsze przenosimy na prawą stronę.

Zmiana zobowiązań

To samo co z należnościami, tylko że to nasza firma musi zapłacić kontrahentowi. Czyli przenosimy zobowiązania na prawą stronę.

Zmiana aktywów trwałych

Aktywa trwałe mają na celu przynieść korzyści finansowe, jednak po pierwsze nie jest pewne ani to, że przyniosą, ani to, kiedy przyniosą, ani to, jak dużo przyniosą. Dlatego nie będziemy mieć oporów, aby pozostawić aktywa trwałe po lewej stronie (zostaną więc odjęte od zysku). Są jednak branże, w których aktywa trwałe mają szczególne znaczenie - budownictwo i nieruchomości. W ich przypadku aktywa trwałe stanowią dużą część wartości firmy, stąd odejmowanie całej tej pozycji jest niesprawiedliwe. Dla uproszczenia w przypadku takich branż potraktujemy aktywa trwałe jak zapasy, a więc przemnożymy zmianę aktywów trwałych przez (1-p1). 

Wypłata dywidend

Dywidendy są nagrodą dla akcjonariuszy, ale dotyczą poprzedniego okresu. Ich wypłata to formalna operacja, która nie ma znaczenia dla działalności tu i teraz. Co więcej, nawet gdyby dywidendy dotyczyły bieżącego okresu, to nie ma to dla nas znaczenia, dlatego że oczyszczanie, jakiego dokonujemy ma na celu właśnie ocenić, czy spółka jest w ogóle w stanie wypłacić dywidendę z zysku. Gdyby po oczyszczeniu miała stratę, to mogłaby wypłacić dywidendę najwyżej z kapitału własnego, a to nas nie interesuje. Dlatego dywidendy przenosimy na prawą stronę, aby nie odejmować ich od zysku.

Emisja akcji

Analogiczna sytuacja jak w przypadku zobowiązań - nowy kapitał akcyjny jest to forma zobowiązań dla akcjonariuszy, a więc przenosimy na prawą stronę.


Oczyszczony zysk netto

Aby zachować równanie (1) dodamy i odejmiemy wszystkie dodatkowe elementy, o których było wyżej, a także przeniesiemy na prawą stronę te elementy, o których byłą mowa:

Zysk netto + amortyzacja - min(amortyzacja, x) + min(amortyzacja, x) - (1-p1)*(zmiana stanu zapasów) - p1*(zmiana stanu zapasów)  - (1-p2)*zmiana aktywów trwałych netto - p2*(zmiana aktywów trwałych netto) 

środki pieniężne + zmiana należności - zmiana zobowiązań (krótkoterminowych i długoterminowych) + wypłacone dywidendy – emisja akcji

Czyli:

(2)

Zysk netto + amortyzacja - min(amortyzacja, x) - (1-p1)*(zmiana stanu zapasów) – (1-p2)*(zmiana aktywów trwałych netto) = środki pieniężne + zmiana należności - zmiana zobowiązań (krótkoterminowych i długoterminowych) + wypłacone dywidendy – emisja akcji - min(amortyzacja, x) + p1*(zmiana stanu zapasów) + p2*(zmiana aktywów trwałych netto)   

Lewa strona (2) będzie stanowić oczyszczony zysk netto:

Oczyszczony zysk netto = Zysk netto + amortyzacja - min(amortyzacja, x) - (1-p1)*(zmiana stanu zapasów) – (1-p2)*(zmiana aktywów trwałych netto)

gdzie 
p1 = min(oczekiwane przychody  ze sprzedaży produktów * wskaźnik operacyjności / (średnie) zapasy ; 1)

gdzie 
wskaźnik operacyjności = koszty wytwarzania / przychody ze sprzedaży produktów.

Jeżeli za oczekiwane przychody przyjmiemy bieżące przychody, wtedy 

p1 = min(rotacja zapasów ; 1)
gdzie:
rotacja zapasów = koszty wytwarzania / zapasy

Definicji rotacji zapasów nie traktujmy zupełnie ściśle. Generalnie "wskazuje" ile razy zapasy są na nowo odtwarzane w ciągu danego okresu. Wynika to z faktu, że zapasy są aktywami  krótkoterminowymi, więc są wykorzystywane do roku czasu i stąd mogą być porównywane z rocznymi kosztami firmy. Koszty wytwarzania będą prawie zawsze większe niż wartość zapasów (która stanowi koszt ich nabycia), ponieważ oprócz części zapasów zawierają m.in. koszty robocizny (pełna definicja tych kosztów znajduje się w art. 28, p. 3 Ustawy o rachunkowości). Są jednak wyjątki, jak branża deweloperska, gdzie rotacja zapasów jest zawsze mniejsza od 1. Od 1 kw 2016 do 1 kw 2022 rotacja wahała się od 0,5 do 1. Budynki buduje się zwykle kilka lat, a że to droga inwestycja, to nie sprzedaje się z dnia na dzień.

Często też spotykana jest rotacja zapasów w postaci:

rotacja zapasów = koszty wytwarzania / średnia zapasów z 4 ostatnich kwartałów.

x = zmiana aktywów trwałych netto - zmiana zobowiązań finansowych - emisja nowych akcji ;

p2 = p1 dla branży budowlanej lub nieruchomości, a dla pozostałych 0 lub inna subiektywna wartość.

Przykłady

a) PKN Orlen. Założenia: p1 = 1, p2 = 0. 
Wyniki 2021: źródło: link
Zysk netto = 11 188
amortyzacja = 5 341
aktywa trwałe netto 2021 = 68 706
aktywa trwałe netto 2020 = 59 433
zmiana aktywów trwałych =  68 706 - 59 433 = 9 273
zobowiązania finansowe 2021:
    długoterminowe: kredyty, pożyczki, leasing i obligacje = 18 618
    krótkoterminowe: kredyty, pożyczki, leasing i obligacje = 2 108
zobowiązania finansowe 2020:
    długoterminowe: kredyty, pożyczki i obligacje = 13 931
    krótkoterminowe: kredyty, pożyczki i obligacje = 5 643

zmiana zobowiązań fin = (18 618 + 2 108) - (13 931 + 5 643) = 1 152

x = 9 273 - 1 152 = 8 121

Oczyszczony zysk = 11 188 + 5 341 - min(5 341; 8 121) – 9 273 = 1915

O ile księgowy zysk netto wynosił ponad 11 mld zł, o tyle "faktyczny", oczyszczony ze zmian aktywów trwałych zysk wyniósł 1,9 mld zł, czyli tylko 17% pierwotnego. Okazuje się więc, że niewielka dywidenda 3,5 zł stanowi prawie 80% całego zysku (3,5 / 26 = 13,5% => 13,5 / 17 = 79%). Zwróćmy jeszcze uwagę, że cena / oczyszczony zysk dawałoby 15-16.


b) Kolejnym przykładem będzie Cavatina, spółka deweloperska. Podobnie jak PKN ma niskie C/Z oraz C/WK, obecnie odpowiednio 3,35 i 0,5, co oznacza ROE = 15% , a więc powyżej średniego kosztu kapitału własnego. Podobnie jak w przypadku PKN wydawałoby się, że jest niedowartościowana. Aby to ocenić, oczyśćmy zysk netto.
 
Źródło: biznesradar.pl (link).
Wyniki 2021:
Założenia: spróbujmy oszacować p1 przyjmując wzór:

p1 = min(rotacja zapasów ; 1)
rotacja zapasów = koszty wytwarzania / zapasy
koszty wytwarzania = 18 201
zapasy 2021 = 64 203
p1 = 18 201 / 64 203 = 0,28

p2 = p1

zapasy 2020 = 32 954
zmiana stanu zapasów = 64 203 - 32 954 = 31 249

Zysk netto = 189 620

amortyzacja = 3 733

aktywa trwałe netto 2021 = 2 150 469
aktywa trwałe netto 2020 = 1 292 025
zmiana aktywów trwałych = 858 444

zobowiązania finansowe 2021:
długoterminowe: kredyty, pożyczki, leasing i obligacje = 1 001 651
krótkoterminowe: kredyty, pożyczki, leasing i obligacje = 180 345
zobowiązania finansowe 2020:
długoterminowe: kredyty, pożyczki, leasing i obligacje = 490 949
krótkoterminowe: kredyty, pożyczki, leasing i obligacje = 74 614
zmiana zobowiązań fin = 616 433

x = 858 444 - 616 433 = 242 011

Oczyszczony zysk = 189 620 + 3 733 - min(3 733; 242 011) - (1-0,28)*31 249 - (1-0,28)*858 444 = -450 959

Komentarz:
Księgowy zysk netto Cavatiny wyniósł w 2021 prawie 190 mln zł, jednak po oczyszczeniu z pozycji niepewnych dostaliśmy stratę w tym sensie, że spółka nie jest w stanie wypłacić dywidendy z zysku. Potwierdzają to też przepływy z działalności oper., które są ujemne.

Dwie modyfikacje:
Biorąc pod uwagę, że w działalności deweloperskiej zapasy wolno schodzą, moglibyśmy posłużyć się wcześniej wspomnianym wzorem na rotację zapasów:

rotacja zapasów = koszty wytwarzania / średnia zapasów z 4 ostatnich kwartałów.

Taki wskaźnik oblicza biznesradar.pl , więc nie musimy się trudzić i to wykorzystamy. Wskaźnik na 2021 rok wyniósł 0,36, więc taki wstawiamy za p1 i p2:

Oczyszczony zysk = 189 620 + 3 733 - min(3 733; 242 011) - (1-0,36)*31 249 - (1-0,36)*858 444 = -379 784

Bardzo sympatyczne ze strony portalu, że podaje od razu wartości wskaźników w sektorze. Dla 2021 sektor miał rotację ponad 2 razy wyższą, bo 0,8. Co ciekawe, gdyby podstawić ten wskaźnik pod p1 i p2, oczyszczony zysk stałby się lekko dodatni. W dotychczasowej historii spółka osiągnęła raz taki poziom w 2020 r. Jeżeli wierzymy, że nastąpi regresja do średniej, to możemy wstawić średnią historyczną, która wynosi 0,72, wtedy zbliżymy się prawie do zera. 

Na 2021 C/Z wyniosło 2,9, a C/WK 0,52, co w tym momencie przestaje być niezrozumiałe czy nieracjonalne. Podobnie dziś nie ma sensu inwestować w takie akcje. Generalnie inwestorzy powinni rozważyć prędzej inwestowanie w jej obligacje, których rentowność znacznie przekracza 10% (zdaje się, że obecnie daje ok. 13%).  

środa, 3 sierpnia 2022

Kreatywna księgowość czy cudotwórca? Może ktoś wie?

Przed chwilą porównałem wyniki finansowe PKN Orlen w różnych okresach i zauważyłem dziwny błąd w sprawozdaniach w latach 2018-2021. Wchodzę najpierw na rok 2021:

https://www.orlen.pl/pl/relacje-inwestorskie/raporty-i-publikacje/sprawozdania/2021

i otwieram plik "Dane liczbowe ze sprawozdania" za cały rok. i wybieram dane skonsolidowane. W kom. D20 jest podany zysk netto  11 188

Następnie otwieram plik za 4 kwartał i w kom. D26 (12 MIESIĘCY ZAKOŃCZONE 31/12/2021) podano zysk netto 10 241. 

Skąd różnica 947 mln ? Normalnie roczny wynik powinien sumować się do kwartalnych.

Dalej, robię to samo dla roku 2020. Analogicznie widzę kwoty:  2 825 oraz  3 383 . Różnica: -558

Sprawdźmy dalej po kolei: 

2019:  4 298 vs 4 487 . Różnica =  -189

2018:  5 604 vs  5 511 . Różnica = 93

2017:  7 173  vs 7 173 . Różnica = 0

2016:   5 740 vs 5 740 . Różnica = 0

2015:   3 233 vs  3 233 . Różnica = 0

2014: (5 828) vs (5 828) . Różnica = 0


Aktualizacja:
Sprawdziłem sprawozdanie roczne. Wygląda na to, że mamy tu do czynienia z odwróceniem odpisów aktualizujących. No tak, miliard w te czy wewte, co za różnica? To przykład, jak bardzo należy być krytycznym w stosunku do zysku netto. Bez tego odwrócenia zysk na akcję wyniesie niecałe 24 zł. Warto też zauważyć, że Orlen ma niedużo środków pieniężnych na akcję, bo niecałe 7 zł. To też jest przyczyna, dlaczego dywidendy są tak niewielkie.

poniedziałek, 1 sierpnia 2022

Czy opłaca się kupować PKN Orlen?

 Ostatnio głośno o rekordowych wynikach PKN Orlen. W 2021 zwiększył zysk netto dla akcjonariuszy o ponad 300%, a w 1 kw 2022 o 50% r/r. Mimo to kurs wcale nie urósł, jak można by się tego spodziewać. Ale czy na pewno? Zauważmy, że zysk na akcję wyniósł na 2021 26 zł, a wypłacona dywidenda w 2022 to zaledwie 3,5 zł. Większość zysku przeznaczono na inwestycje. 

Problemem przy wycenie PKN jest duża zmienność jego wyników. Łatwiej podejść do ich prognozy od strony rentowności, która jest bardziej stabilna - ok. 11-11,5%. Jest to poziom długoterminowego kosztu kapitału własnego dla WIGu i powiedzmy, że mniej więcej też na PKN. Spółka przeznacza ostatnio coraz mniej na dywidendę, poniżej 25%. Powiedzmy, że za rok będzie to 20%. W takim razie oczekiwana stopa wzrostu wynosi 80%*11% = 8,8%. W takim razie dywidenda wzrosłaby za rok do ok. 3,8 zł. W optymistycznym wariancie, gdyby spółka potrafiła utrzymać taki wzrost, wartość akcji wyniosłaby 3,8/(0,11-0,088) = 173 zł. Jak widać jednak inwestorzy nie są tak optymistyczni. Długoterminowy nominalny wzrost gosp. Polski to ok. 6% i gdyby go użyć wtedy 3,8/(0,11-0,06) = 76 zł. Czyli tyle, ile obecnie wynosi kurs. Jednak poprawniejsza wycena powinna uwzględnić wyższy krótkoterminowy koszt kapitału. W dużej mierze wyższa cena ropy poprawiła wyniki PKN (wyższe przychody minus wyższe koszty o tyle samo daje ciągle wyższe zyski dzięki różnicy procentowej), ta przyczyniła się do wyższej inflacji (chociaż nie przesadzajmy z jej efektem, bo kosztuje tyle co 10 lat temu), rządy przestraszyły się inflacji (opóźniona reakcja o 2 lata - normalka) i zaczęły podnosić stopy. Wyższa stopa procentowa to wyższa rentowność obligacji, a więc i akcji, gdyż inwestorzy wymagają co najmniej takiego samego wzrostu na ryzykownych aktywach. To zatrzymuje wzrost kursów. Dodatkowo pojawiają się sygnały o nadchodzącej recesji. W sumie krótkoterminowy koszt kapitału wynosi obecnie co najmniej 15%, co można ustawić na kilka lat wprzód. Wtedy wycena wyjdzie na poziomie ok. 75-75,5 zł, tyle co dziś w środku sesji.

Możemy spróbować wyznaczyć 3-wartościowy quasi-kanał, w jakim kurs powinien obecnie oscylować. Wykorzystam do tego swoje wzory na wycenę inwestora krótkoterminowego i długoterminowego. 

Dla krótkoterminowego (zob. ten artykuł):

Dla długoterminowego (zob. ten artykuł):


Dla dywidendowego (zwykły wzór Gordona):

gdzie:

P - wartość wewnętrzna akcji

D - oczekiwana dywidenda brutto w kolejnym okresie

x - koszt transakcji od wartości akcji

k - koszt kapitału własnego (bez odejmowania podatków)

g - oczekiwana stopa wzrostu dywidendy brutto


Krótkoterminowy inwestor trzyma akcje rok i sprzedaje innemu krótkoterminowemu inwestorowi. Długoterminowy inwestor trzyma akcje wiele lat i sprzedaje innemu długoterminowemu. Inwestor dywidendowy nie sprzedaje nigdy akcji - stąd nie wliczamy kosztów transakcyjnych (od zakupu pomijamy).

Wycena dla krótkoterminowego:

P = 3.8 / (0.11 - (0.06 - 2*0.0039*(1+0.06)/(1+0.0039))) = 65 zł

Wycena dla długoterminowego:

P = 3.8/(0.11-(0.06*(1-0.0039)-0.0039)) = 70

Wycena dla dywidendowego:

P = 3.8/(0.11-0.06) = 76

Dla pierwszych lat koszt kapitału wyniesie 15-17%, ale jest to kompensowane wyższą oczekiwaną stopą wzrostu, stąd można od razu posługiwać się długoterminowymi parametrami. Dlatego kanał między 65 a 76 jest jak najbardziej racjonalny.

Oczywiście to powierzchowna analiza i wielu rzeczy nie uwzględnia, ale na ten moment ja nie widzę sensu kupowania akcji. 


P. S. Co do fuzji PKN Orlen i PGNiG ( (źródło)

- analityk Trigon uważa, że "dziś trudno oczekiwać istotnych efektów synergii wynikających z przeprowadzonego połączenia. Już przy fuzji Orlenu z Lotosem były one mgliste, a w przypadku fuzji z PGNiG są tym bardziej wątpliwe – twierdzi Kozak. Zauważa, że Orlen i PGNiG działają w nieco innych obszarach rynku, a to mocno utrudnia uzyskanie prostych synergii."

- analityk mbanku potwierdza: "połączenie Orlenu z PGNiG prawdopodobnie nie przyniesie istotnych synergii. Trudno będzie o nie w sytuacji, gdy obie grupy koncentrują swoją działalność w innych obszarach rynku".

Analitycy widzą też plusy w postaci większych aktywów, w tym większych możliwości regulowania zobowiązań. To jednak nie ma dużego znaczenia, a jedynie synergia przeniesie się na wartość akcji. Warto przypomnieć, że nieuzasadniona dywersyfikacja obniża cenę akcji (zob. art. Czy firma powinna dywersyfikować produkt? ). To przecież miało miejsce w przypadku kupna Polska Press przez Orlen. Wydano pieniądze, a żadna z tego poprawna inwestycja. Była ta jedynie polityczna inwestycja, niekorzystna dla inwestorów (zob. np. art. Wyborczej Orlen kupił wydmuszkę. Grupa Polska Press od dawna pozbywała się majątku). Biorąc pod uwagę te polityczne wydatki, należy się dwa razy zastanowić czy prezes Obajtek bardziej nie dba o interesy swojej partii niż firmy.

środa, 27 lipca 2022

Jak sprawdzić w języku R czy wariancja się zmieniła?

Porównując zyski z akcji, funduszy lub firm powinniśmy nie tylko sprawdzić stabilność średniej, ale też wariancji. Jeżeli mamy dwie spółki, które średnio biorąc zarabiają tyle samo, ale jedna może mieć zarówno stratę jak i duży zysk, a druga generuje nieduży, ale zawsze stabilny zysk, to jasne, że w kategoriach inwestycji ocenimy wyżej tę drugą. Ale to tylko pierwsza część analizy. Bo gdyby to było jedyne kryterium oceny (tego co by można nazwać jakością zysku), to wystarczyłoby zastosować współczynnik zmienności czy Sharpe'a. A przecież może się zdarzyć, że całkowita wariancja będzie równa w obu firmach, ale jedna spółka będzie miała wyższą zmienność w przeszłości, a obecnie niższą, natomiast druga historycznie ma taką samą. Możliwe, że wtedy inwestor wolałby pierwszą, która ma niższe bieżące odchylenia. Chociaż to zależy od innych czynników, to ktoś inny mógłby wybrać drugą, którą uznałby za bardziej przewidywalną. 

Przetestować wariancję można na różne sposoby, np. testem Goldfelda-Quandta (GQ), Breuscha-Pagana (BP) czy Harrisona-McCabe (HM). Te i jeszcze parę innych znajdują się w pakiecie lmtest, z którego skorzystamy. Oczywiście można także spróbować zbadać efekt ARCH, ale dotyczy on warunkowej wariancji, która na ten moment nas nie interesuje.

Kontynuujemy 3 poprzednie artykuły dotyczące programowania w języku R. Jeżeli pewne kwestie są niezrozumiałe, to prawdopodobnie wyjaśnione są tam właśnie. Dobry test na stałość wariancji powinien działać zarówno na danych skorelowanych, jak i nieskorelowanych. Przetestujmy moc BP, GQ i HM. Aby ich użyć, zainstalujemy pakiet lmtest. Więc skrypt zaczynamy od:

if (require(lmtest)==FALSE) {

  install.packages("lmtest")

  library("lmtest")

}

Zanim przejdę do testowania mocy obydwu testów, pokażę prosty przykład ich zastosowania.

Symulujemy proces ARMA(1, 1). Odchylenie standardowe (sd) wzrasta w połowie próby z 1 na 2. Przykładowo:

set.seed(199)
sim1 = rnorm(n=dlugosc_proby, sd=1)
sim2 = rnorm(n=dlugosc_proby, sd=2)
sim12 = c(sim1, sim2)
plot(sim12, type="l")



Test BP, GQ i HM wykonamy przy użyciu kodu:

czas = c(1:length(sim12))
bptest(sim12 ~ czas)
gqtest(sim12 ~ czas)
hmctest(sim12 ~ czas)

Efekt:


studentized Breusch-Pagan test

data:  sim12 ~ czas
BP = 13.843, df = 1, p-value = 0.0001987

> gqtest(sim12 ~ czas)

Goldfeld-Quandt test

data:  sim12 ~ czas
GQ = 5.323, df1 = 48, df2 = 48, p-value = 2.222e-08
alternative hypothesis: variance increases from segment 1 to 2

Harrison-McCabe test

data:  sim12 ~ czas
HMC = 0.15859, p-value < 2.2e-16

Trzy testy wskazały prawidłowo p < 0,05.

Zaczynamy badanie mocy dla różnych przypadków.

Część stała kodu będzie następująca - nie będę już jej powtarzał potem, ale wykonujemy to za każdym razem : 

set.seed(NULL)

liczba_prob = 1000

dlugosc_proby=50

sukces1=0

sukces2=0

sukces3=0

Część 1) Stała wariancja

A) bez zmiany średniej i autokorelacji

 Zaczniemy od standardowej sytuacji błądzenia losowego bez zmian wariancji. Tworzymy próbkę losową z rozkładu normalnego (50 danych). Najprostszy możliwy kod byłby następujący:

# brak zmian wariancji, brak autokorelacji

czas = c(1:(dlugosc_proby))

for (i in 1:liczba_prob) {

  sim1 = rnorm(n=dlugosc_proby)

  p1 = as.numeric(bptest(sim1 ~ czas)[4])

  if (p1>0.05) {

    sukces1 = sukces1 + 1

  }

  p2 = as.numeric(gqtest(sim1 ~ czas)[5])

  if (p2>0.05) {

    sukces2 = sukces2 + 1

  }

  p3 = as.numeric(hmctest(sim1 ~ czas)[3])

  if (p3>0.05) {

    sukces3 = sukces3 + 1

  }

moc_bp = sukces1/liczba_prob

moc_bp

moc_gq = sukces2/liczba_prob

moc_gq

moc_hm = sukces3/liczba_prob

moc_hm


BP dał 96%, GQ 95% i HM 95%. mocy.

B) z autokorelacją

Sprawdźmy to samo, ale w przypadku AR(1), np. phi 0.5. Zmieniamy tylko 1 linię

# brak zmian wariancji, stały AR(1)

  sim1 = arima.sim(list(order=c(1,0,1), ar = 0.5, ma=0.00001), n=dlugosc_proby)


Tym razem wynik pogorszył się: wszystkie testy dały po 91-92%, Co by się stało, gdyby phi wyniosło 0,9? Okazuje się, że BP pogorszył się do 71%, GQ do 80% i HM 81%. 

Co by się stało, gdyby w połowie próby AR(1) phi wyniosło 0.5, a od połowy próby wróciła do 0? Musimy już stworzyć dwie próby, więc kod się zmieni:

# brak zmian wariancji, zmienny AR(1)

  czas = c(1:(dlugosc_proby*2))

  for (i in 1:liczba_prob) {

   sim1 = arima.sim(list(order=c(1,0,1), ar = 0.5, ma=0.00001), n=dlugosc_proby)

    sim2 = arima.sim(list(order=c(1,0,1), ar=0.00001, ma=0.00001), n=dlugosc_proby)

    sim12 = c(sim1,sim2)

    p1 = as.numeric(bptest(sim12 ~ czas)[4])

    if (p1>0.05) {

      sukces1 = sukces1 + 1

    }

    p2 = as.numeric(gqtest(sim12 ~ czas)[5])

    if (p2>0.05) {

      sukces2 = sukces2 + 1

    }

    p3 = as.numeric(hmctest(sim12 ~ czas)[3])

        if (p3>0.05) {

      sukces3 = sukces3 + 1

    }

  } 

  moc_bp = sukces1/liczba_prob

  moc_bp

  moc_gq = sukces2/liczba_prob

  moc_gq

  moc_hm = sukces3/liczba_prob

  moc_hm


BP dał 89%, GQ 98% i HM 98%. Podobne efekty zobaczymy w przypadku MA(1). Jednak ARMA(1, 1) przy zmianie phi z 0,5 do 0 i theta 0,3 do 0:

# brak zmian wariancji: zmienny ARMA(1,1)
  sim1 = arima.sim(list(order=c(1,0,1), ar = 0.5, ma=0.3), n=dlugosc_proby)
  
  sim2 = arima.sim(list(order=c(1,0,1), ar=0.00001, ma=0.00001), n=dlugosc_proby)


znowu pogarsza wyniki BP (67%), ale poprawia GQ (99%) i HM (99%).

-----------------------------------------------------------------------------------------
Uwaga. Powyższy kod jest prawidłowy tylko dla sytuacji, gdy mamy spadek wartości parametrów do zera. Gdybyśmy chcieli zmienić parametry ARMA na inne niż zerowe, to już nie można tak prosto zapisywać obu symulacji. W poprzednim wpisie pokazywałem, jak należy to robić dla samego AR(1). A teraz pokażę, jak należy zapisać dla ARMA. Czyli np. chcemy np. zmienić parametry odpowiednio z 0,5 i 0,3 na na 0.4 i 0.2. Nie możemy po prostu ich zastąpić i resztę zostawiać bez zmian, bo proces zaczyna się nie od wartości losowej, tylko od zakończenia poprzedniej symulacji - bo ciągle występuje autokorelacja. W przypadku gdy współczynniki spadają do zera, to autokorelacja też znika, więc nie trzeba się tym martwić. Dlatego badając zmiany najłatwiej zacząć od niezerowych, a skończyć na zerowych parametrach, a nie odwrót. Jeżeli jest to niemożliwe, to trzeba wykonać trudniejszy kod. 

Przede wszystkim trzeba sobie przyswoić, że symulacje są wykonywane często z tzw. wypalaniem (ang. burn-in). Chodzi tu o pominięcie pierwszych wartości po to, aby startowe wartości nie wpływały na przebieg symulowanego procesu. Powiedzmy, że standardowy rozkład normalny i przypadkowo pierwsze wartości mogą być bliskie 3 odchyleniom standardowym, czyli wynieść ok. 3. Teraz druga wartość ARMA powstaje poprzez przemnożenie 3 przez phi, np. 0,5, co daje 1.5 i dodatkowo theta, np. 0.3, jest znów mnożona przez 3, co daje 0,9. Do tego dodajemy jakieś znów odchylenie, np. 1.  Czyli druga obserwacja wynosi razem 3.4. Powiedzmy, że kolejna losowa wartość wynosi 2. Znowu sporo, tak że trzecia obserwacja wyniesie 0,5*3,4 + 0,3*1 + 2 = 4. Czyli zauważmy, że powstaje ciąg 3, 3.4, 4, którego wartości rosną. W końcu zaczną spadać do zera, bo taka jest średnia, ale zanim to nastąpi może minąć trochę czasu, innymi słowy dane muszą się "wygrzać". Proces wygrzewania jest czysto przypadkowy, więc ucinanie danych nie jest konieczne, ale jest pomocne w przypadku symulowania niewielkiej liczby danych; wtedy bowiem mamy pewność, że nawet mała próba jest reprezentatywna. Poza tym wypalanie czy wygrzewanie będzie miało istotne znaczenie przy rozkładach z grubymi ogonami - jeśli akurat trafi się dalekie odchylenie, to powrót do średniej będzie trwał za długo. Proces będzie wydawał się niestacjonarny. Przy okazji warto zauważyć tutaj związek między grubymi ogonami a długą pamięcią - ten powrót do średniej będzie powodował autokorelacje dalekiego zasięgu. Te są uwzględnione dopiero w procesie ARFIMA.

W każdym razie jeśli ręcznie nie ustawimy liczby opuszczonych wartości losowych, to algorytm zrobi to sam. W przypadku AR(1) należy pominąć co najmniej 1 wartość - zapisujemy to przez n.start = 1. W przypadku ARMA(1, 1) będą to 2 wartości. 

Jak to się wiąże z poprzednim akapitem? Jeżeli chcemy mieć powiązanie drugiej próby z pierwszą, to musimy algorytmowi ustawić wartości startowe za pomocą start.innov. Wynika z tego, że w drugiej próbie wartością startową jest ostatnia wartość pierwszej próby. Problem polega na tym, że algorytm wypalania danych wymaga - jak już wyżej napisałem - pominięcia co najmniej dwóch wartości w przypadku ARMA(1, 1). Stąd musimy ustawić nie jedną, a dwie pierwsze wartości startowe, tzn. dwie ostatnie z pierwszej próby. Czyli aby precyzyjnie powiązać obie próby, drugą powinniśmy zapisać następująco:

sim2 = arima.sim(list(order=c(1,0,1), ar=0.4, ma=0.2), n=dlugosc_proby, start.innov=c(sim1[dlugosc_proby-1], sim1[dlugosc_proby]), n.start=2)

Co więcej, co pierwszej także powinniśmy dopisać opcję n.start=2 , aby wypalanie danych z drugiej próby synchronizowało z wypalaniem w pierwszej próbie. Czyli sim1 będzie zapisana tak:

sim1 = arima.sim(list(order=c(1,0,1), ar=0.5, ma=0.3), n=dlugosc_proby, n.start=2) 

Dla eksperymentu porównamy ARMA z błądzeniem przypadkowym z tego samego ziarna. Uruchomijmy poniższy kod:

dlugosc_proby=20

set.seed(7003)

sim1 = arima.sim(list(order=c(1,0,1), ar=0.1, ma=0.1), n=dlugosc_proby, n.start=2)

sim2 = arima.sim(list(order=c(1,0,1), ar=0.1, ma=0.1), n=dlugosc_proby, start.innov=c(sim1[dlugosc_proby-1], sim1[dlugosc_proby]), n.start=2)

sim12 = c(sim1,sim2)

plot(sim12,type="l")

set.seed(7003)

r1 = rnorm(dlugosc_proby+2)[3:(dlugosc_proby+2)]

r2 = rnorm(dlugosc_proby)

r12 = c(r1, r2)

lines(r12, lty=2)

legend("topright", legend=c("ARMA(1,1)", "błądzenie losowe"), lty=c(1,2))

Efekt:



Rysunek wskazuje, że ARMA(1,1) o niskiej autokorelacji niemal pokrywa się z błądzeniem losowym. Oczywiście, jeśli ustawimy phi i theta na zero, to nie będzie między nimi żadnych różnic. Następnie możemy już drugą próbę zmienić - podwyższyć phi i theta do 0,2:


A więc od połowy sumarycznej próby następują zwiększone różnice spowodowane występowaniem autokorelacji. A dodatkowo do drugiej próby dodamy wyraz wolny, np. 5:


Rysunek ten dowodzi, że przeszliśmy płynnie od jednej próby do drugiej. 

Zwrócę jeszcze uwagę na dwie linie:

r1 = rnorm(dlugosc_proby+2)[3:(dlugosc_proby+2)]

r2 = rnorm(dlugosc_proby)

pierwsza jest oczywista: chcę dostać zmienną z rozkładu normalnego od 3-go punktu, bo 2 pierwsze pomijam, aby dopasować do ARMA. Dlaczego więc w drugiej nie muszę już tego robić? Bo druga próba ARMA nie zaczyna się od wartości losowych - ustawiłem dwie pierwsze wartości startowe jak chciałem i te dwie wartości właśnie algorytm pominął... Pamiętajmy, że nie chodzi o pominięcie w obliczeniach, bo one są uwzględniane, a jedynie pomijamy je w ostatecznej symulacji. Gdy algorytm pominął już dwie stałe, to dalej do odchyleń bierze już zmienną gaussowską od początku danego ziarna. Kiedy się to czyta pierwszy raz, to wydaje się to pewnie mocno skomplikowane, ale tak naprawdę, jest to zupełnie logiczne. 

-----------------------------------------------------------------------------------------


C) ze zmianą średniej
To zmieńmy samą średnią (wyraz wolny), np. o 3:

# brak zmian wariancji, zmiana średniej

  sim1 = arima.sim(list(order=c(0,0,0)), n=dlugosc_proby, mean = 0)

  sim2 = arima.sim(list(order=c(0,0,0)), n=dlugosc_proby, mean = 3)

BP -> 98%  , GQ -> 96%, HM -> 97,5%.



Część 2) Zmiana wariancji

Zmienimy odchylenie standardowe (sd) z 1 na 2.
A) bez zmiany średniej i autokorelacji

# zmiana wariancji, brak autokorelacji

czas = c(1:(dlugosc_proby*2))

for (i in 1:liczba_prob) {
  sim1 = arima.sim(list(order=c(0,0,0)), n=dlugosc_proby, sd=1)
  sim2 = arima.sim(list(order=c(0,0,0)), n=dlugosc_proby, sd=2)
  
  sim12 = c(sim1, sim2)
  
  p1 = as.numeric(bptest(sim12 ~ czas)[4])
  
  if (p1<0.05) {
    sukces1 = sukces1 + 1
  }
  
  p2 = as.numeric(gqtest(sim12 ~ czas)[5])
  
  if (p2<0.05) {
    sukces2 = sukces2 + 1
  }

    p3 = as.numeric(hmctest(sim12 ~ czas)[3])

        if (p3>0.05) {

      sukces3 = sukces3 + 1

    }

  } 

  moc_bp = sukces1/liczba_prob

  moc_bp

  moc_gq = sukces2/liczba_prob

  moc_gq

  moc_hm = sukces3/liczba_prob

  moc_hm



Odwrotny sprawdzian okazuje się być znów na korzyść GQ i HM: BP ->93,5% , GQ -> 99,9%, HM -> 99,9%. 

B) z autokorelacją

 Dodajmy autokorelacje, AR(1), phi = 0.5

# zmiana wariancji + stały AR(1):
sim1 = arima.sim(list(order=c(1,0,0), ar=0.5), n=dlugosc_proby, sd=1)
sim2 = arima.sim(list(order=c(1,0,0), ar=0.5), n=dlugosc_proby, start.innov=sim1[dlugosc_proby], n.start=1, sd=2)

BP -> 91%, GQ -> 99%, HM -> 99,5%. 


Ze zmianą phi od połowy próby

# zmiana wariancji + zmienny AR(1):
  sim1 = arima.sim(list(order=c(1,0,0), ar=0.5), n=dlugosc_proby, sd=1)
  sim2 = arima.sim(list(order=c(1,0,0), ar=0.00001), n=dlugosc_proby, sd=2)

BP -> 83%, GQ -> 98%, HM -> 97%

Zmiana ARMA(1, 1): spadek phi i theta do 0 od połowy próby:

# zmiana wariancji + zmienny ARMA(1,1):
  sim1 = arima.sim(list(order=c(1,0,1), ar=0.5, ma=0.3), n=dlugosc_proby, sd=1)
  sim2 = arima.sim(list(order=c(1,0,1), ar=0.00001, ma=0.00001), n=dlugosc_proby, sd=2)

BP -> 63%, GQ -> 89%, HM -> 87%.


C) z autokorelacją i zmianą średniej

To na koniec zróbmy ten najbardziej złożony przykład: wzrost phi z 0.5 do 0,6 i theta z 0.1 do 0,2 oraz zmianą średniej z 1 do 3:

# zmiana wariancji + zmiana ARMA(1, 1) + zmiana średniej:

  sim1 = arima.sim(list(order=c(1,0,1), ar=0.5, ma=0.1), n=dlugosc_proby, n.start=2, mean=1, sd=1)

  sim2 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.2), n=dlugosc_proby, start.innov=c(sim1[dlugosc_proby-1], sim1[dlugosc_proby]), n.start=2, mean=3, sd=2)   


BP -> 68%, GQ -> 99,8%, HM -> 83%.

Najsłabszym testem okazał się BP, najlepszy GQ.


Część 3) Testowanie rzeczywistych danych

A) Stopa zmian PKB USA w latach 1798-2020. Źródło: https://www.measuringworth.com/




To są te same dane, które analizowałem w poprzednim artykule. Już wtedy test niestabilności modelu MA(1) wykrył zmianę w roku 1929. Zobaczmy czy mogłoby być to spowodowane zmianą wariancji. Zmienną nazwałem stopa i wtedy robimy:

czas = c(1:length(stopa))
bptest(stopa ~ czas)
gqtest(stopa ~ czas)
hmctest(stopa~ czas)
 

Efekt:

studentized Breusch-Pagan test

data:  stopa ~ czas
BP = 0.64387, df = 1, p-value = 0.4223

Goldfeld-Quandt test

data:  stopa ~ czas
GQ = 1.688, df1 = 110, df2 = 109, p-value = 0.003301
alternative hypothesis: variance increases from segment 1 to 2

Harrison-McCabe test

data:  stopa ~ czas
HMC = 0.36963, p-value = 0.005


Rzeczywiście GQ i HM wykryły zmianę w wariancji, zaś BP nie poradził sobie.


B) WIG - miesięczne stopy zwrotu (01.1995-06.2022). Źródło: stooq.pl


 
Wykres pozwala intuicyjnie ocenić, że wariancja zmieniła się mniej więcej w 2010 r. Zobaczmy co mówią testy. Jest to 330 obserwacji, więc wystarczająco dużo, aby obiektywnie ocenić.

Wykonujemy kod:

czas = c(1:length(stopa))
bptest(stopa ~ czas)
gqtest(stopa ~ czas)
hmctest(stopa~ czas)

studentized Breusch-Pagan test

data:  stopa ~ czas
BP = 20.154, df = 1, p-value = 7.145e-06

> gqtest(stopa ~ czas)

Goldfeld-Quandt test

data:  stopa ~ czas
GQ = 0.38905, df1 = 163, df2 = 163, p-value = 1
alternative hypothesis: variance increases from segment 1 to 2

Harrison-McCabe test

data:  stopa ~ czas
HMC = 0.71949, p-value = 1


Spore zaskoczenie: BP wskazuje na silną heteroskedastyczność, podczas gdy GQ i HM odrzucają ją. Pierwszy pokazuje zerowe p-value, pozostałe dokładnie przeciwne. Co to może oznaczać? Jest kilka możliwości. Przypomnijmy, że BP dość słabo radził sobie z autokorelacjami. Jeżeli to co widzimy na wykresie jest nimi spowodowane, a nie zmiennością, to GQ i HM mają rację. W dodatku autokorelacje nie muszą ograniczać się tylko do stóp zwrotu, ale dotyczą też ich warunkowych wariancji. Jeżeli występuje efekt ARCH, to niewarunkowa wariancja pozostaje stała, a zmienia się jedynie warunkowa (zob. ten art.). Zakładam więc, że BP będzie wrażliwy na takie zmiany, a GQ i HM nie, ale to trzeba sprawdzić. Formalnie GQ i HM mogą mieć rację, ale BP poinformuje, że autokorelacje wpływają na nasz proces. Co więcej, autokorelacje o których tu mówimy dotyczą krótkiej pamięci procesu, podczas gdy może występować także długa pamięć w postaci ARFIMA (dla stóp zwrotu) lub FIGARCH (dla warunkowych wariancji). W końcu powinno się także zbadać, czy długie ogony nie mają przypadkiem wpływu na wyniki. Nasze symulacje ograniczały się do rozkładu normalnego, a jak wiadomo stopy zwrotu znacznie się od niego różnią. 

niedziela, 3 lipca 2022

Trzy testy F na stabilność modelu ARMA i porównanie ich z KPSS. Przykłady ekonomiczne i giełdowe

Testowanie stabilności modelu ekonometrycznego przeprowadza się w dwóch celach. Po pierwsze poszukuje się momentu zmiany wartości parametrów modelu, a po drugie odpowiedzi czy występuje zmiana (auto)korelacji w czasie, ewentualnie wariancji i kowariancji. Te dwa cele pozwalają odróżnić testy stabilności od testów stacjonarności. Przeprowadzając te drugie nie dowiemy się, w którym momencie mogła zajść zmiana np. trendu. Ponadto test stacjonarności bada przebieg jednej zmiennej, jednowymiarowego szeregu, podczas gdy test stabilności koncentruje się na relacji między zmienną zależną a niezależną, co właśnie definiuje model.

Powyższe prowadzi do wniosku, że gdy chcemy jedynie sprawdzić czy średnia stopa zwrotu zmieniła się po czasie, to powinniśmy stosować test stacjonarności i nie potrzebujemy tak naprawdę testu stabilności. Należy nadmienić, że powstały testy, które oddzielają załamania strukturalne od stacjonarności (czy przed i po załamaniu proces jest stacjonarny), ale jest to kwestia konwencji i przyjętego celu.    

Powstaje ważne pytanie jak testować zmienną strukturę parametrów w modelu ARIMA(p, d, q). Rozbiję ten problem na 3 części. Na koniec zrobię 3 przykłady na prawdziwych danych. Przedstawiona analiza jest kontynuacją dwóch poprzednich artykułów (tutu).

Część 1)

Powiedzmy, że ignorujemy fakt, że w analizowanych szeregach występują autokorelacje. Chcemy sprawdzić czy i w którym momencie ich średnia zmienia się w czasie. W art. na temat różnicy między testem t-Studenta i testem Chowa pokazałem, że można użyć zarówno t-Studenta (Welcha), jak i Chowa, ale w praktyce lepszą metodą jest test QLR. Ponieważ w R będą potrzebne dwa zewnętrzne pakiety - strucchange i forecast, to na początku możemy wykonać kod:


if (require(strucchange)==FALSE) {

  install.packages("strucchange")

  library("strucchange")

}

if (require(forecast)==FALSE) {

  install.packages("forecast")

  library(forecast)

}

Dodatkowo, instalujemy pakiet tseries, który przyda się do testowania stacjonarności:

if (require(tseries)==FALSE) {

  install.packages("tseries")

  library("tseries")

}


Zacznijmy od tych samych przykładów co w art. Symulacja i automatyczne modelowanie ARIMA w języku R. Przykład AR(1) z phi = 0,7:



# brak zmian (tak jest źle)
set.seed(1)
sim1 = arima.sim(list(order=c(1,0,0), ar=0.7), n=100)
qlr = Fstats(sim1 ~ 1)
plot(qlr, alpha=0.1)
sctest(qlr, type="supF")
breakpoints(qlr)

Wyniki:



Pierwszy wynik (brak istotności) sugerowałby, że QLR nadaje się do badania AR. To weźmy kolejne ziarno. Ponieważ jedynie zmienia się ten element, to nie powtarzam już kodu, jedyną zmianą jest pierwsza linia: set.seed(2). Wyniki:



Weźmy jeszcze przykładowo ziarno nr 20

    


Zatem już dostaliśmy sygnał, że test błędnie wskazuje na zmianę struktury. Test źle działa, bo nie uwzględniamy autokorelacji, które stanowią część AR lub MA. W symulowanym przykładzie wiemy, że jest to proces AR(1). Gdyby były to rzeczywiste dane, musielibyśmy znaleźć najbardziej dopasowany model, co uzyskamy dzięki funkcji auto.arima (z pakietu forecast). Dokładniej wpisujemy wtedy auto.arima(x, ic = "bic"), gdzie x to nasze dane. Później zrobię przykład. W każdym razie widzimy, że za pomocą QLR (a więc też Chowa i Welcha) nie można sprawdzać samych stóp zwrotu z akcji, jeśli występują jakieś autokorelacje.

Część 2)

Uwzględniamy autokorelacje. Poprawa testu jest prosta, ale trzeba stworzyć zmienną z poprzedniego okresu i następnie wpisać do testu model, tutaj AR(1). Dla ziarna 2 będzie to kod:

# brak zmian (podejście prawidłowe)
set.seed(2)
sim0 = arima.sim(list(order=c(1,0,0), ar=0.7), n=100)
sim1_1 = sim0[1:(length(sim0)-1)]
sim1 = sim0[2:length(sim0)]

qlr = Fstats(sim1 ~ 1+sim1_1)
plot(qlr, alpha=0.05)
sctest(qlr, type="supF")
breakpoints(qlr)


Statystyka F znalazła się poniżej progu istotności, czyli rzeczywiście nie znaleziono dowodów na zmianę struktury.

Dla ziarna 20 wykres będzie podobny:


Część 3)

To teraz sprawdźmy odwrotnie - czy test uchwyci faktyczne zmiany. Zauważmy, że w przypadku modelu ARMA mogą się ogólnie zmienić trzy parametry - wyraz wolny, wyraz AR i wyraz MA. 

A) zmiana wyrazu wolnego (średniej)

Przykład 1. Załóżmy, że wyraz wolny zmieni się  w połowie o 1 (powiedzmy dla seed 999):

# zmiana średniej (stopy zwrotu itp.)
dlugosc_proby=200
set.seed(999)
sim0 = arima.sim(list(order=c(1,0,0), ar=0.7), n=dlugosc_proby)
sim1=sim0[1:100]
sim2 = sim0[101:200]+1
sim12 = c(sim1, sim2)
sim12_1 = sim12[1:(length(sim12)-1)]
sim121 = sim12[2:length(sim12)]

par(mfrow=c(2,1), mar=c(2,5,2,2))
plot(sim12,type="l")

qlr = Fstats(sim121 ~ 1+sim12_1)
plot(qlr, alpha=0.05)
sctest(qlr, type="supF")
breakpoints(qlr)
par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)



Pierwszy wykres stanowi symulowany proces AR(1), a drugi statystykę F dla każdego punktu. Test QLR prawidłowo wskazał zmianę mniej więcej w połowie zakresu.  

Przykład 2. set.seed(13000), reszta bez zmian. Efekt:


QLR nie wykrył zmiany, dostałem p-value = 0,094. Poszukiwanie maksimum (supremum) F nie jest jednak jedynym kierunkiem oceny istotności statystyki F. Jak podaje A, Zeileis et al. [1] wypracowano 3 różne testy F:



gdzie i można interpretować jako punkt w czasie, a F(i) to statystyka F w punkcie i.


Dotąd przyjmowałem tylko pierwszy z nich. Drugi to średnia między kresem górnym i dolnym, a trzeci - najbardziej skomplikowany - średnia eksponencjalna. Chociaż tu trochę zgaduję, to sądzę, że można ją utożsamiać z medianą; zauważmy bowiem, że uśredniane są tu połowy wartości, a każde przekształcenie mediany z oryginalnego rozkładu staje się medianą przekształconego rozkładu [2].

Dlaczego obliczenie średniej statystyki może być lepsze niż poszukiwanie max? Bo inaczej jest obliczane p-value. W supF sprawdzamy czy największe F znajduje się powyżej poziomu uznanego jeszcze za losowy (np. szansa 5%). W aveF i expF sprawdzamy czy pewne średnie są istotnie większe niż wynikające z teoretycznego rozkładu zmiennej czysto losowej.  

Wróćmy do przykładu. Teraz staje się jasne, dlaczego dotąd stosowałem zapis sctest(qlr, type="supF"). Zamiast supF możemy wpisać aveF lub expF. Nie będzie już prawidłowe używanie zapisu qlr - bo to już nie będzie test Quandta - dlatego, aby uogólnić zapis, zamienię nazwę qlr na fs. Dodatkowo type zmienię na aveF:

# zmiana średniej (stopy zwrotu itp.)
dlugosc_proby=200
set.seed(999)
sim0 = arima.sim(list(order=c(1,0,0), ar=0.7), n=dlugosc_proby)
sim1=sim0[1:100]
sim2 = sim0[101:200]+1
sim12 = c(sim1, sim2)
sim12_1 = sim12[1:(length(sim12)-1)]
sim121 = sim12[2:length(sim12)]

fs = Fstats(sim121 ~ 1+sim12_1)
sctest(fs, type="aveF")
breakpoints(fs)

Tym razem nie tworzymy wykresów (stąd brak plot i par), bo chodzi tylko o zbadanie istotności wg innego algorytmu niż supF. W powyższym przypadku uzyskałem p-value = 0,023, a więc odrzucamy hipotezę o niezmiennej strukturze. Dla przypomnienia supF dał 0,094. Czyli widzimy, że aveF lepiej spełnił zadanie niż supF. 

Przykład 3. Zachowajmy aveF, ale zmieńmy ziarno na set.seed(1555. Tym razem dostaniemy p-value = 0.05748, a więc nieco brakuje do istotności. To wróćmy do supF. Tym razem supF dał p-value = 0.01394, czyli lepiej poradził sobie. Ponadto możemy użyć jeszcze expF - dostaniemy p 0.029, czyli też dobrze. 

Generalnie zaproponowałbym stosowanie następującej metody: jeśli którykolwiek z omówionych 3 testów wskazuje zmianę struktury, to uznajemy, że taka zachodzi. Wtedy dopiero możemy posłużyć się wykresem fs lub funkcją breakpoints do znalezienia momentu zmiany. 

Moc testu zależy jak dużą zmianę wpiszemy i jak dużo danych mamy. Np. powyżej zakładałem, że średnia stopa zwrotu zmieni się zaledwie o 1 pkt proc., więc test nie będzie silny. Przy 200 danych dostaniemy zaledwie 40% prawidłowych wyników. Jednak dla samego supF dostałem tylko 33% mocy. Nie ma co się dziwić, nawet w realnych warunkach 1 pkt proc. jest mało znaczący. Jednak wszystko się zmieni, jeśli ustalimy 2 pkty zmiany. Wykonajmy więc test mocy tych trzech testów F:

# sprawdzam moc - zmiana średniej
set.seed(NULL)
sukces = 0
dlugosc_proby=200
liczba_prób = 1000
zmiana = 2

for (i in 1:liczba_prób) {
  sim0 = arima.sim(list(order=c(1,0,0), ar=0.7), n=dlugosc_proby)
  sim1=sim0[1:100]
  sim2 = sim0[101:200]+zmiana
  sim12 = c(sim1, sim2)
  sim12_1 = sim12[1:(length(sim12)-1)]
  sim121 = sim12[2:length(sim12)]

  fs = Fstats(sim121 ~ 1+sim12_1)
  
  sc1 = sctest(fs, type="supF")
  p1 = as.numeric(sc1[2])
  
  sc2 = sctest(fs, type="aveF")
  p2 = as.numeric(sc2[2])

  sc3 = sctest(fs, type="expF")
  p3 = as.numeric(sc3[2])
  
  
  if (p1<0.05 || p2<0.05 || p3<0.05) {
    sukces = sukces + 1
  }

moc = sukces/liczba_prób
moc

set.seed(NULL) usuwa z pamięci ustawienia ziarna.

Taki test okazuje się już niezwykle silny - dostaniemy ponad 90% prawidłowych wyników. 

Jednak trzeba znowu przypomnieć, że test stacjonarności, jak KPSS, będzie dużo lepiej się sprawdzał w tym przypadku. Wykonamy ten kod:

# dla porównania test kpss (zmiana średniej)
set.seed(NULL)
sukces = 0
dlugosc_proby=200
liczba_prob = 1000
zmiana = 2
for (i in 1:liczba_prob) {
  
  sim0 = arima.sim(list(order=c(1,0,0), ar=0.7), n=dlugosc_proby)
  sim1=sim0[1:100]
  sim2 = sim0[101:200]+zmiana
  sim12 = c(sim1, sim2)
  
  p = kpss.test(sim12)[3]
  
  if (p<0.05) {
    sukces=sukces+1
  }
  

moc = sukces/liczba_prob
moc

Otóż dostaniemy moc bliską 100%.  Jest więc niemal pewne, że test wykryje każdą większą zmianę w stopie zwrotu. Nawet dla poprzedniej niewielkiej zmiany o 1, dostaniemy moc ok. 70% (testy F dawały 40%).

Można zapytać, po co robić test stabilności, skoro testy na stacjonarność dają lepsze rezultaty? Są dwie odpowiedzi. Po pierwsze taki test pokazuje jak na dłoni, w którym momencie występuje zmiana średniej, czego testy stacjonarności nie pokazywały. Po drugie testy na stacjonarność badają jedynie zmianę średniej i najwyżej wariancji. Test stabilności może nam pokazać czy i w którym momencie następują zmiany w korelacjach i autokorelacjach.

B) zmiana części AR

Ograniczymy się do 1 rzędu. 

Przykład 4. Powiedzmy, że do połowy próby występuje autokorelacja, a od drugiej połowy znika (phi spada z 0,5 do 0):

# zmiana nachylenia AR(1)
dlugosc_proby=200
set.seed(7003)
sim1 = arima.sim(list(order=c(1,0,0), ar=0.5), n=dlugosc_proby)
set.seed(7003)
sim2 = arima.sim(list(order=c(1,0,0), ar=0), n=dlugosc_proby)
sim12 = c(sim1, sim2)
sim12_1 = sim12[1:(length(sim12)-1)]
sim121 = sim12[2:length(sim12)]

par(mfrow=c(2,1), mar=c(2,5,2,2))
plot(sim12,type="l")

qlr = Fstats(sim121 ~ 1+sim12_1)
qlr
plot(qlr, alpha=0.05)

breakpoints(qlr)
par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)





Sprawdźmy moc tego testu:

# sprawdzam moc - zmiana nachylenia AR(1)
set.seed(NULL)
sukces = 0
dlugosc_proby=100
liczba_prob = 1000
zmiana = 0

for (i in 1:liczba_prob) {
  sim1 = arima.sim(list(order=c(1,0,0), ar=0.5), n=dlugosc_proby)
  sim2 = arima.sim(list(order=c(1,0,0), ar=0), n=dlugosc_proby)
  sim12 = c(sim1, sim2)
  sim12_1 = sim12[1:(length(sim12)-1)]
  sim121 = sim12[2:length(sim12)]

  fs = Fstats(sim121 ~ 1+sim12_1)
  
  sc1 = sctest(fs, type="supF")
  p1 = as.numeric(sc1[2])
  
  sc2 = sctest(fs, type="aveF")
  p2 = as.numeric(sc2[2])

  sc3 = sctest(fs, type="expF")
  p3 = as.numeric(sc3[2])
  
  
  if (p1<0.05 || p2<0.05 || p3<0.05) {
    sukces = sukces + 1
  }

moc = sukces/liczba_prób
moc


Taki test ma wysoką moc - dostałem ponad 85% prawidłowych wyników.


--------------------------------------------------------------------------------------
Uwaga: gdybyśmy chcieli zmienić wartość parametru na inną niż zero, to w powyższym kodzie nie można po prostu zmienić phi z zera na inną wartość. Np. gdybyśmy chcieli zmienić część AR z 0,5 na -0.5 , to poniższy kod na sim2:

sim2 = arima.sim(list(order=c(1,0,0), ar=-0.5), n=dlugosc_proby)

będzie błędny, dlatego że wartość początkowa tego szeregu będzie zmienną losową, więc kolejna modelu AR(1) powstanie przez pomnożenie tej wartości losowej przez -0,5, a powinna być przemnożona przez ostatni wyraz sim1. Prawidłowy kod na sim2 w tej sytuacji będzie następujący:

sim2 = arima.sim(list(order=c(1,0,0), ar=-0.5), n=dlugosc_proby, start.innov=sim1[dlugosc_proby], n.start=1)

W naszym przypadku te dodatkowe opcje są niepotrzebne, bo i tak ustalamy, że pierwsza wartość nie koreluje z poprzednim okresem.

--------------------------------------------------------------------------------------


C) zmiana części MA

Tak samo ograniczymy się do 1 rzędu. Testowanie MA jest bardziej skomplikowane, bo wymaga wyłuskania reszt modelu (tzn. odchyleń od niego), dlatego dla każdej podpróby należy najpierw wygenerować dany proces, a następnie utworzyć z niego model - tutaj MA. Dopiero mając faktyczny model będziemy mogli znaleźć od niego odchylenia, czyli reszty (tzn. symulowany proces minus model). Całą procedurę wykonujemy dwukrotnie, dla obu podprób.  

Przykład 5.
# zmiana nachylenia MA(1)
dlugosc_proby=100
set.seed(1099)
sim1 = arima.sim(list(order=c(0,0,1), ma=0.5), n=dlugosc_proby)
arma_stopa1 = arima(sim1, order=c(0,0,1))
reszty1 = arma_stopa1$residuals

set.seed(500)
sim2 = arima.sim(list(order=c(0,0,1), ma=0), n=dlugosc_proby)
arma_stopa2 = arima(sim2, order=c(0,0,1))
reszty2 = arma_stopa2$residuals

reszty12 = c(reszty1, reszty2)
reszty12_1 = reszty12[1:(length(reszty12)-1)]
reszty121 = reszty12[2:length(reszty12)]

sim12 = c(sim1, sim2)
sim12_1 = sim12[1:(length(sim12)-1)]
sim121 = sim12[2:length(sim12)]

par(mfrow=c(2,1), mar=c(2,5,2,2))
plot(sim12,type="l")
qlr = Fstats(sim121 ~ 1+reszty12_1)
plot(qlr, alpha=0.05)
sctest(qlr, type="supF")
breakpoints(qlr)
par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)





Zwróćmy uwagę, że naocznie byśmy nie byli w stanie zauważyć jakiejkolwiek zmiany w sim12.

Analogicznie do poprzednich przykładów, sprawdźmy moc za pomocą supF, aveF i expF:

# sprawdzam moc - zmiana nachylenia MA(1)
set.seed(NULL)
sukces = 0
dlugosc_proby=100
liczba_prob = 1000
zmiana = 0

for (i in 1:liczba_prob) {
  sim1 = arima.sim(list(order=c(0,0,1), ma=0.5), n=dlugosc_proby)
  arma_stopa1 = arima(sim1, order=c(0,0,1))
  reszty1 = arma_stopa1$residuals
  
  sim2 = arima.sim(list(order=c(0,0,1), ma=0), n=dlugosc_proby)
  arma_stopa2 = arima(sim2, order=c(0,0,1))
  reszty2 = arma_stopa2$residuals
  
  reszty12 = c(reszty1, reszty2)
  reszty12_1 = reszty12[1:(length(reszty12)-1)]
  reszty121 = reszty12[2:length(reszty12)]
  
  sim12 = c(sim1, sim2)
  sim121 = sim12[2:length(sim12)]
  

  fs = Fstats(sim121 ~ 1+reszty12_1)
  
  sc1 = sctest(fs, type="supF")
  p1 = as.numeric(sc1[2])
  
  sc2 = sctest(fs, type="aveF")
  p2 = as.numeric(sc2[2])
  
  sc3 = sctest(fs, type="expF")
  p3 = as.numeric(sc3[2])
  
  if (p1<0.05 || p2<0.05 || p3<0.05) {
    sukces = sukces + 1
  }

moc = sukces/liczba_prob
moc

Uzyskałem podobną moc co dla AR(1) - 85%. Oczywiście mówimy tu o próbie 200 obserwacji.


D) zmiana części AR i MA

Łączymy zmiany AR i MA. Jest tu najwięcej do napisania, ale jest to już mechaniczne, bo wystarczy powtórzyć to co było dotychczas.

Przykład 6.
# zmiana - phi i theta
dlugosc_proby=100
set.seed(61000)
sim1 = arima.sim(list(order=c(1,0,1), ar=0.5, ma=0.5), n=dlugosc_proby)
arma_stopa1 = arima(sim1, order=c(1,0,1))
reszty1 = arma_stopa1$residuals

set.seed(8899)
sim2 = arima.sim(list(order=c(1,0,1), ar=0.00001, ma=0.00001), n=dlugosc_proby)
arma_stopa2 = arima(sim2, order=c(1,0,1))
reszty2 = arma_stopa2$residuals

reszty12 = c(reszty1, reszty2)
reszty12_1 = reszty12[1:(length(reszty12)-1)]
reszty121 = reszty12[2:length(reszty12)]

sim12 = c(sim1, sim2)
sim12_1 = sim12[1:(length(sim12)-1)]
sim121 = sim12[2:length(sim12)]

par(mfrow=c(2,1), mar=c(2,5,2,2))

plot(sim12,type="l")
qlr = Fstats(sim121 ~ sim12_1 + reszty12_1)
plot(qlr, alpha=0.05)
sctest(qlr, type="supF")
breakpoints(qlr)
par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)




Aby nie pojawiały się ostrzeżenia, wpisałem zamiast zer w nowych parametrach liczby bliskie zer.

I test mocy:

# sprawdzam moc - zmiana phi i theta
set.seed(NULL)
sukces = 0
dlugosc_proby=100
liczba_prob = 1000
zmiana = 0
for (i in 1:liczba_prob) {
  sim1 = arima.sim(list(order=c(1,0,1), ar=0.5, ma=0.5), n=dlugosc_proby)
  arma_stopa1 = arima(sim1, order=c(1,0,1))
  reszty1 = arma_stopa1$residuals
  
  sim2 = arima.sim(list(order=c(0,0,0)), n=dlugosc_proby)
  arma_stopa2 = arima(sim2, order=c(0,0,0))
  reszty2 = arma_stopa2$residuals
  
  reszty12 = c(reszty1, reszty2)
  reszty12_1 = reszty12[1:(length(reszty12)-1)]
  reszty121 = reszty12[2:length(reszty12)]
  
  sim12 = c(sim1, sim2)
  sim12_1 = sim12[1:(length(sim12)-1)]
  sim121 = sim12[2:length(sim12)]
  

  fs = Fstats(sim121 ~ sim12_1 + reszty12_1)
  
  sc1 = sctest(fs, type="supF")
  p1 = as.numeric(sc1[2])
  
  sc2 = sctest(fs, type="aveF")
  p2 = as.numeric(sc2[2])
  
  sc3 = sctest(fs, type="expF")
  p3 = as.numeric(sc3[2])
  
  if (p1<0.05 || p2<0.05 || p3<0.05) {
    sukces = sukces + 1
  }

moc = sukces/liczba_prob
moc

Taki test ma już moc 100%. Natomiast test KPSS niemal w każdym przypadku będzie prawidłowo wskazywał na zachowanie stacjonarności próby, czyli "moc" wyniesie blisko zera.



Przeprowadźmy testy na prawdziwych danych.

Przykład 7. Stopa zmian PKB USA w latach 1798-2020. Źródło: https://www.measuringworth.com/

#wczytujemy dane, np. mam plik z latami w kolumnie 1 i stopami zmian w kolumnie 2, 
#plik znajduje się w katalogu roboczym (który ustawiamy w opcjach albo za pomocą funkcji setwd(), 
# powiedzmy Pobrane / Downloaded), więc od razu wczytujemy dzięki read.csv z pomocnym warunkiem:

plik = read.csv(nazwa_pliku, sep=";")
if (ncol(plik)==1) {
  plik = read.csv(nazwa_pliku, sep=",")
}

daty = plik[,1]
stopa = plik[,2]

# tworzymy szeregi czasowe (każda obserwacja ma okres)

stopa = ts(stopa, start=rok[1])

Po wczytaniu danych zaczynamy od testowania stacjonarności, np. za pomocą KPSS. Warto jednak przypomnieć, że wynik testu dość mocno zależy od wyboru okna spektralnego. W art. Jakie opóźnienie stosować w KPSS i ADF-GLS? wskazywałem, że autorzy sugerowali inne opóźnienia dla mniejszych i większych prób. W R te mniejsze opóźnienia są zmienną lshort, które ustawione są jako poziom domyślny (czyli jako TRUE). Jeżeli ustawimy je jako FALSE, dostaniemy opóźnienia dla większych prób. Chociaż nie ma wyraźnej granicy, kiedy stosować jedną, a kiedy drugą, to KPSS podali przykład 200 danych, dla których dobrze spełnia rolę lshort. Dlatego zrobiłem warunek:

# test stacjonarności
if (length(stopa)>200) {
  kpss.test(stopa, lshort=FALSE)
}  else {
    kpss.test(stopa)
  }


Efekt:


Przyznam, że nie rozumiem tych ostrzeżeń, które pojawiają się niemal za każdym razem. W dokumentacji do tej funkcji (pakietu forecast) zawarto krótkie zdanie, że jeśli obliczona statystyka znajduje się poza tabelą wartości krytycznych, pojawi się właśnie ostrzeżenie. Jaki to ma sens, skoro bez przerwy się pojawia? Ostrzeżenie powinno wskazywać na rodzaj wyjątku sugerując badaczowi, że może popełnił gdzieś błąd w danych. A to nie o to chodzi. Np. w Gretlu nie było takich komunikatów.  Uważam to za mały błąd tej funkcji, którą mam nadzieje, że autorzy kiedyś naprawią.

W każdym razie p-value powyżej 0.1 oznacza brak istotności, tzn. brak dowodów na niestacjonarność. Średnia stopa zwrotu się nie zmienia.

Kolejnym etapem będzie badanie stabilności. Aby to zrobić, najpierw poszukajmy najlepszego modelu arima:

#szukanie najlepszego modelu, wg BIC
arma_stopa = auto.arima(stopa,ic="bic")
arma_stopa

Efekt:


Najlepszym modelem jest MA(1). Aby go zastosować szukamy reszt i tworzymy model:

reszty = arma_stopa$residuals
reszty0 = reszty[-1]
reszty1 = reszty[-length(reszty)]

stopa0 = window(stopa, start=rok[2])
stopa1 = window(stopa, end=rok[length(rok)-1])

fs = Fstats(stopa0 ~ 1 + stopa1 + reszty1)


Aby  pokazać wykres zmian PKB i statystyki F powinniśmy dopasować osie x obu wykresów. Chodzi o to, że test F bierze do obliczeń 15% danych przed i po każdym okresie, więc jego zakres zawiera się w przedziale 0,15-0,85. A więc graficznie PKB też musimy tak zawęzić. Stąd tworzę nową zmienną o nazwie stopa_gr, którą definiuję tak:

stopa_gr = window(stopa, start=rok[ceiling(0.15*length(rok))], end=rok[length(rok)-floor(0.15*length(rok))])

To ją pokażę w porównaniu z fs. Oczywiście, jak dotychczas stosuję ustawienia graficzne okien za pomocą funkcji par, aby pokazać dwa wykresy obok siebie, a na koniec przywracam domyślne ustawienia:

fs = Fstats(stopa0 ~ 1 + stopa1 + reszty1)
par(mfrow=c(2,1), mar=c(2,5,2,2))
plot(stopa_gr, ylab="real GDP growth")
plot(fs, alpha=0.05)
par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Efekt:



Czyli mamy do czynienia z niestabilnością modelu MA(1). To sprawdźmy dokładnie, w którym roku nastąpiła zmiana struktury, stosując kod:

daty[as.numeric(breakpoints(fs)[1])]


Dostałem rok 1929. Czyli zauważmy, że QLR wskazał precyzyjnie, że w całym okresie badań PKB USA aż od 1797 r. powstało jedno pęknięcie - był to rok 1929. Dla przypomnienia QLR to supF, który wystarczył do wykazania zmiany; nie stosujemy już aveF ani expF, bo jest zbędne.

Trzeba zwrócić uwagę, że skoro KPSS nie wykrył żadnej zmiany, to nie nastąpiła realna zmiana średniej. A to logicznie znaczy, że nastąpiła zmiana w części MA. Tu z kolei dwie możliwości: albo nastąpiła zmiana autokorelacji albo wariancji. Patrząc na te stopy zmian, podejrzewam, że chodzi raczej o wariancję, bo właśnie w 1929 nastąpiło tąpnięcie. Do roku 1929 warunkowa wariancja ewidentnie rosła, osiągnęła max w 1929 i od tego roku zaczęła spadać.

Jeżeli jednak tak właśnie się stało, to powinniśmy posługiwać się przynajmniej modelem MA-GARCH, a nie zwykłym MA. To pokazuje, że ARMA nie jest wystarczające, a więc auto.arima także nie jest. Mamy też dowód, że zwykły ARMA nie może służyć do prognozy. Może jedynie służyć jako wstęp do dalszych badań. 


Przykład 8. Ropa naftowa (Crude Oil WTI): miesięczne stopy zwrotu (ostatnie 10 lat). Źródło ze stooq.pl.

#Wczytanie pliku csv:
nazwa = 'ropa.csv'

plik = read.csv(nazwa, sep=";")

if (ncol(plik)==1) {
  
  plik = read.csv(nazwa, sep=",")
  
}


daty = plik[,1]
daty = as.Date(daty)

#jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

if (daty[length(daty)]>Sys.Date()) {
  
  plik = plik[-nrow(plik),]
  
  daty = plik[,1]
  
  daty = as.Date(daty)
  
}


#stopy zwrotu
stopa = diff(cena)/cena[-length(cena)]
daty = daty[-1]

rok = as.numeric(format(daty, "%Y"))
mc = as.numeric(format(daty, "%m"))
stopa = ts(stopa, start=c(rok[1], mc[1]), frequency=12)
cena = ts(cena, start=c(rok[1], mc[1]), frequency=12)


Powtarzamy kod na KPSS:

#Analiza stóp zwrotu:
# test stacjonarności
if (length(stopa)>200) {
  kpss.test(stopa, lshort=FALSE)
}  else {
    kpss.test(stopa)
  }


Efekt:



Brak istotności.

Stabilność:
Ustawiam seasonal = FALSE, bo może się zdarzyć, że program wskaże sezonowość, na której jednak nie chcę się teraz skupiać.

arma_stopa = auto.arima(stopa,ic="bic", seasonal=FALSE)
arma_stopa 

Efekt:


auto.arima pokazał, że zgodnie BIC miesięczne stopy zwrotu nie są procesem ARMA, a mogą być białym szumem. Oczywiście jest to nieprawda, bo jest pewne, że ropa jest pewnym rodzajem modelu GARCH, a więc posiada zmienną wariancję. Na razie mamy jednak zwykłe błądzenie przypadkowe i dlatego testowany model będzie ARMA(0,0,0):

#testy stabilności
fs = Fstats(stopa ~ 1)
par(mfrow=c(2,1), mar=c(2,5,2,2))
cena_gr = window(cena, start=c(rok[ceiling(0.15*length(rok))], mc[ceiling(0.15*length(mc))]), end=c(rok[length(rok)-floor(0.15*length(rok))], mc[length(mc)-floor(0.15*length(mc))]))
plot(cena_gr, ylab="cena ropy")
plot(fs, alpha=0.05)
par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Efekt:


Jest to pierwszy test - supF. Wskazuje, że na początku 2020, a więc na początku pandemii, nastąpiła zmiana struktury modelu. Dokładną datę wyciągniemy dzięki:

daty[as.numeric(breakpoints(fs)[1])]

, co wskaże na kwiecień 2020.
Aby uzyskać p-value tego testu wpiszemy
sctest(fs, type="supF")[2]

Dostałem 0,0195. Podobnie jak poprzednio, nie używamy już aveF i expF.


Rodzi się tu istotne jedno pytanie. Skoro: 
- KPSS nie wykazał zmiany
- najlepszym modelem jest błądzenie losowe
- model jest niestabilny

to co to oznacza? Tutaj właśnie wkracza GARCH. Jeżeli (warunkowa) wariancja się zmienia, to KPSS tego nie uchwyci, natomiast testy stabilności tak, ponieważ założyliśmy błądzenie losowe - a GARCH nim nie jest. Widzimy więc, że dochodzimy do nowego problemu. Po pierwsze należy zbadać występowanie efektów (G)ARCH, a po drugie uwzględnić je w modelu i dopiero wtedy testować stabilność.  


Przykład 9.  sWIG80 - miesięczne stopy zwrotu od początku notowań do czerwca 2022 (12.1994-06.2022). Źródło: stooq.pl

Wszystko powtarzamy. Zwrócę może tylko uwagę, że przy badaniu stabilności ustawiam tutaj seasonal = FALSE, tak jak na ropie. Bez tego ustawienia program zdaje się szuka najpierw najlepszego modelu z sezonowością, mimo że bez niej mógłby być jeszcze lepszy. Gdyby jednak nie zmieniać tego ustawienia, to faktycznie program wskazał... sezonowość na małym WIGu. Można będzie się nad tym potem zastanowić, ale na razie nie chcę tego komplikować.

KPSS nie odrzuca stacjonarności, natomiast auto.arima pokazała, że najlepiej stosować AR(1) ze średnią zero (phi = 0,34 +/- 0,05). To zero trochę niepokoi, bo jednak powinna być niewielka (bo miesięczna) dodatnia stopa zwrotu z powodu awersji do ryzyka, szczególnie na małych spółkach. Może to sugerować pewne niedowartościowanie. Istotnie, obecne ROE na tym indeksie wynosi prawie 16%, a przecież taka rentowność dotyczyła S&P 500, gdzie korespondowała historycznie ze średnim C/WK aż 2,3. Tymczasem obecne C/WK naszych maluchów wynosi 1,3. Trzeba jednak pamiętać, że koszt kapitału samej Polski będzie nieco wyższy niż USA, a do tego należy dodać wyższy koszt kapitału małych i niepłynnych spółek. Ostatnie wydarzenia na Ukrainie, a także ryzyko polityczne (najniższe inwestycje prywatne od lat) dodatkowo hamują napływ kapitału. Jeśli spojrzeć na C/WK sprzed 6 miesięcy, to dostaniemy 1,7, czyli już poziom racjonalny.

W każdym razie zerowa średnia to zerowy wyraz wolny, czyli używamy modelu:

fs = Fstats(stopa0 ~ stopa1)

I wykres



Okazuje się, że wg supF model AR(1) dla swig pozostał stabilny. Najwyższa wartość pojawiła się w maju 2007, jednak daleko jej do przekroczenia czerwonej linii. W takim razie zróbmy wszystkie 3 testy:

sc1 = sctest(fs, type="supF")
p1 = as.numeric(sc1[2])
noquote(c("supF: p-value:", p1))

sc2 = sctest(fs, type="aveF")
p2 = as.numeric(sc2[2])
noquote(c("aveF: p-value:", p2))

sc3 = sctest(fs, type="expF")
p3 = as.numeric(sc3[2])
noquote(c("expF: p-value:", p3))


Funkcja noquote usuwa niepotrzebne cudzysłowy wyświetlanego tekstu przy stosowaniu print(). I tak po kolei dostałem:

supF: p-value: 0.518
aveF: p-value: 0.516
expF: p-value: 0.567

Zatem wszystkie 3 testy zgodnie podtrzymują hipotezę zerową o braku zmian w AR(1). Wygląda na to, że sWIG80 zachowuje się stabilnie -występuje stabilna autokorelacja 1 rzędu. Nie znaczy to jeszcze, że rynek jest nieefektywny (zob. Autokorelacja i efektywność rynku), ponieważ może się wiązać z niską całkowitą stopą zwrotu (w tym sporymi kosztami ekonomicznymi - czasochłonność, ryzyko), ale dostarcza pewną przewidywalność na kolejny miesiąc. 



[1] Zeileis, A., Leisch, F., Hornik, J., Kleiber, C., strucchange: An R Package for Testing for Structural Change in Linear Regression Models ;
[2]. Miller, D. M., Reducing Transformation Bias in Curve Fitting, May, 1984.