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
