niedziela, 3 lipca 2022

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

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

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

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

Część 1)

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


if (require(strucchange)==FALSE) {

  install.packages("strucchange")

  library("strucchange")

}

if (require(forecast)==FALSE) {

  install.packages("forecast")

  library(forecast)

}

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

if (require(tseries)==FALSE) {

  install.packages("tseries")

  library("tseries")

}


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



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

Wyniki:



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



Weźmy jeszcze przykładowo ziarno nr 20

    


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

Część 2)

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

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

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


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

Dla ziarna 20 wykres będzie podobny:


Część 3)

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

A) zmiana wyrazu wolnego (średniej)

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

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

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

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



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

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


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



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


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

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

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

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

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

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

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

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

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

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

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

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

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

moc = sukces/liczba_prób
moc

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

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

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

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

moc = sukces/liczba_prob
moc

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

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

B) zmiana części AR

Ograniczymy się do 1 rzędu. 

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

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

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

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

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





Sprawdźmy moc tego testu:

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

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

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

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

moc = sukces/liczba_prób
moc


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


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

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

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

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

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

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


C) zmiana części MA

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

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

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

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

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

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





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

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

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

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

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

moc = sukces/liczba_prob
moc

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


D) zmiana części AR i MA

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

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

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

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

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

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

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




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

I test mocy:

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

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

moc = sukces/liczba_prob
moc

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



Przeprowadźmy testy na prawdziwych danych.

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

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

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

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

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

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

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

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


Efekt:


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

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

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

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

Efekt:


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

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

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

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


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

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

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

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

Efekt:



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

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


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

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

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


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

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

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

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


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

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

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


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

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


Powtarzamy kod na KPSS:

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


Efekt:



Brak istotności.

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

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

Efekt:


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

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

Efekt:


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

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

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

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


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

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


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

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

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

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

fs = Fstats(stopa0 ~ stopa1)

I wykres



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

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

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

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


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

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

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



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

środa, 25 maja 2022

Symulacja i automatyczne modelowanie ARIMA w języku R

W języku R można symulować i modelować procesy ARIMA(p, d, q) za pomocą wielu pakietów, ale zaczniemy od funkcji wbudowanej w sam R. 

Symulacja

Symulacje wykonamy za pomocą funkcji arima.sim. 

Wybór ziarna

Do odtwarzania zawsze tych samych wyników na liczbach pseudolosowych służy tzw. ziarno, czyli wartość startowa, która determinuje cały ciąg liczb udających losowe. Możemy więc wybrać dowolną liczbę całkowitą i uzyskamy nowe ziarno, a tym samym inny przebieg losowy, który jednak możemy za każdym razem odtwarzać. Podstawowa funkcja to set.seed(x), gdzie x to dowolna liczba całkowita.  Czy liczba ziarn jest nieskończona? Nie. W pomocy do funkcji set.seed dowiemy się, że domyślnie używany jest stosunkowo nowy algorytm z 2002 r., którego okres wynosi 2^19937 - 1, czyli tyle jest unikatowych ziarn, a potem powtarzają się od nowa. Tyle w teorii. A w praktyce taka liczba ziarn jest niemożliwa dla programu, gdyż to oznaczałoby liczbę składającą się z 6002 cyfr. Z moich testów wynika, że obecna wersja programu (64bit) jest ograniczona do paru miliardów. W wersji której używam dokładnie jest 4 294 967 295 ziaren. Ponieważ argumenty funkcji seed() nie muszą być dodatnie, to tę liczbę należy rozbić:  2147483647*2 + 1, z czego 2147483647 dotyczy liczb dodatnich oraz ujemnych, stąd mnożymy razy 2, a dodanie 1 dotyczy argumentu 0. Oznacza to ni mniej ni więcej, że pierwsze ziarno otrzymamy za pomocą funkcji seed(-2147483647), a największe seed(2147483647). Możemy to sprawdzić sami. Stwórzmy próbkę z rozkładu normalnego o średniej 0 i wariancji 1:

set.seed(2147483646)
sim1 = rnorm(100, 0, 1)
plot(sim1, type="l")

To samo, ale zwiększmy seed o 1:

set.seed(2147483647)
sim1 = rnorm(100, 0, 1)
plot(sim1, type="l")

Wykres zmienił się. A teraz:

set.seed(2147483648)
sim1 = rnorm(100, 0, 1)
plot(sim1, type="l")

Otrzymamy informację, że wystąpił błąd. W zależności od wersji, którą posiadamy, program może się inaczej zachować i albo w ogóle nie zmienić wykresu, albo zresetować ziarno. To samo dotyczy liczb ujemnych. We wszystkich symulacjach musimy więc ograniczyć argument ziarna od -2147483647 do 2147483647. Jednak w poniższych przykładach będę używał zawsze tego samego ziarna - seed(1) - które przy zmienianych parametrach modelu i tak będzie przynosić zupełnie różne szeregi. 

Przykład 1: brak autokorelacji

AR(0) to proces czysto losowy o rozkładzie normalnym. Aby go zasymulować, stosowałem wcześniej funkcję rnorm. Dzięki ziarnu zobaczymy, że arima i rnorm dadzą dokładnie to samo. Przykładowo bierzemy pierwsze ziarno i generujemy 100 danych za pomocą arima.sim

set.seed(1)

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

plot(sim1, type="l")




Otrzymamy dokładnie ten sam wykres, gdy wykonamy kod z rnorm:

set.seed(1)

sim1 = rnorm(100,0,1)

plot(sim1, type="l")


Za pomocą funkcji acf i pacf możemy sprawdzić autokorelacje próby. O różnicy między ACF i PACF pisałem kiedyś w art. Jak sprawdzić stabilność modelu ekonometrycznego w czasie?. Zatem wiemy, że ACF służy do oceny stacjonarności, a PACF do oceny autokorelacji danego rzędu. W tym przypadku wpiszmy najpierw

acf(sim1)

Samo takie wpisanie wywoła wykres:

Zerowy rząd, który zawsze będzie równy 1 pomijamy, więc widać, że proces jest stacjonarny.

Potem 

pacf(sim1)

Brak autokorelacji - tak jak być powinno.

Przykład 2: autoregresja rzędu 1: AR(1). Powiedzmy, że współczynnik AR(1) - phi - wynosi 0,7. Stosujemy kod:

set.seed(1)

sim2 = arima.sim(list(order=c(1,0,0), ar=0.7), n=100)

plot(sim2, type="l")

Wykres:


Aby pokazać ACF i PACF obok siebie zastosujemy funkcję par(). Funkcja ustawia parametry grafiki w programie. Ponieważ robi to na stałe w sesji, to użyjemy jej dwukrotnie: raz do ustawienia dwóch wykresów, dwa do powrotu do ustawień domyślnych. Dodatkowo ustawimy tylko 5 rzędów acf i pacf, bo więcej nie potrzeba. To się robi za pomocą parametru lag.max:  

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

acf(sim2, lag.max=5)

pacf(sim2,lag.max=5)

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




Jak widać funkcja acf jest niespójna z pacf, bo w tej drugiej prawidłowo od razu zaczynamy od rzędu 1, a nie bez sensu od 0. W każdym razie ACF wskazuje na występowanie dwóch rzędów autokorelacji, co jednak nie przekłada się na liczbę rzędów AR, która zostaje wykryta dopiero przez pacf. Im silniejszy współczynnik AR, tym większa będzie liczba rzędów na ACF, nie wpływając na liczbę rzędów PACF. Gdy współczynnik przekroczy 1, proces staje się niestacjonarny. Popatrzmy na kolejny przykład.

Przykład 3. AR(1) z siłą phi = 1,3. Najpierw wpiszmy kod:

sim3 = arima.sim(list(order=c(1,0,0), ar=1.3), n=100)

Wynik:



Otrzymaliśmy błąd, ponieważ AR nie może być niestacjonarny w modelu ARIMA. Aby uzyskać AR(1), musimy go sprowadzić do stacjonarności i do tego właśnie służy drugi parametr, d. Otrzymamy model autoregresyjny rzędu 1 z integracją stopnia 1, tj. ARI(1, 1) lub ogólniej ARIMA(1, 1, 0). Ze względu na integrację zamiast siły 1,3 wpisujemy różnicę o 1, czyli 0,3:

sim3 = arima.sim(list(order=c(1,1,0), ar=0.3), n=100)

plot(sim3, type="l")

I wykres będzie typowo rynkowy:


Po wykonaniu jeszcze raz kodu:

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

acf(sim3, lag.max=5)

pacf(sim3,lag.max=5)

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

 dostaniemy coś ciekawego:


