Prognozowanie w pakiecie rugarch jest dość skomplikowane, a słabo napisana instrukcja (bo niektórych rzeczy trzeba się domyślać) nie ułatwia sprawy. Dotychczas pokazywałem jedynie wstęp, bo model był całkowicie dopasowany do próby. Prognoza na podstawie takiego dopasowania może mieć sens jedynie dla bardzo dużych prób, ponieważ wtedy optymalne dopasowanie do błądzenia przypadkowego (czy jakiegoś podobnego procesu losowego) staje się na tyle trudne, że czysta statystyka wymusza sprowadzenie optymalnych rzędów ARMA do zera. Stąd jeżeli tych zerowych rzędów nie otrzymujemy (poprzez funkcję autoarfima), możemy przypuszczać, że jakiś nielosowy proces tutaj występuje. Aby jednak zyskać więcej pewności, model należy przetestować. Można to zrobić "na żywo", a więc sprawdzić czy nadchodzące zmiany zostaną potwierdzone przez model. To jest jednak zbyt czasochlonne i nie zawsze możliwe, dlatego musimy przeprowadzić badanie poza próbą treningową, czyli out-of-sample. Przez out-of-sample rozumiemy specjalnie wydzieloną próbę do testowania. Generalnie prognozowanie jako proces składa się z 4 etapów:
1. treningu,
2. walidacji, czyli testowania poza próbą treningową i ewentualnej poprawy modelu,
3. testu, czyli ostatniego testowania poza próbą treningową i walidacyjną, już bez poprawiania,
4. prognozy na przyszły okres.
W rugarch mamy 3 ściśle powiązane ze sobą parametry: out.sample, n.ahead oraz n.roll:
- out.sample - jest argumentem funkcji dopasowania modelu: arfimafit() i ugarchfit(). Dzięki niemu nie musimy dzielić próby na treningową i testową, bo on załatwia to za nas. Podajemy w nim liczbę, ile chcemy odjąć ostatnich obserwacji z próby.
- n.ahead - argument funkcji prognozy, czyli arfimaforecast() i ugarchforecast(). Wskazuje horyzont prognozy. Normalnie zawsze n.ahead = 1. Sprawa się komplikuje, gdy ustawione są out.sample i n.roll. Jeśli nie są one ustawione, to dla okresu powyżej 1 prognoza jest szacowana w ten sposób, że do modelu podstawiane są wartości prognoz z poprzednich okresów odpowiadające określonym rzędom modelu.
- n.roll - tak samo argument funkcji prognozy, czyli arfimaforecast() i ugarchforecast(). Wskazuje liczbę prognoz kroczących (rolling forecast), a dokładniej rzecz biorąc ile razy prognoza z n.ahead ma zostać "zrolowana", tj. powtórzona w następnych okresach. Bardziej szczegółowo: mówi ile razy wykorzystać do prognozy te obserwacje z próby testowej (out of sample), które wcześniej były prognozowane zgodnie z n.ahead. Nie może być mniejszy od out.sample, dlatego że bierze do prognozy obserwacje pochodzące zarówno z próby treningowej, jak i z próby testowej, i w ten sposób pozwala porównać prognozę z faktyczną obserwacją w tym samym okresie.
Ponieważ niniejsza analiza jest kontynuacją art. Jak wykorzystać model GARCH w trejdingu? , to wykonam ten sam przykład co ostatnio, tj. na tygodniowych stopach zwrotu WIG20. Ponieważ jednak mamy dodatkowe 2 tygodnie nowych danych (stąd zmienna Obs = 202, a nie 200), to te dwie obserwacje będą "out-of-sample", tj. przypisane do próby testowej.
# Przedwstęp
if (require("xts")==FALSE) {
install.packages("xts")
}
library("xts")
if (require("rugarch")==FALSE) {
install.packages("rugarch")
}
library("rugarch")
# Wstęp
nazwaPliku = "wig20_d.csv"
Obs <- 202
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
Poprzednio używałem funkcji autoarfima, aby znaleźć optimum. Teraz jednak nie mogę tego zrobić bez ucinania danych, bo w tej funkcji nie ma opcji out.sample. Zamiast tego musimy stworzyć obiekt arfimafit, który jest analogią do wcześniej stosowanego ugarchfit. Zatem podobnie jak w tamtej najpierw potrzebny jest model teoretyczny (specyfikacja modelu). To właśnie tutaj wpisujmy rzędy AR i MA za pomocą funkcji arfimaspec():
mojaArima <- arfimaspec(mean.model = list(armaOrder = c(rzadAR, rzadMA)), distribution.model = "sged")
Poprzednio podstawiłem rzadAR = 5, a rzadMA = 3, bo tak wskazała autoarfima. czyli
mojaArima <- arfimaspec(mean.model = list(armaOrder = c(5, 3)), distribution.model = "sged")
Następnie szukamy modelu empirycznego, czyli parametrów do teoretycznego. To właśnie tutaj ustawimy out.sample = 2:
arimaDopas <- arfimafit(spec = mojaArima, data = stopaZwrotu, out.sample = 2)
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.75Możemy sprawdzić i zobaczymy, że wszystkie wartości są identyczne jak dla tamtego modelu.
Jeżeli chodzi o prognozę tego modelu, to pokaże on poprawny kierunek na kolejny okres, mimo że modele GARCH pokazały błędny. Może to sugerować, że sama ARMA jest lepsza od ARMA-GARCH. Ustawiamy n.ahead = 1 i tworzymy wykres:
prognoza <- arfimaforecast(arimaDopas, n.ahead = 1)
prognozaSrednia <- as.numeric(fitted(prognoza))
# tworzenie dat dla wykresu
czest <- "week"
dopas <- fitted(arimaDopas)
datyWykres1 <- index(tail(stopaZwrotu, nrow(stopaZwrotu) - (nrow(dopas) - 1)))
if (length(prognozaSrednia) > length(datyWykres1) - 1) {
datyWykres2 <- seq(from=last(datyWykres1), by=czest, length.out=length(prognozaSrednia) - (length(datyWykres1) - 1) + 1)[-1]
datyWykres <- c(datyWykres1, datyWykres2)
} else {
datyWykres <- head(datyWykres1, length(prognozaSrednia) + 1)
}
# Wykres
prognoza_wykres <- xts(c(coredata(last(dopas)), prognozaSrednia), order.by=datyWykres)
plot(tail(stopaZwrotu, 30), main = "log-stopy zwrotu WIG20", extend.xaxis = TRUE)
lines(dopas, col="lightblue", lwd=3)
lines(prognoza_wykres, col="green", lwd=3)
addLegend(legend.loc = "topleft", legend.names = c("WIG20", "ARMA(5, 3)", "Prognoza"),
lwd=3, bty="o", cex=1.2, col=c("black", "lightblue", "green"))
Powiedzmy, że teraz chcemy sprawdzić prognozę na dwa okresy wprzód. Moglibyśmy więc ustawić n.ahead = 2. Okazuje się jednak, że jest to bez sensu w tym wypadku. Wydaje się, że funkcja powinna automatycznie tworzyć zrolowane prognozy, skoro już ustawiliśmy out.sample > 0. Zamiast tego podstawia same prognozy z poprzednich okresów (tzn. nie widzi prawdziwych obserwacji pochodzących z próby testowej i używa najwyżej prognozy jako zmiennej niezależnej). Niestety trzeba dodać specjalny parametr - n.roll - który będzie to robił. Trochę to głupie, ale gdy wpiszemy n.ahead = 1 i n.roll = 1, to nie dostaniemy wcale jednej prognozy naprzód jakby się wydawało, ale właśnie dwie prognozy. Popatrzmy:
prognoza <- arfimaforecast(arimaDopas, n.ahead = 1, n.roll = 1)
Tak więc n.roll bierze obserwacje z próby testowej i oblicza prognozę n.ahead = 1 naprzód. Czy ma to w ogóle uzasadnienie, żeby odróżniać n.roll od out.sample? W pewnym sensie ma, ale żeby to zobaczyć musimy zwiększyć out.sample, powiedzmy do 5.
arimaDopas <- arfimafit(spec = mojaArima, data = stopaZwrotu, out.sample = 5)
prognoza <- arfimaforecast(arimaDopas, n.ahead = 1, n.roll = 1)
Dalej to samo co wcześniej.
Model się już zmienił, więc i prognoza także. Wpiszmy teraz n.roll = 2:
prognoza <- arfimaforecast(arimaDopas, n.ahead = 1, n.roll = 2)
a wtedy dostaniemy 3 prognozy naprzód:
A więc maksymalnie n.roll może się tutaj równać 4, bo wtedy dostaniemy 5 prognoz:
Tak więc o ile na początku prognoza była prawidłowa, to potem się pogorszyła, co wskazywałoby, że model wcale nie jest tak dobry. Zauważmy jednak, że dla out.sample = 2 (oraz n.ahead = 1 i n.roll = 1), prognoza była poprawna. Być może więc model musi być przez cały czas "refitowany", czyli być modelowany dynamicznie co 1 lub 2 okresy. Problem polega na tym, że takiego modelu nie da się przetestować. Jedyną możliwością jest określenie długości minimalnej próby testowej i na jej podstawie podjąć decyzję, czy zostawić model, czy modelować na nowo. W określeniu długości tej próby pomaga sam pakiet, bo kryteria dokładności prognozy, jak MSE, MAE i DAC, można szybko uzyskać przy pomocy funkcji fpm [forecast performance measures] - jeżeli out.sample wynosi minimum 5:
fpm(prognoza)
MSE MAE DAC
1 0.0009823 0.02598 0.4DAC = 0.4 oznacza, że jedynie 2 na 5 prognoz się sprawdziło. Oznacza to, że konieczna jest poprawa modelu. Możemy przyjąć kryterium min. DAC = 0.7. Stosujemy więc autoarfimę, ale bez ostatnich 5 obserwacji:
klaster <- makeCluster(detectCores() - 1)
arimaDopas <- autoarfima(head(stopaZwrotu, -5), ar.max = 10, ma.max = 10, criterion = "HQIC",
distribution.model = "sged", cluster = klaster,
method="partial")
stopCluster(klaster)
Otrzymałem ARMA(2, 2). Ale okazało się, że model jest jeszcze gorszy (DAC = 0.2):
Są dwie rzeczy, które możemy zmienić: rozkład reszt albo kryterium dokładności. Zacznijmy od rozkładu. Oczywiście najlepiej by było sprawdzić w ogóle dopasowanie rozkładu empirycznego do teoretycznego dla każdego przypadku, ale na razie jest to mniej istotne. Kod sprawdzający każdy rozkład:
klaster <- makeCluster(detectCores() - 1)
wyrzuc <- 5
rozklady = c("norm", "snorm", "std", "sstd", "ged", "sged", "nig", "ghyp", "jsu")
wyniki <- data.frame()
clusterExport(klaster, varlist = c("stopaZwrotu", "wyrzuc", "rozklady"))
for (rozklad in rozklady) {
clusterExport(klaster, varlist = "rozklad")
tryCatch({
arimaDopas <- autoarfima(head(stopaZwrotu, -wyrzuc), ar.max = 10, ma.max = 10, criterion = "HQIC",
distribution.model = rozklad, cluster = klaster,
method="partial")
rzadAR <- arimaDopas$rank.matrix[1,1]
rzadMA <- arimaDopas$rank.matrix[1,2]
arimaDopas <- arfimafit(spec = arfimaspec(mean.model = list(armaOrder = c(rzadAR, rzadMA)),
distribution.model = rozklad), data = stopaZwrotu, out.sample = wyrzuc)
# Prognoza
prognoza <- arfimaforecast(arimaDopas, n.ahead = 1, n.roll = wyrzuc-1)
wynikiNowe <- data.frame(rzad_AR = rzadAR, rzad_MA = rzadMA,
rozklad_reszt = rozklad, DAC = fpm(prognoza)$DAC)
wyniki <- rbind(wyniki, wynikiNowe)
}, error = function(e) {
message(paste("Błąd: ", e$message))
})
}
stopCluster(klaster)
wyniki
rzad_AR rzad_MA rozklad_reszt DAC
1 6 7 norm 0.6
2 6 6 snorm 0.4
3 5 4 std 0.4
4 6 10 sstd 0.4
5 2 2 sged 0.2
6 6 8 nig 0.4
7 6 7 ghyp 0.4
8 5 5 jsu 0.4Okazało się, że rozkład normalny przyniósł największy DAC = 0.6, co może dziwić, ale zauważmy, że używamy tu logarytmicznych stóp zwrotu, które "normalizują" duże wahania. W każdym razie 0.6 to ciągle za mało.
W takim razie musimy zmienić kryterium dopasowania. HQIC stanowi kompromis między AIC a BIC, dlatego wydaje się najlepszy dla modelu z autoregresją. Raczej uznaje się, że AIC jest lepszy do prognozowania niż BIC [zob. Cavanaugh, J. E., & Neath, A. A. (2019). "The Akaike Information Criterion: Background, derivation, properties, application, interpretation, and refinements". WIREs computational statistics, 11]. Problem polega na długości próby - żeby AIC miało sens, próba musi być znacznie dłuższa niż 200 obserwacji. HQIC jest tu dobrym rozwiązaniem. Natomiast oczywistą pokusą jest bezpośrednie zastosowanie DAC jako kryterium wyboru modelu. Problem z nim jest taki sam jak z MAE i MSE, które nie chronią przed przeuczeniem (overfitting), co robią AIC, BIC i HQC [ibidem]. Gdy będziemy dodawać kolejne zmienne, MAE i MSE będą się zmniejszać, dając pozory poprawiania się modelu (im więcej zmiennych, tym większa możliwość dopasowania do dowolnej próby). Dopiero na próbie testowej widać, że to złudzenie. Z drugiej strony, jeśli interesuje nas bardziej kierunek niż wartość, DAC może stanowić dodatkowy warunek obok kryterium informacyjnego.
Wprowadzimy więc zasadę ad hoc, że robimy dwa testy DAC - na próbie treningowej (na podstawie modelu uzyskanego z HQC) oraz na próbie testowej. To znaczy, otrzymaliśmy teraz ARMA(6, 7), więc sprawdzamy DAC traktując dopasowanie modelu jako prognozę, którą porównamy z próbą treningową:
wynikiSort <- wyniki[order(wyniki$DAC, decreasing = TRUE), ]
rzadAR <- wynikiSort$rzad_AR[1]
rzadMA <- wynikiSort$rzad_MA[1]
rozklad <- wynikiSort$rozklad_reszt[1]
# Sprawdzamy DAC dla optymalnego modelu:
mojaArma <- arfimaspec(mean.model = list(armaOrder = c(rzadAR, rzadMA)),
distribution.model = rozklad)
arimaDopas <- arfimafit(spec = mojaArma, data = stopaZwrotu,
out.sample = wyrzuc, solver = "hybrid")
dopas <- fitted(arimaDopas)
# część wspólna danych i modelu
polaczone <- merge(stopaZwrotu, dopas, join="inner")
stopaZwrotuZgodna <- polaczone[, 1]
dopasZgodna <- polaczone[, 2]
DACTest(forecast = as.numeric(dopasZgodna),
actual = as.numeric(stopaZwrotuZgodna),
test = "PT")
$Test
[1] "Pesaran and Timmermann"
$Stat
[1] 7.091
$p.value
[1] 0.0000000000006655
$H0
[1] "Independently Distributed"
$Decision
[1] "Reject H0"
$DirAcc
[1] 0.7551Czyli prawidłowy kierunek w próbie treningowej wyniósł 75,5%. Przypomnę, że w próbie testowej już tylko 0.6, chociaż jej długość była trochę za mała, by tak porównywać. Dlatego zwiększmy ją do 10. Wyznaczmy jeszcze raz dopasowanie dla ARMA(6, 7), ale z out.sample = 10:
mojaArma <- arfimaspec(mean.model = list(armaOrder = c(6, 7)),
distribution.model = "norm")
arimaDopas <- arfimafit(spec = mojaArma, data = stopaZwrotu,
out.sample = 10)
Wówczas ustawiamy n.ahead = 1 i n.roll = 9:
prognoza <- arfimaforecast(arimaDopas, n.ahead = 1, n.roll = wyrzuc-1)
I sprawdzamy dokładność prognozy:
fpm(prognoza)
MSE MAE DAC
1 0.001483 0.03301 0.6DAC nadal wynosi 0.6, co pokazuje, że DAC z próby testowej będzie faktycznie gorszy niż w treningowej. Model jest za słaby, by go wykorzystać, zresztą popatrzmy jak to wygląda na wykresie:
Pomijając diagnostykę, która powinna parę rzeczy podpowiedzieć, można sprawdzić, czy dodanie GARCH wpłynie pozytywnie na model. Dodajmy EGARCH(1, 1), który na tygodniowych zmianach wydaje się być niezły (zob. tu). Rozkładu nie zmieniamy (zostajemy przy normalnym).
# Dodanie GARCH
mojGarch <- ugarchspec(variance.model = list(model = "eGARCH", submodel = NULL,
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(6, 7), archm = FALSE),
distribution.model = "norm")
I tak samo zachowujemy out.sample = 10, który dodajemy do ugarchfit():
dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu, out.sample = 10)
dopasGarch
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : eGARCH(1,1)
Mean Model : ARFIMA(6,0,7)
Distribution : norm
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.001214 0.000007 171.982 0
ar1 -0.786039 0.000026 -30461.837 0
ar2 -0.541886 0.000042 -12762.906 0
ar3 0.647248 0.000050 13058.421 0
ar4 -0.545099 0.000044 -12263.428 0
ar5 -0.725927 0.000056 -13005.354 0
ar6 -0.850517 0.000063 -13557.083 0
ma1 0.811287 0.000050 16196.014 0
ma2 0.621477 0.000051 12141.818 0
ma3 -0.632186 0.000026 -24094.748 0
ma4 0.641535 0.000050 12857.551 0
ma5 0.828067 0.000056 14806.597 0
ma6 0.960052 0.000056 17113.362 0
ma7 0.000509 0.000015 34.761 0
omega -0.420318 0.000703 -597.585 0
alpha1 -0.213766 0.000329 -649.429 0
beta1 0.942937 0.000133 7099.713 0
gamma1 -0.078948 0.000168 -470.651 0
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.001214 0.000008 149.87 0
ar1 -0.786039 0.000116 -6759.62 0
ar2 -0.541886 0.000107 -5054.70 0
ar3 0.647248 0.000155 4182.40 0
ar4 -0.545099 0.000152 -3592.58 0
ar5 -0.725927 0.000152 -4784.95 0
ar6 -0.850517 0.000197 -4326.03 0
ma1 0.811287 0.000230 3523.51 0
ma2 0.621477 0.000202 3082.78 0
ma3 -0.632186 0.000057 -10997.27 0
ma4 0.641535 0.000174 3677.13 0
ma5 0.828067 0.000252 3292.35 0
ma6 0.960052 0.000263 3646.23 0
ma7 0.000509 0.000025 20.16 0
omega -0.420318 0.000905 -464.42 0
alpha1 -0.213766 0.001755 -121.78 0
beta1 0.942937 0.000788 1196.43 0
gamma1 -0.078948 0.000449 -175.70 0
LogLikelihood : 436.7
Information Criteria
------------------------------------
Akaike -4.3846
Bayes -4.0781
Shibata -4.4004
Hannan-Quinn -4.2604
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 4.033 0.04463
Lag[2*(p+q)+(p+q)-1][38] 24.887 0.00000
Lag[4*(p+q)+(p+q)-1][64] 39.001 0.06727
d.o.f=13
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.5674 0.4513
Lag[2*(p+q)+(p+q)-1][5] 2.9719 0.4122
Lag[4*(p+q)+(p+q)-1][9] 5.7146 0.3323
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 1.056 0.500 2.000 0.3042
ARCH Lag[5] 4.408 1.440 1.667 0.1405
ARCH Lag[7] 5.962 2.315 1.543 0.1438
Nyblom stability test
------------------------------------
Joint Statistic: 3.772
Individual Statistics:
mu 0.01613
ar1 0.03226
ar2 0.04187
ar3 0.01256
ar4 0.02268
ar5 0.04005
ar6 0.02935
ma1 0.04190
ma2 0.02966
ma3 0.01215
ma4 0.03250
ma5 0.04756
ma6 0.01888
ma7 0.01562
omega 0.17387
alpha1 0.32652
beta1 0.16804
gamma1 0.29282
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 3.83 4.14 4.73
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.3620 0.7177
Negative Sign Bias 0.6644 0.5072
Positive Sign Bias 1.0742 0.2841
Joint Effect 1.6372 0.6510
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 11.62 0.9013
2 30 14.60 0.9879
3 40 31.20 0.8087
4 50 39.63 0.8280
Ewidentnie model nie jest doskonały, bo występują autokorelacje w resztach (test Ljung-Boxa), ale za to jest stabilny. DAC dla próby treningowej wyniósł ok. 0.7.
W funkcji ugarchforecast() ustawiamy n.ahead = 1, n.roll = 9:
prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = 1, n.roll = 9)
fpm(prognozaGarch)
Okazuje się, że DAC dla próby testowej spadł do 0.5, a więc GARCH nie pomógł. Niemniej model może mieć wartość z innej perspektywy, co zobaczymy na wykresie.
Cała reszta to już techniczne dostosowanie w celu utworzenia wykresu:
dopas <- fitted(dopasGarch)
ryzyko <- sigma(dopasGarch)
prognozaSrednia <- as.numeric(fitted(prognozaGarch))
prognozaRyzyko <- as.numeric(sigma(prognozaGarch))
# tworzenie dat dla wykresu
czest <- "week"
datyWykres1 <- index(tail(stopaZwrotu, nrow(stopaZwrotu) - (nrow(dopas) - 1)))
if (length(prognozaSrednia) > length(datyWykres1) - 1) {
datyWykres2 <- seq(from=last(datyWykres1), by=czest, length.out=length(prognozaSrednia) - (length(datyWykres1) - 1) + 1)[-1]
datyWykres <- c(datyWykres1, datyWykres2)
} else {
datyWykres <- head(datyWykres1, length(prognozaSrednia) + 1)
}
# Wykres
prognoza_wykres <- xts(c(coredata(last(dopas)), prognozaSrednia), order.by=datyWykres)
prognozaRyzyko_wykres <- xts(c(coredata(last(ryzyko)), prognozaRyzyko), order.by=datyWykres)
plot(tail(stopaZwrotu, 30), main = "log-stopy zwrotu WIG20", extend.xaxis = TRUE)
lines(dopas, col="lightblue", lwd=3)
lines(prognoza_wykres, col="green", lwd=3)
#addLegend(legend.loc = "topleft", legend.names = c("WIG20", paste("ARMA(", rzadAR,",", rzadMA,")"), "Prognoza"),
lines(tail(ryzyko, 30), col="red", lwd=2)
lines(tail(-ryzyko, 30), col="red", lwd=2)
lines(prognozaRyzyko_wykres, col="darkred", lwd=2)
lines(-prognozaRyzyko_wykres, col="darkred", lwd=2)
addLegend(legend.loc = "topleft", legend.names = c("WIG20", paste("ARMA(", 6,",", 7,")"), "Prognoza", "Ryzyko - EGARCH(1,1)", "Prognoza ryzyka"),
lwd=3, bty="o", cex=1, col=c("black", "lightblue", "green", "red", "darkred"))
Jeżeli chcemy to połączyć z prognozą na okres poza próbą testową, wystarczy zwiększyć n.roll o 1, czyli wpisujemy n.ahead = 1 oraz n.roll = 10:
prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = 1, n.roll = 10)
Dzięki warunkowej funkcji tworzącej daty oraz dzięki argumentowi extend.xaxis = TRUE wewnątrz funkcji plot() z pakietu xts pozostały kod się nie zmienia, więc nie będę go powtarzał. Wykres:
