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.