ACF wskazuje na niestacjonarność (powolne spadki autokorelacji), podczas gdy PACF informuje nadal tylko o 1 rzędzie. Można by zapytać dlaczego, skoro mamy tu do czynienia z procesem niestacjonarnym? No właśnie dzieje się tak, bo to co tutaj mamy to AR(1) z p >= 1, więc całkowita autokorelacja wskazuje na bardzo dużą liczbę rzędów (cena z dawnej przeszłości koreluje z ceną teraźniejszą), ale żeby dostać autokorelację tylko dla jednego rzędu, trzeba usunąć wszystkie inne, które na nią wpływają. Tak działa PACF. Czyli zauważmy, że mimo iż przeszliśmy z procesu stacjonarnego do niestacjonarnego tylko ACF się bardzo zmieniło, a PACF zmieniło jedynie poziom samej autokorelacji rzędu 1.

Przykład 4: średnia ruchoma rzędu 1: MA(1). Powiedzmy, że współczynnik MA(1) - theta  - wynosi 0,7. Stosujemy kod:

set.seed(1)

sim4 = arima.sim(list(order=c(0,0,1), ma=0.7), n=100)

plot(sim4, type="l")


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

acf(sim4, lag.max=5)

pacf(sim4,lag.max=5)

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





Po pierwsze - mimo że ciągle używamy seed(1), to wynik znowu jest inny dla MA(1) niż dla AR(1). Po drugie zmieniła się struktura PACF - korelacja 1-go rzędu jest nadal wysoka, ale pojawiła się drugiego i czwartego ujemna. To świadczy o tym, że MA jest dużo mniej intuicyjna niż AR, ale z drugiej strony obejmuje bardziej skomplikowane procesy. Fakt, że mamy tu ujemną autokorelację 2-go rzędu nie jest przypadkiem. Dla lepszej ilustracji oraz dalszej analizy zróbmy to samo, ale zwiększmy próbę do 1000 oraz liczbę rzędów autokorelacji do 15. Za n wstawiamy do kodu 1000, reszta taka sama. Efekt:



Jak się okazuje PACF realizuje się w postaci na przemian dodatniej i ujemnej autokorelacji. Z czego to wynika? Popatrzmy, MA(1) przy zerowym wyrazie wolnym ma postać:

(1)


W okresie 0 składnik losowy z okresu -1 jest równy zero i wtedy:



Dla okresu 1 podstawiamy do (1) i wyznaczamy składnik losowy:




W okresie 2 i 3 analogicznie:




W końcu w okresie 4 na podstawie poprzednich danych już tylko wyznaczamy y(4):




Zauważmy, że dostaliśmy właśnie na przemian znaki plus i minus dla kolejnych rzędów zależnych od y z okresów poprzednich. Ten wzór będzie się powtarzał do nieskończoności, a to znaczy, że MA(1) można zapisać w postaci procesu AR, w którym następuje jakby częściowe "wygaszanie" wpływu poprzednich obserwacji. 

MA jest procesem ergodycznym i stacjonarnym (w przeciwieństwie do AR, który być nie musi). Zarówno ergodyczność jak i stacjonarność wiążą się prawem regresji do średniej. Sama stacjonarność nie wystarcza dla niego, bo może być proces stacjonarny bez żadnych zmian, jakaś stała [1]. Dopiero ergodyczność wymusza zmienność, gdyż wiąże się z rozkładem częstości. Ergodyczność zapewnia, że prędzej czy później zmienna zapełni cały swój rozkład, tak że zawsze uzyskamy np. kształt dzwonu. W każdym razie zawsze dostaniemy taki obraz PACF dla MA(1) z powodu powrotu do średniej, która zostaje jakby wzmocniona zarówno w stosunku do AR jak i procesu czysto losowego.

Można by jeszcze zapytać dlaczego MA nazywana jest średnią ruchomą albo średnią kroczącą. Pomijając względy historyczne, bo nie jest to rzeczywiście średnia, to zauważmy, że:

* dostaliśmy wyżej model zależny liniowo od wartości poprzednich, a współczynniki można porównać z wagami średniej ważonej (chociaż nie muszą się sumować do 1, dlatego to nie jest prawdziwa średnia),

* mimo że występuje zamiana znaków plus na minus, co wydaje się psuć ideę średniej kroczącej, to zauważmy, że MA dotyczy zmiennej stacjonarnej, a więc zawsze przyspiesza prawo regresji do średniej. Co więcej, możemy utworzyć MA kolejnych rzędów. To znaczy, że dla MA(2) drugi rząd będzie się zachowywał tak samo jak pierwszy z przesunięciem o 1, a to oznacza, że dla okresu 2 znak minus z MA(1) będzie jakby redukowany przez plus z MA(2). 


