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

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 (tu i tu).

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.