wtorek, 1 października 2024

Jak używać argumentów n.ahead, out.sample i n.roll w rugarch (język R)?

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.75


Moż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.4


DAC = 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.4


Okazał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.7551


Czyli 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.6


DAC 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:


 
Ta prognoza sugeruje spadek ryzyka w tym tygodniu, czyli mniejsze odchylenia. Dodatkowo wskazuje spadek WIG20, chociaż trzeba pamiętać, że ten model nie ma mocy prognostycznej w sensie dodatniej lub ujemnej stopy zwrotu. Natomiast model stworzył parotygodniową cykliczność na stopie zwrotu. Jeśli rozszerzymy n.ahead, to jaśniej to zobaczymy - najlepiej w ogóle ustawić n.roll = 0. Np. ustawienie:

prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = 40, n.roll = 0)

daje taki obraz:



Zauważmy, że pomimo braku rolowania prognozy, prawie się ona nie zmieniła w próbie testowej.  Oczywiście ta cykliczność może być także złudzeniem i nie należy do tego modelu się przywiązywać. Jest to raczej przykład, punkt wyjścia do dalszych analiz. Najważniejsze było tu pokazać, jak testować model i weryfikować jego zdolność prognostyczną w rugarch przy pomocy trzech tytułowych parametrów.

Brak komentarzy:

Prześlij komentarz