Przykład 5. ARMA(1, 0, 1). phi(1) = 0,7 i theta(1) = 0,7. Kod:

set.seed(1)

sim5 = arima.sim(list(order=c(1,0,1), ar=0.7, ma=0.7), n=1000)

plot(sim5, type="l")



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

acf(sim5, lag.max=15)

pacf(sim5,lag.max=15)

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






Znowu możemy dokładnie obserwować zmiany procesu, bo ziarno nie zmieniło się. Dostaliśmy mix AR i MA, dzięki czemu proces ma własności jednego i drugiego.

Przykład 6. ARIMA(1, 1, 1). Siła AR(1) = 0,7 i siła MA(1) = 0,7 oraz d = 1. Kod:

set.seed(1)

sim6 = arima.sim(list(order=c(1,1,1), ar=0.7, ma=0.7), n=1000)

plot(sim6, type="l")



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

acf(sim6, lag.max=15)

pacf(sim6,lag.max=15)

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



Wprowadzenie jednoczesne AR i MA wywołało autokorelację rzędu 1 praktycznie równą jedności. I co więcej, nawet jeśli zmniejszymy parametr MA niewiele się zmieni. W ogóle wykres się praktycznie nie zmieni. Zobaczymy to w kolejnym przykładzie.

Przykład 7. ARIMA(1, 1, 1): phi(1) = 0,7 i theta(1) = 0,1 oraz d = 1. Kod:

set.seed(1)

sim7 = arima.sim(list(order=c(1,1,1), ar=0.7, ma=0.1), n=1000)

plot(sim7, type="l")

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

acf(sim7, lag.max=15)

pacf(sim7,lag.max=15, plot=FALSE)

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





Przykład 8Ciekawszy wynik dostaniemy, jeśli ustawimy phi na 0.1 thetę na -0.95. Funkcja PACF występuje dwukrotnie, pierwsza pokazuje wykres, druga - same współczynniki korelacji:

set.seed(1)

sim8 = arima.sim(list(order=c(1,1,1), ar=0.1, ma=-0.95), n=1000)

plot(sim8, type="l")


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

acf(sim8, lag.max=15)

pacf(sim8,lag.max=15)

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





Mimo iż wykres stał się bardzo poszarpany, to nadal występuje dodatnia korelacja  i to wielu rzędów. Nawet bardzo ujemna część MA nie przykrywa małej dodatniej części AR.


Optymalne modelowanie

Kolejnym etapem jest fitting, czyli dopasowanie modelu. Najłatwiej użyć optymalizacji za pomocą auto.arima(). Jest to funkcja zawarta w pakiecie forecast, który należy zainstalować:
install.packages("forecast")
i włączyć
library(forecast)

Sprawdzimy, jak auto.arima radzi sobie ze wszystkimi podanymi wyżej przykładami. 

Ad 1: brak autokorelacji. Dla sim1:

dopas=auto.arima(sim1)

dopas

Otrzymamy:

Series: sim1
ARIMA(0,0,0) with zero mean


Prawidłowo.

Ad 2: phi(1) = 0,7. Za sim1 wstawiamy sim2, więc nie ma sensu powtarzać kodu. Wyniki:

Series: sim2
ARIMA(1,0,0) with zero mean

Coefficients:
ar1
0.6956
s.e. 0.0712

Prawidłowo.

Ad 3phi = 1,3.

Series: sim3
ARIMA(1,1,0)


Coefficients:
ar1
0.3116
s.e. 0.0958

Prawidłowo.


Ad 4theta = 0.7

Series: sim4 

ARIMA(0,0,2) with zero mean 

Coefficients:

         ma1      ma2

      0.6609  -0.0724

s.e.  0.0318   0.0345


Pierwszy raz błąd. Powinno być MA(1), a dostaliśmy MA(2), pomimo, iż było 1000 danych. Gdy wrócimy do PACF dla przykładu 4, to zobaczymy, że występują istotne ujemne autokorelacje. Jest to niewątpliwie przyczyna błędu. Mimo wszystko powinniśmy uzyskać optimum dla MA(1). Co w takiej sytuacji zrobić? Warto przyjrzeć się p-value. Niestety pakiet tego nie oferuje. Żeby je uzyskać możemy użyć kodu:

