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ą.