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.
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
# 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
B) z autokorelacją
sim1 = arima.sim(list(order=c(1,0,1), ar = 0.5, ma=0.00001), n=dlugosc_proby)
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
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:
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.
-----------------------------------------------------------------------------------------
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)
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
Brak komentarzy:
Prześlij komentarz