(1-pnorm(abs(dopas$coef)/sqrt(diag(dopas$var.coef))))*2

Wynik:

       ma1        ma2 

0.00000000 0.03583408  

p-value dla ma1 = 0 i ma2 = 0.0358. Nie można więc odrzucić modelu MA(2). 

Tak więc już napotykamy pierwsze problemy. Pakiet auto.arima domyślnie głównie opiera się na AIC, ale wg moich testów BIC lepiej się sprawdza dla ARIMA. Algorytm można jednak ustawić tak, aby kierował się tylko BIC. Kod musimy zmodyfikować. Domyślny kod miał ustawienia następujące:

dopas=auto.arima(sim1, ic = c("aicc", "aic", "bic"))

Zastąpimy go:

dopas=auto.arima(sim1, ic = c("bic"))

dopas

Wynik:

Series: sim4 

ARIMA(0,0,1) with zero mean 

Coefficients:

         ma1

      0.7039

s.e.  0.0258


Okazuje się, że kryterium BIC prawidłowo pokazało, że najlepszym modelem jest MA(1).


Ad 5. ARIMA(1, 0, 1). phi(1) = 0,7 i theta(1) = 0,7. Ponowne zastosowanie z BIC daje wyniki

ARIMA(1,0,1) with zero mean 

Coefficients:

         ar1     ma1

      0.6404  0.7329

s.e.  0.0266  0.0268


Czyli prawidłowo. Dla domyślnych ustawień ponownie algorytm nie działał dobrze, dając ARIMA(2, 0, 1).

Ad 6. ARIMA(1, 1, 1). phi(1) = 0,7 i theta(1) = 0,7. To samo, ok.

Ad 7. ARIMA(1, 1, 1). Siła AR(1) = 0,7 i siła MA(1) = 0,1. Dla BIC:

Series: sim7 

ARIMA(1,1,0) 

Coefficients:

         ar1

      0.7030

s.e.  0.0225


Tym razem BIC nie poradził sobie, najwyraźniej theta = 0,1 była za małą wartością. Ale phi jest właściwa. Sprawdźmy jeszcze domyślne kryteria:

Series: sim7 

ARIMA(2,1,2) 

Coefficients:

          ar1     ar2     ma1     ma2

      -0.2678  0.6085  1.0397  0.0686

s.e.   0.0485  0.0346  0.0554  0.0470


Domyślny algorytm poradził sobie jeszcze gorzej.


Ad 8. ARIMA(1, 1, 1) z phi 0.1, theta -0.95. Zarówno domyślny jak i kryterium BIC dało to samo:

Series: sim8 

ARIMA(1,1,2) 

Coefficients:

          ar1      ma1      ma2

      -0.9124  -0.0011  -0.9222

s.e.   0.0401   0.0304   0.0285


Obydwa kryteria błędnie wskazały zarówno liczbę rzędów, jak i wartości współczynników. Proces okazał się dla nich zbyt złożony. Zwiększmy więc liczbę danych z 1000 na 10000:

set.seed(1)

sim8 = arima.sim(list(order=c(1,1,1), ar=0.1, ma=-0.95), n=10000)

plot(sim8, type="l")




Co się okazuje, obydwa kryteria pokazują to samo i radzą sobie dość dobrze:

Series: sim8 

ARIMA(2,1,1) 

Coefficients:

         ar1     ar2      ma1

      0.1205  0.0157  -0.9582

s.e.  0.0105  0.0105   0.0032

Chociaż model błędnie dodaje phi(2), to jego p-value okazuje się nieistotne (powtarzamy kod na p-value):

(1-pnorm(abs(dopas$coef)/sqrt(diag(dopas$var.coef))))*2

      ar1       ar2       ma1 

0.0000000 0.1331158 0.0000000 

Wobec tego odrzucamy ar2, tak że automat znajduje prawidłowe rozwiązanie.


Podsumowując część optymalnego modelowania, dostajemy dowód, że w ARIMA należy posługiwać się BIC - przynajmniej stosując auto.arima tak jak wyżej. Oczywiście są jeszcze inne możliwe kryteria do sprawdzenia, ale pakiet oferuje jedynie AIC i BIC. 


Literatura:

[1] Hansen, B., Econometrics, 2022, link: https://www.ssc.wisc.edu/~bhansen/econometrics/

P. S. Bardzo dobre, zrozumiałe opracowanie.