Model GARCH jako taki nie służy do prognozowania cen czy stóp zwrotu, ale ich zmienności. Można go jednak rozumieć w szerszy sposób, tj. ARIMA-GARCH. Ten artykuł jest kontynuacją wcześniejszych: Dlaczego oddzielne badanie autokorelacji jest błędem? oraz Modelowanie ryzyka z rugarch (język R). Omówię te dwa sposoby z punktu widzenia trejdingu.
Przygotowanie danych
Na początek ładuję dwa pakiety: xts oraz rugarch
# Przedwstęp
if (require("xts")==FALSE) {
install.packages("xts")
}
library("xts")
if (require("rugarch")==FALSE) {
install.packages("rugarch")
}
library("rugarch")
Analizę wykonam na WIG20, plik .csv - dzienne dane - pobrane ze stooq.pl. Przekształcam je w tygodniowe i ucinam do 200 ostatnich notowań (ostatnie to 13.09.2024). Przekształcam w stopy zwrotu.
# Wstęp
nazwaPliku = "wig20_d.csv"
Obs <- 200
plik <- read.csv(nazwaPliku)
plik$Data <- as.Date(plik$Data)
tabela <- tail(to.weekly(as.xts(plik), indexAt="endof", OHLC=FALSE), Obs)
cena <- tabela$Zamkniecie
stopaZwrotu <- diff.xts(log(cena))
tabela$stopaZwrotu <- stopaZwrotu
tabela <- na.omit(tabela)
stopaZwrotu <- tabela$stopaZwrotu
1. ARIMA-GARCH
Szersze ujęcie GARCH dzieli się na dwie części: najpierw szukamy optymalnego ARIMA za pomocą funkcji autoarfima(). Ustawiam max 10 rzędów zarówno dla części AR, jak i MA. Na początek zakładam rozkład skośno-normalny głównie z powodu szybkości obliczania jego parametrów:
# 1. część ARMA
arimaDopas <- autoarfima(stopaZwrotu, ar.max = 10, ma.max = 10, criterion = "HQIC",
distribution.model = "snorm",
method="partial")
Po wpisaniu zmiennej dostajemy statystyki
arimaDopas
*----------------------------------* * ARFIMA Model Fit * *----------------------------------* Mean Model : ARFIMA(5,0,3) Distribution : snorm Optimal Parameters ------------------------------------ Estimate Std. Error t value Pr(>|t|) mu 0.002644 0.000001 2376.47 0 ar1 -0.616537 0.000031 -19828.09 0 ar2 1.054657 0.000034 31037.94 0 ar3 0.533467 0.000028 19356.08 0 ar4 -0.173088 0.000010 -17629.16 0 ar5 0.053956 0.000006 8470.41 0 ma1 1.038622 0.000009 118359.33 0 ma2 -1.083726 0.000008 -139217.06 0 ma3 -1.131477 0.000060 -18831.84 0 sigma 0.024462 0.000096 254.15 0 skew 0.833062 0.002360 352.95 0 Robust Standard Errors: Estimate Std. Error t value Pr(>|t|) mu 0.002644 0.000000 6484.164 0 ar1 -0.616537 0.000119 -5183.800 0 ar2 1.054657 0.000090 11723.967 0 ar3 0.533467 0.000040 13248.164 0 ar4 -0.173088 0.000041 -4179.572 0 ar5 0.053956 0.000019 2915.708 0 ma1 1.038622 0.000358 2897.527 0 ma2 -1.083726 0.000012 -88962.282 0 ma3 -1.131477 0.000060 -18854.153 0 sigma 0.024462 0.000949 25.785 0 skew 0.833062 0.003603 231.216 0 LogLikelihood : 456.9 Information Criteria ------------------------------------ Akaike -4.4819 Bayes -4.2999 Shibata -4.4876 Hannan-Quinn -4.4082 Weighted Ljung-Box Test on Standardized Residuals ------------------------------------ statistic p-value Lag[1] 2.763 0.09649 Lag[2*(p+q)+(p+q)-1][23] 11.182 0.91984 Lag[4*(p+q)+(p+q)-1][39] 17.964 0.71904 H0 : No serial correlation Weighted Ljung-Box Test on Standardized Squared Residuals ------------------------------------ statistic p-value Lag[1] 2.091 0.1482 Lag[2*(p+q)+(p+q)-1][2] 2.238 0.2262 Lag[4*(p+q)+(p+q)-1][5] 4.670 0.1816 ARCH LM Tests ------------------------------------ Statistic DoF P-Value ARCH Lag[2] 2.347 2 0.3093 ARCH Lag[5] 9.088 5 0.1056 ARCH Lag[10] 13.467 10 0.1987 Nyblom stability test ------------------------------------ Joint Statistic: 7.683 Individual Statistics: mu 0.00450 ar1 0.21211 ar2 0.20711 ar3 0.19583 ar4 0.09643 ar5 0.12618 ma1 0.19549 ma2 0.19475 ma3 0.16916 sigma 0.27608 skew 0.10861 Asymptotic Critical Values (10% 5% 1%) Joint Statistic: 2.49 2.75 3.27 Individual Statistic: 0.35 0.47 0.75
Najlepszy to ARMA(5, 3). Zapisujemy optymalne rzędy:
rzadAR <- arimaDopas$rank.matrix[1,1]
rzadMA <- arimaDopas$rank.matrix[1,2]
Model jest niestabilny (Joint Statistics > 2,75). Prawdopodobnie wybraliśmy zły rozkład. Zmieńmy na skośny uogólniony rozkład błędów, tj. "sged". Ponieważ wybór innego rozkładu niż normalny i skośno-normalny wydłuża znacząco cały proces obliczeń, dodamy parametr cluster, który wykorzystuje pakiet parallel do równoczesnego przetwarzania danych przez pozostałe rdzenie procesora (na koniec wyłączamy tę opcję); to przyspieszy proces:
klaster <- makeCluster(detectCores() - 1)
arimaDopas <- autoarfima(stopaZwrotu, ar.max = 10, ma.max = 10, criterion = "HQIC",
distribution.model = "sged", cluster = klaster,
method="partial")
stopCluster(klaster)
arimaDopas
*----------------------------------*
* ARFIMA Model Fit *
*----------------------------------*
Mean Model : ARFIMA(5,0,3)
Distribution : sged
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.002994 0.000001 2333.1750 0
ar1 -0.735343 0.000231 -3185.0574 0
ar2 1.065479 0.000335 3185.0574 0
ar3 0.715202 0.000225 3185.0574 0
ar4 -0.164068 0.000052 -3185.0574 0
ar5 -0.012163 0.000005 -2435.7292 0
ma1 1.043899 0.000029 36290.9045 0
ma2 -1.073710 0.000050 -21428.3317 0
ma3 -1.130274 0.000008 -139344.0095 0
sigma 0.024216 0.001212 19.9784 0
skew 0.872133 0.075043 11.6218 0
shape 1.711239 0.172880 9.8984 0
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.002994 0.000005 556.1676 0.000000
ar1 -0.735343 0.021895 -33.5850 0.000000
ar2 1.065479 0.053812 19.8000 0.000000
ar3 0.715202 0.020475 34.9306 0.000000
ar4 -0.164068 0.001227 -133.6667 0.000000
ar5 -0.012163 0.000004 -2972.8839 0.000000
ma1 1.043899 0.000706 1479.4296 0.000000
ma2 -1.073710 0.001064 -1009.1810 0.000000
ma3 -1.130274 0.000055 -20500.1353 0.000000
sigma 0.024216 0.001543 15.6918 0.000000
skew 0.872133 0.112199 7.7731 0.000000
shape 1.711239 0.705025 2.4272 0.015216
LogLikelihood : 459.4
Information Criteria
------------------------------------
Akaike -4.4965
Bayes -4.2979
Shibata -4.5033
Hannan-Quinn -4.4161
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.002713 0.9585
Lag[2*(p+q)+(p+q)-1][23] 8.977574 1.0000
Lag[4*(p+q)+(p+q)-1][39] 16.731900 0.8385
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.5398 0.4625
Lag[2*(p+q)+(p+q)-1][2] 0.6712 0.6184
Lag[4*(p+q)+(p+q)-1][5] 3.1718 0.3766
ARCH LM Tests
------------------------------------
Statistic DoF P-Value
ARCH Lag[2] 0.8306 2 0.6601
ARCH Lag[5] 7.7465 5 0.1708
ARCH Lag[10] 12.4933 10 0.2534
Nyblom stability test
------------------------------------
Joint Statistic: 1.744
Individual Statistics:
mu 0.05618
ar1 0.03768
ar2 0.03560
ar3 0.03446
ar4 0.03642
ar5 0.03377
ma1 0.06555
ma2 0.06355
ma3 0.07824
sigma 0.24067
skew 0.25045
shape 0.27590
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 2.69 2.96 3.51
Individual Statistic: 0.35 0.47 0.75Rzeczywiście model staje się stabilny, co oznacza, że problemem był rozkład. Nadal najlepszy jest ARMA(5, 3). Zaskakujące, że diagnostyka nie wykrywa ani autokorelacji w kwadratach reszt, ani efektu ARCH.
2. GARCH
Skoro diagnostyka nie wychwyciła występowania ARCH, można by pominąć ten krok i zająć się od razu prognozą. Jednak zastosowanie modelu wariancji pozwoli sprawdzić również efekt Arch-In-Mean. W tym wpisie pokazałem, że dzienne stopy zwrotu mogą być modelowane za pomocą TGARCH-ARCH-In-Mean, a miesięczne za pomocą EGARCH-ARCH-In-Mean. Sugeruje to, że tygodniowe będą czymś pomiędzy. Dla przypomnienia model ARCH-in-Mean oznacza, że stopa zwrotu jest zależna od wariancji warunkowej. Wówczas jednak rzędy ARMA były pominięte. Co by się stało, gdybyśmy je dodali? Oczywiście możemy to zrobić poprzez dodanie rzędów ARMA(5, 3), ale wcale nie jest pewne, że będzie to model najlepszy. Dlatego warto zrobić własną funkcję szukającą najlepszego modelu ARMA-GARCH. Funkcja wykorzysta otrzymane z autoarfimy optymalne rzędy jako maksymalne rzędy przeszukiwania, tzn. będziemy sprawdzać rzędy od 0 do rzadAR oraz od 0 do rzadMA. Co więcej będzie uwzględniał zarówno Arch-In-Mean, jak i bez niego. Kod dla TGARCH ma postać:
# 2. część GARCH
kryteriaInfor <- list()
for (k in 0:1) {
for (i in 0:rzadAR) {
for (j in 0:rzadMA) {
mojGarch <- ugarchspec(variance.model = list(model = "fGARCH", submodel = "TGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(i, j), archm = as.logical(k)),
distribution.model = "sged")
dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu)
nr <- paste0("ARMA(", i, ", ", j, ") - ARCH In Mean")
# usuwamy ARCH in mean dla k = 1
nr <- ifelse(k == 1, nr, sub(" - ARCH In Mean", "", nr))
kryteriaInfor[[nr]] <- infocriteria(dopasGarch)
}
}
}
kryteriaDf <- do.call(cbind, kryteriaInfor)
colnames(kryteriaDf) <- names(kryteriaInfor)
minKryteria <- apply(kryteriaDf, 1, function(wiersz) {
wartoscMin = min(wiersz)
kolumnaMin = colnames(kryteriaDf)[which.min(wiersz)]
return(c(Minimum = wartoscMin, Kolumna = kolumnaMin))
})
minKryteriaDf <- as.data.frame(t(minKryteria))
minKryteriaDf
Otrzymany wynik:
Minimum Kolumna
Akaike -4.45083487203567 ARMA(5, 3)
Bayes -4.2045670547213 ARMA(3, 3)
Shibata -4.46117146205014 ARMA(5, 3)
Hannan-Quinn -4.35036611102375 ARMA(5, 3)
Oznacza to, że wynik autoarfimy koresponduje z TGARCH. Dla EGARCH kod jest oczywiście taki sam poza fragmentem:
list(model = "fGARCH", submodel = "TGARCH"
który zamienia się w:list(model = "eGARCH", submodel = NULL
I z niego otrzymałem:
Minimum Kolumna
Akaike -4.43263086741263 ARMA(2, 2)
Bayes -4.25058889217661 ARMA(2, 2)
Shibata -4.43832575443847 ARMA(2, 2)
Hannan-Quinn -4.3589537760039 ARMA(2, 2)Sytuacja okazuje się niejednoznaczna. Z punktu widzenia AIC TGARCH jest lepszy. Z punktu widzenia BIC i HQIC lepszy jest EGARCH. Jeśli chodzi o SHC, czyli kryterium Shibata, to ja to pomijam, po pierwsze niemal zawsze wskazuje ten sam kierunek co AIC, a po drugie jest słabo opisany w literaturze.
Należy zauważyć, że oba modele nie zawierają elementu Arch-In-Mean. Można powiedzieć, że ARMA go zastąpiła.
Otrzymane modele wymagają diagnostyki. Zacznijmy od ARMA(5, 3)-TGARCH(1, 1):
mojGarch <- ugarchspec(variance.model = list(model = "fGARCH", submodel = "TGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(rzadAR, rzadMA)), distribution.model = "sged")
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : fGARCH(1,1)
fGARCH Sub-Model : TGARCH
Mean Model : ARFIMA(5,0,3)
Distribution : sged
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.001222 0.000004 287.4520 0.000000
ar1 -0.646631 0.000094 -6851.4518 0.000000
ar2 1.106298 0.000084 13106.6297 0.000000
ar3 0.627234 0.000090 6953.5872 0.000000
ar4 -0.102406 0.000040 -2533.7961 0.000000
ar5 0.099523 0.000039 2578.2384 0.000000
ma1 1.035966 0.000078 13332.0845 0.000000
ma2 -1.079756 0.000171 -6301.9226 0.000000
ma3 -1.128152 0.000256 -4412.9976 0.000000
omega 0.002007 0.001332 1.5071 0.131781
alpha1 0.054355 0.024212 2.2450 0.024768
beta1 0.876710 0.053657 16.3391 0.000000
eta11 0.996087 0.508698 1.9581 0.050217
skew 0.864175 0.062759 13.7697 0.000000
shape 1.601437 0.251582 6.3655 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.001222 0.000018 66.728331 0.00000
ar1 -0.646631 0.004215 -153.419278 0.00000
ar2 1.106298 0.007908 139.895269 0.00000
ar3 0.627234 0.004506 139.206023 0.00000
ar4 -0.102406 0.000244 -419.427330 0.00000
ar5 0.099523 0.000060 1667.478136 0.00000
ma1 1.035966 0.000476 2177.437996 0.00000
ma2 -1.079756 0.000322 -3357.818162 0.00000
ma3 -1.128152 0.004830 -233.579659 0.00000
omega 0.002007 0.006530 0.307381 0.75855
alpha1 0.054355 0.686095 0.079224 0.93685
beta1 0.876710 0.741582 1.182216 0.23712
eta11 0.996087 8.183445 0.121720 0.90312
skew 0.864175 0.042192 20.481868 0.00000
shape 1.601437 0.680752 2.352454 0.01865
LogLikelihood : 457.9
Information Criteria
------------------------------------
Akaike -4.4508
Bayes -4.2026
Shibata -4.4612
Hannan-Quinn -4.3504
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.3174 0.5732
Lag[2*(p+q)+(p+q)-1][23] 10.3414 0.9983
Lag[4*(p+q)+(p+q)-1][39] 18.5933 0.6475
d.o.f=8
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.00004328 0.9948
Lag[2*(p+q)+(p+q)-1][5] 3.39042514 0.3404
Lag[4*(p+q)+(p+q)-1][9] 5.51638020 0.3576
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 4.431 0.500 2.000 0.03530
ARCH Lag[5] 5.177 1.440 1.667 0.09349
ARCH Lag[7] 6.140 2.315 1.543 0.13209
Nyblom stability test
------------------------------------
Joint Statistic: 3.741
Individual Statistics:
mu 0.01584
ar1 0.01865
ar2 0.02030
ar3 0.01897
ar4 0.02063
ar5 0.01946
ma1 0.01911
ma2 0.01778
ma3 0.01841
omega 0.06211
alpha1 0.04989
beta1 0.08077
eta11 0.31789
skew 0.07723
shape 0.04888
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 3.26 3.54 4.07
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.9408 0.3480
Negative Sign Bias 0.3877 0.6987
Positive Sign Bias 0.2705 0.7871
Joint Effect 1.0407 0.7914
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 18.99 0.4575
2 30 34.82 0.2105
3 40 44.82 0.2409
4 50 40.95 0.7865Model jest lekko niestabilny (przez to parametry (G)ARCH nieistotne). Teraz ARMA(2, 2)-EGARCH(1,1).
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : eGARCH(1,1)
Mean Model : ARFIMA(2,0,2)
Distribution : sged
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.000872 0.000007 1.3068e+02 0.0e+00
ar1 -1.695986 0.000000 -1.6489e+07 0.0e+00
ar2 -0.735113 0.000000 -9.7011e+06 0.0e+00
ma1 1.874230 0.000648 2.8906e+03 0.0e+00
ma2 0.868271 0.001679 5.1706e+02 0.0e+00
omega -0.384768 0.000038 -1.0045e+04 0.0e+00
alpha1 -0.189127 0.046551 -4.0628e+00 4.8e-05
beta1 0.947281 0.001075 8.8104e+02 0.0e+00
gamma1 -0.075395 0.000542 -1.3914e+02 0.0e+00
skew 0.836575 0.112010 7.4688e+00 0.0e+00
shape 2.646523 0.666196 3.9726e+00 7.1e-05
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.000872 0.000542 1.6082e+00 0.107795
ar1 -1.695986 0.000013 -1.3200e+05 0.000000
ar2 -0.735113 0.000008 -8.7846e+04 0.000000
ma1 1.874230 0.091483 2.0487e+01 0.000000
ma2 0.868271 0.315734 2.7500e+00 0.005959
omega -0.384768 0.107277 -3.5867e+00 0.000335
alpha1 -0.189127 4.249978 -4.4501e-02 0.964505
beta1 0.947281 0.106032 8.9339e+00 0.000000
gamma1 -0.075395 0.025790 -2.9234e+00 0.003462
skew 0.836575 0.494241 1.6926e+00 0.090523
shape 2.646523 61.861149 4.2782e-02 0.965876
LogLikelihood : 452.0468
Information Criteria
------------------------------------
Akaike -4.4326
Bayes -4.2506
Shibata -4.4383
Hannan-Quinn -4.3590
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 1.428 0.2320
Lag[2*(p+q)+(p+q)-1][11] 3.101 1.0000
Lag[4*(p+q)+(p+q)-1][19] 9.156 0.6119
d.o.f=4
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.1386 0.7097
Lag[2*(p+q)+(p+q)-1][5] 0.8697 0.8884
Lag[4*(p+q)+(p+q)-1][9] 2.8003 0.7914
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 0.2651 0.500 2.000 0.6066
ARCH Lag[5] 0.8789 1.440 1.667 0.7692
ARCH Lag[7] 2.6093 2.315 1.543 0.5901
Nyblom stability test
------------------------------------
Joint Statistic: 2.4751
Individual Statistics:
mu 0.03485
ar1 0.03350
ar2 0.03469
ma1 0.04462
ma2 0.03174
omega 0.15840
alpha1 0.14496
beta1 0.16233
gamma1 0.20550
skew 0.09238
shape 0.09576
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 2.49 2.75 3.27
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.6633 0.5080
Negative Sign Bias 0.1653 0.8689
Positive Sign Bias 0.3918 0.6956
Joint Effect 0.5216 0.9141
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 15.57 0.6856
2 30 39.04 0.1009
3 40 39.99 0.4258
4 50 54.52 0.2728EGARCH jest dużo bardziej stabilny od TGARCH, ale efekt GARCH tylko częściowo występuje (alpha1 nieistotna, beta1 istotna). W tym miejscu powinniśmy zadać sobie pytanie, czy w ogóle nie należy pominąć GARCH i nie wybrać samego modelu ARMA(5, 3), który został uznany za najlepszy przez autoarfimę, a jej diagnostyka w ogóle nie wykryła autokorelacji w wariancjach składnika losowego. Dla ARMA(5, 3) dostaliśmy wartości kryteriów
Akaike -4.4965
Bayes -4.2979
Shibata -4.5033
Hannan-Quinn -4.4161Dla ARMA(5, 3)-TGARCH(1,1) dostaliśmy
Akaike -4.4508
Bayes -4.2026
Shibata -4.4612
Hannan-Quinn -4.3504Dla ARMA(2, 2)-EGARCH(1, 1):
Akaike -4.4326
Bayes -4.2506
Shibata -4.4383
Hannan-Quinn -4.3590Wydawałoby się, że należy po prostu zastosować zwykły ARMA(5, 3). Jednakże takie porównanie 1:1 nie wydaje się rozsądne. W przypadku modelu ARIMA prognozujemy jedynie stopy zwrotu, podczas gdy w GARCH (czy precyzyjniej ARIMA-GARCH) nie tylko stopy zwrotu, ale i odchylenie od średniej.
W sumie należy oddzielić dokładność prognozy stopy zwrotu od prognozy jej odchylenia. Zwykła korelacja nie nadaje się, bo skupia się na kierunku zmian, a nie ich znaku. Gdyby nasza prognoza miała zawsze taki sam znak jak prawdziwa stopa zwrotu, ale różny kierunek (tzn. gdy prawdziwa zmiana rośnie, to prognoza spada), to dostalibyśmy ujemną korelację, a więc zupełnie nie ma to sensu. Dlatego potrzebna jest miara badająca bardziej wspólność znaków niż kierunku. Oczywiście najbardziej naturalnym podejściem jest obliczenie częstości poprawności prognozy i porównanie jej z teoretycznym prawdopodobieństwem co najmniej tylu sukcesów w rozkładzie dwumianowym - tak jak pokazałem w tym artykule. Test dwumianowy wskaże czy liczba sukcesów jest istotna statystycznie, ale nie poda już faktycznej skuteczności - tj. prawdopodobieństwa po uwzględnieniu przypadkowych odchyleń.
Taką miarą jest zaś tzw. dokładność kierunkowa. rugarch oferuje specjalny test DAC (ang. Directional Accuracy Test - test kierunkowej dokładności), który ją oblicza i testuje jej istotność. Zastosowanie DAC jest proste. jego konstrukcja jest następująca:
DACTest(forecast, actual, test="PT")
Symbol "PT" oznacza test Pesaran-Timmermann. Zamiast tego można wpisać "AG" - test Anatolyev-Gerko. Oba dają bardzo zbliżone rezultaty.
Znajdujemy więc prognozy w każdym modelu. Dla ARMA(5, 3):
# prognoza w próbie
stopaZwrotuDopas <- arimaDopas$fit@fit$fitted.values
Testujemy zdolność prognostyczną modelu:
DACTest(forecast = stopaZwrotuDopas, actual = as.numeric(stopaZwrotu), test="PT")
$Test
[1] "Pesaran and Timmermann"
$Stat
[1] 5.390263
$p.value
[1] 3.51773e-08
$H0
[1] "Independently Distributed"
$Decision
[1] "Reject H0"
$DirAcc
[1] 0.6884422Prognoza zdaje test na 69% i jest istotna.
Następnie sprawdzamy to samo dla ARIMA(5, 3)-TGARCH(1, 1) oraz ARIMA(2, 2)-EGARCH(1, 1). Prognoza w próbie:
# prognoza w próbie
stopaZwrotuDopas <- fitted(dopasGarch)
Test DAC wskazał dla TGARCH dokładność kierunkową = 0.693, a dla EGARCH 0.658. Wynika z tego, że TGARCH jest najlepszy spośród wszystkich trzech modeli. Dostaliśmy dowód, że dodanie do ARMA modelu GARCH może poprawić prognozę samej części ARMA, choć mogłoby się wydawać, że są to części niezależne.
Przed utworzeniem wykresu z prognozą, warto jeszcze spojrzeć na dopasowanie reszt TGARCH do rozkładu SGED:
plot(dopasGarch, which=8)
Wygląda w miarę dobrze.
Możemy robić prognozę: dwa tygodnie do przodu, pokażemy 30 ostatnich obserwacji na wykresie:
prognoza_plus_wykres <- function(stopaZwrotu, dopasGarch, ileDoPrzodu, Ost) {
# warunkowe odchylenie std w próbie
ryzyko <- sigma(dopasGarch)
# prognoza poza próbą
prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = ileDoPrzodu)
prognozaStopaZwrotu <- fitted(prognozaGarch)
prognozaRyzyko <- sigma(prognozaGarch)
# zamiana na ts i łączenie z prognozą dla wykresu
stopaZwrotuWykres <- tail(as.numeric(stopaZwrotu), Ost)
stopaZwrotuDopasWykres <- tail(as.numeric(stopaZwrotuDopas), Ost)
ryzykoWykres <- tail(as.numeric(ryzyko), Ost)
stopaZwrotuWykres <- c(stopaZwrotuWykres, rep(NA, ileDoPrzodu))
stopaZwrotuDopasWykres <- c(stopaZwrotuDopasWykres, rep(NA, ileDoPrzodu))
prognozaStopaZwrotuWykres <- c(rep(NA, Ost-1), stopaZwrotuWykres[Ost],
as.numeric(prognozaStopaZwrotu))
# Wyłuskaj nazwę modelu wariancji
modelWariancji <- dopasGarch@model$modeldesc$vmodel
podmodelWariancji <- dopasGarch@model$modeldesc$vsubmodel
modelWariancji <- ifelse(!is.null(podmodelWariancji), podmodelWariancji, modelWariancji)
# Wyłuskaj nazwę modelu średniej
modelSredniej <- paste0("ARMA(", dopasGarch@model$modelinc["ar"], ",", dopasGarch@model$modelinc["ma"], ")")
# zmienne do wykresu
ryzykoWykres <- c(ryzykoWykres, rep(NA, ileDoPrzodu))
prognozaRyzykoWykres <- c(rep(NA, Ost-1), ryzykoWykres[Ost], as.numeric(prognozaRyzyko))
okres <- index(stopaZwrotu)
okresPrognoza <- seq.Date(from=max(okres), by="week", length.out=ileDoPrzodu+1)[-1]
okresWykres <- c(tail(okres, Ost), okresPrognoza)
# wykres
zakresX <- c(min(okresWykres), max(okresWykres))
zakresY <- c(min(na.omit(-ryzykoWykres))*3, max(na.omit(ryzykoWykres))*3)
plot(x=okresWykres, y=stopaZwrotuWykres, main="stopa zwrotu WIG20", xaxt="n", xlim=zakresX, ylim=zakresY, lwd=2, type="l", xlab="")
axis(side=1, at=okresWykres, labels=okresWykres, las=2, cex.axis=0.7)
lines(x=okresWykres, y=stopaZwrotuDopasWykres, col="lightblue", lwd=2)
lines(x=okresWykres, y=prognozaStopaZwrotuWykres, lwd=2, col="green")
lines(x=okresWykres, y=ryzykoWykres, lwd=2, col="red")
lines(x=okresWykres, y=-ryzykoWykres, lwd=2, col="red")
lines(x=okresWykres, y=prognozaRyzykoWykres, lwd=2, col="brown")
lines(x=okresWykres, y=-prognozaRyzykoWykres, lwd=2, col="brown")
abline(h=0, col="grey")
abline(v=okresWykres, col="grey")
legend(x="topleft", legend=c("stopa zwrotu",
paste("modelowana stopa zwrotu -", modelSredniej),
paste("prognoza stopy zwrotu -", modelSredniej),
paste("modelowane ryzyko -", modelWariancji),
paste("prognoza ryzyka -", modelWariancji)),
bg = "transparent", cex = 0.8,
lwd=c(2,2,2,2,2), col=c("black", "lightblue",
"green",
"red",
"brown"))
legend("bottomleft", legend=paste("korelacja stopa zwrotu / prognoza stopy zwrotu:", korelacja),
bty="o", bg="white", box.lwd=1, text.font=2)
}
prognoza_plus_wykres(stopaZwrotu, dopasGarch, ileDoPrzodu = 2, Ost = 30)
Tak więc dostaliśmy prognozę dodatnią na ten tydzień i ujemną na następny.
Pozostaje kwestia stabilności. Przypomnę, że ARMA(5,3)-TGARCH(1,1) jest niestabilny, ale sam ARMA(5, 3) był stabilny. Logiczne więc, że cała niestabilność bierze się z części modelu wariancji. To dobra wiadomość dla prognozowania stopy zwrotu. Z kolei ARMA(2, 2)-EGARCH(1, 1) jest stabilny. Wobec tego możemy użyć tego ostatniego do prognozy wariancji. Dla niego dostałem:
Zauważmy, że prognozowane odchylenia (ryzyko) lepiej się tu dopasowują do faktycznych niż przy zastosowaniu TGARCH. Jednak końcówka i prognoza odchyleń - w zasadzie są identyczne. Natomiast nie sugerowałbym się już prognozą stóp zwrotu, chociaż w tym wypadku kierunek jest ten sam.
W ogólnym przypadku zmienność, a więc i sam GARCH (bez ARMA) można wykorzystać w trejdingu na różne sposoby. Jeżeli występuje efekt Arch-In-Mean, to wzrost ryzyka będzie powiązany z oczekiwaną stopą zwrotu, co pokazałem w tym art. Jeżeli Arch-In-Mean nie występuje, tak jak tutaj, to nadal możemy wykorzystać zmienność ryzyka. W przypadku EGARCH rosnące ryzyko oznacza, że w kolejnych okresach nadal będzie ono duże, a co za tym idzie możemy spodziewać się większych zmian. W takim razie, przykładowo, czekamy na odpowiednio duży spadek, aby w tym momencie wejść na rynek. Ze względu na większą zmienność, rośnie szansa, że nastąpi odbicie. Sprzedawać można, gdy stopa zbliża się do krawędzi prognozy odchylenia. Obecnie ryzyko raczej pozostaje podwyższone. W połączeniu z dodatnią prognozą stopy spodziewałbym się większego impulsu w górę, a za tydzień w dół - oczywiście na ten moment. W przypadku spadającego ryzyka można dłużej utrzymywać stratną pozycję i czekać na zajęcie pozycji przeciwnej - gdy stopa przekroczy odchylenie - w dół - pozycja długa, w górę - pozycja krótka, mówiąc umownie.
