czwartek, 17 października 2024

Dwa rodzaje prognoz kroczących w rugarch (język R)

Każdy model prognostyczny może być rozumiany na dwa sposoby. Pierwszy wynika z założenia, że istnieje "prawdziwy model", którym można prognozować pewne wielkości w kolejnych okresach. Nie chodzi tu o prawdziwy model w ścisłym znaczeniu tego słowa, ale o model wystarczająco dobry w pewnym oknie czasowym. Z tego założenia wywodzi się sama możliwość użycia próby testowej (out-of-sample) w ten sposób, że do modelu podstawiamy albo prognozy z przeszłości (argument n.ahead w rugarch omówiony w poprzednim artykule) albo faktyczne wartości, których model "nie widział" podczas estymacji parametrów (n.roll dla prognozy kroczącej). W przypadku prognozy kroczącej (rolling forecast) przedłużamy niejako model w przyszłość, używając możliwie najbardziej aktualnych danych, ale przy zachowaniu starych parametrów, i to się czasem nazywa rolowaniem prognozy. 

Drugi sposób modelowania wynika z założenia, że nie istnieje prawdziwy model i to czym dysponujemy należy przez cały czas reestymować. W świecie ekonomii i giełdy jest to naturalne założenie, a nawet stanowi coś w rodzaju aksjomatu. Problem z nim jest taki, że dużo trudniej jest testować taki model. Próba testowa jest równa treningowej, bo każda kolejna obserwacja jest dodawana do próby uczącej, która z kolei wykorzystywana jest do nowej estymacji parametrów. 

W praktyce obydwa sposoby mieszane są ze sobą lub używane równorzędnie. I tak najczęściej prognoza krocząca zaczyna się dopiero od jakiegoś dalszego okresu, a reestymacja nie odbywa się co okres, ale np. co kilka okresów, co przyspiesza proces obliczeń. Równoległe zaś użycie obu metod umożliwia ocenę pogarszania się modelu w czasie. Jeżeli bowiem zakładamy, że parametry modelu są zmienne w czasie - są dynamiczne, to jeśli nasza prognoza też jest dynamiczna, to powinna być lepsza niż statyczna (czyli niż bez reestymacji). Jeżeli zatem widzimy, że prognoza dynamiczna jest gorsza od statycznej, to dostajemy dowód, że model jest błędny - przypadkowo się dostroił do próby testowej. 

Inaczej mówiąc, o ile użycie kryterium informacyjnego (AIC, BIC, HQC) może dać złudzenie uzyskania dobrego modelu z powodu sztucznego dopasowania go do próby uczącej, o tyle prognoza krocząca bez reestymacji może dać złudzenie uzyskania dobrego modelu z powodu sztucznego dopasowania do próby testowej. I tak jak prognoza krocząca bez reestymacji jest sposobem, żeby odkryć pierwsze złudzenie, o tyle prognoza z reestymacją jest sposobem, by odkryć drugie złudzenie. 

Prognozę bez reestymacji omówiłem wcześniej. Teraz zajmę się drugim przypadkiem. Prognozę kroczącą z reestymacją można oczywiście samodzielnie zakodować, ale rugarch oferuje specjalną funkcję: arfimaroll() i ugarchroll(). Ta druga jest bardziej ogólna, dlatego to nią się posłużymy.

Funkcja ma następujące argumenty:

spec - teoretyczny model GARCH,

data - dane najlepiej w xts,

forecast.length - długość prognozy dla próby testowej (tj. out-of-sample), która liczona jest końca całej próby

n.start - można użyć zamiast forecast.length, zamiast liczyć od końca można po prostu określić punkt startowy dla prognozy

refit.every - określa, co którą obserwację w próbie testowej ponownie estymować; jeśli więc chcemy aby dla każdej następnej obserwacji reestymować, wpisujemy 1,

refit.window - wybór czy reestymować od początku próby do ostatniej dostępnej obserwacji ("recursive"), czy też reestymować w oknie o stałej liczebności, czyli tak jak np. średnia 30-okresowa ("moving"),

window.size - jeżeli wybraliśmy refit.window = "moving", wtedy ustawiamy tu długość okna, w którym estymowane są parametry; w przeciwnym wypadku dajemy NULL,

solver - zazwyczaj wybieram "hybrid", bo jest najbardziej wszechstronny

cluster - również tu stosujemy przetwarzanie równoległe, aby przyspieszyć obliczenia (na podstawie pakietu parallel)


Mimo dynamiczności kroczącej prognozy, trzeba mieć na względzie, że model teoretyczny w sensie jego rodzaju nie zmienia postaci podczas rolowania prognozy. Szukamy więc najlepszego modelu. Dotychczas posługiwałem się funkcją auto.arfima(). Ponieważ jednak zdecydowałem posłużyć się ogólniejszym modelem, tj. ARMA-GARCH(1, 1), musiałem zbudować funkcję, która będzie iteracyjnie obliczać statystyki dla każdego rzędu ARMA. jak widać rzędy GARCH pozostawiam na poziomie (1, 1). Sprawdzanie kolejnych rzędów do samego GARCH bardzo wydłużyłoby analizę. Z moich dotychczasowych prób wynikało, że GARCH(1, 1) jest całkowicie wystarczający, ale kiedyś można będzie to lepiej ocenić. 

Funkcja będzie poszukiwać najlepszego modelu wg trzech kryteriów: najniższego HQIC, najwyższego DAC w próbie uczącej i najwyższego DAC w próbie testowej. Aby cały proces przyspieszyć trzeba użyć przetwarzania równoległego. O ile rugarch importuje sam pakiet parallel, ja użyję innego - doParallel. Przydatny będzie też w tym celu pakiet foreach. Poza tym zastosuję funkcję bind_rows() z pakietu dplyr. Pakiet ten jest częścią większego pakietu tidyverse. Reszta jak wcześniej (czyli xts i rugarch). Cała analiza jest kontynuacją poprzednich (1, 2, 3, 4). W sumie zaczynamy od:

# Przedwstęp

if (require("tidyverse")==FALSE) {

  install.packages("tidyverse")

}

library("dplyr")

if (require("xts")==FALSE) {

  install.packages("xts")

}

library("xts")

if (require("rugarch")==FALSE) {

  install.packages("rugarch")

}

library("rugarch")

if (require("foreach")==FALSE) {

  install.packages("foreach")

}

library("foreach")

if (require("doParallel")==FALSE) {

  install.packages("doParallel")

}

library("doParallel")

#precyzja i sposób zapisu liczb: 

options(scipen = 100, digits = 4)


Kontynuuję analizę tygodniowych (logarytmicznych) stóp zwrotu WIG20, ale tym razem 230 ostatnich danych - tak żeby mieć 200 poza próbą testową (źródło: stooq.pl). Ostatnie notowanie 11.10.2024:

# Wstęp

nazwaPliku = "wig20_d.csv"

Obs <- 231

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


 Używam modelu EGARCH(1, 1) o rozkładzie skośno-normalnym. 

Może najpierw pokażę prosty i krótki kod na szukanie optimum, bez obliczeń równoległych:

# stare

rozklad <- "snorm"

hqc <- list()

dacDopas <- list()

dacPrognoza <- list()

model_Garch <- "eGARCH"

submodel_Garch <- NULL

wyrzuc <- 30

for (k in 0:1) {

  for (i in 0:15) {

    for (j in 0:15) {

      tryCatch({

        dopasGarch <- NULL

        rzad_modelu <- sprintf("ARMA(%d, %d)%s", i, j, ifelse(k == 1, "-ARCH_In_Mean", ""))

        mojGarch <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,

                                                     garchOrder = c(1, 1)), 

                               mean.model = list(armaOrder = c(i, j), archm = as.logical(k)), 

                               distribution.model = rozklad)

        dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu, out.sample=wyrzuc, 

                                #solver = "nlminb")

                                solver = "hybrid")

        # HQC

        hqc[[rzad_modelu]] <- infocriteria(dopasGarch)[4,]

        # DAC in sample

        stopaZwrotuDopas <- fitted(dopasGarch)

        stopaZwrotuZgodna <- stopaZwrotuDopas + residuals(dopasGarch)

        dacDopas[[rzad_modelu]] <- DACTest(forecast = as.numeric(stopaZwrotuDopas), 

                                actual = as.numeric(stopaZwrotuZgodna), 

                                test = "PT")$DirAcc

        # wg DAC out of sample

        prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = 1, n.roll = wyrzuc-1)

        dacPrognoza[[rzad_modelu]] <- fpm(prognozaGarch)[3]

        names(dacPrognoza[[rzad_modelu]]) <- rzad_modelu

      }, error = function(e) {

        message("Błąd w ARMA(", i, ", ", j, ") z k = ", k, ": ", conditionMessage(e))

        return(list(error = conditionMessage(e), rzad_modelu = rzad_modelu))

      })

    }

  }

}


Gdy jednak ustawiamy max AR i MA na więcej niż 10, czas mocno się wydłuża, szczególnie dla rozkładu niegaussowskiego. Zoptymalizowałem więc kod, tak by używał przetwarzania równoległego, jak już było wspomniane. Powyższy kod może wyglądać wtedy następująco:

# to samo co wyżej z parallel

rozklad <- "snorm"

hqc <- list()

dacDopas <- list()

dacPrognoza <- list()

model_Garch <- "eGARCH"

submodel_Garch <- NULL

wyrzuc <- 30

wynik <- list()


# Ustawienie równoległego backendu do użycia wielu procesorów

ileRdzeni <- detectCores() - 1  # Użyj jednego mniej niż całkowita liczba rdzeni

klaster <- makeCluster(ileRdzeni)

registerDoParallel(klaster)


# Utworzenie siatki parametrów dla pętli

siatka_parametrow <- expand.grid(k = 0:1, i = 0:15, j = 0:15)


# Użycie foreach do przetwarzania równoległego

wyniki <- foreach(param = 1:nrow(siatka_parametrow), .packages = c("rugarch")) %dopar% {

  k <- siatka_parametrow$k[param]

  i <- siatka_parametrow$i[param]

  j <- siatka_parametrow$j[param]

  

  rzad_modelu <- sprintf("ARMA(%d, %d)%s", i, j, ifelse(k == 1, "-ARCH_In_Mean", ""))

  

  tryCatch({

    dopasGarch <- NULL

    

    mojGarch <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,

                                                 garchOrder = c(1, 1)), 

                           mean.model = list(armaOrder = c(i, j), archm = as.logical(k)), 

                           distribution.model = rozklad)

    dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu, out.sample=wyrzuc,

                            solver = "hybrid")

    

    stopaZwrotuDopas <- fitted(dopasGarch)

    stopaZwrotuZgodna <- stopaZwrotuDopas + residuals(dopasGarch)

    

    prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = 1, n.roll = wyrzuc - 1)

    

    # Inicjalizacja wartości DAC

    dacDopas_wartosc <- NA

    dacPrognoza_wartosc <- NA

    

    # Sprawdzenie ARMA(0, 0) i pominięcie DACTest jeśli prawda

    if (!(i == 0 && j == 0 && k == 0)) {

      dacDopas_wartosc <- DACTest(forecast = as.numeric(stopaZwrotuDopas), 

                                  actual = as.numeric(stopaZwrotuZgodna), 

                                  test = "PT")$DirAcc

      dacPrognoza_wartosc <- as.numeric(fpm(prognozaGarch)[3]) # Liczenie DAC dla prognozy

    }

    

    return(list(

      rzad_modelu = rzad_modelu,

      hqc = infocriteria(dopasGarch)[4,],

      dacDopas = dacDopas_wartosc,

      dacPrognoza = dacPrognoza_wartosc

    ))

    

  }, error = function(e) {

    message("Błąd w ARMA(", i, ", ", j, ") z k = ", k, ": ", conditionMessage(e))

    return(list(error = conditionMessage(e), rzad_modelu = rzad_modelu))

  })

  

  return(wynik)

}


# Zatrzymaj klaster po przetwarzaniu

stopCluster(klaster)


# Przetwórz wyniki i przechowaj je w oryginalnych listach

for (wynik in wyniki) {

  if (!is.null(wynik)) {

    for (nazwa_modelu in names(wynik)) {

      hqc[[wynik$rzad_modelu]] <- wynik$hqc

      dacDopas[[wynik$rzad_modelu]] <- wynik$dacDopas

      dacPrognoza[[wynik$rzad_modelu]] <- wynik$dacPrognoza

    }

  }

}


Tak więc staje się dużo bardziej złożony. Dalsza część to już tylko odpowiednie wyciągnięcie danych:

# kryterium inf HQC
hqcDf <- as.data.frame(do.call(cbind, hqc))
colnames(hqcDf) <- names(hqc)
rownames(hqcDf) <- "HQC" 

# DAC w próbie treningowej (uczącej)
dacDopasDf <- as.data.frame(do.call(cbind, dacDopas))
rownames(dacDopasDf) <- "DAC (in sample)" 

# DAC w próbie testowej
dacPrognozaDf <- as.data.frame(do.call(cbind, dacPrognoza))
rownames(dacPrognozaDf) <- "DAC (out of sample)" 

porownanie <- bind_rows(hqcDf, dacDopasDf, dacPrognozaDf)

kryterium_hqc <- porownanie[, which.min(porownanie[1, ]), drop=F]
kryterium_dacDopas <- porownanie[, which.max(porownanie[2, ]), drop=F]
kryterium_dacPrognoza <- porownanie[, which.max(porownanie[3, ]), drop=F]

# Znajdź wszystkie kolumny z minimalną wartością w pierwszym wierszu (HQC)
min_hqc <- min(round(porownanie[1, ], 4), na.rm = TRUE)
kryterium_hqc <- porownanie[, which(round(porownanie[1, ], 4) == min_hqc), drop = FALSE]

# Znajdź wszystkie kolumny z maksymalną wartością w drugim wierszu (DAC w próbie treningowej)
max_dacDopas <- max(round(porownanie[2, ], 4), na.rm = TRUE)
kryterium_dacDopas <- porownanie[, which(round(porownanie[2, ], 4) == max_dacDopas), drop = FALSE]

# Znajdź wszystkie kolumny z maksymalną wartością w trzecim wierszu (DAC w próbie testowej)
max_dacPrognoza <- max(round(porownanie[3, ], 4), na.rm = TRUE)
kryterium_dacPrognoza <- porownanie[, which(round(porownanie[3, ], 4) == max_dacPrognoza), drop = FALSE]

kryterium_hqc
kryterium_dacDopas
kryterium_dacPrognoza

Efekt:


Najważniejsze jest pierwsze i ostatnie. Natomiast DAC w próbie uczącej (in sample) służy do porównania z DAC w próbie testowej (out of sample)  - czy się pogorszyło, czy poprawiło i o ile.

Rodzą się dwa spostrzeżenia. Po pierwsze wszystkie 3 kryteria dały podobne optima, co jest dobrym znakiem. W przypadku HQC i DAC dla próby uczącej dostajemy to samo optimum, ARMA(13, 10). Szybko jednak zauważymy, że ten model daje zupełnie przypadkowe wyniki z punktu widzenia DAC w próbie testowej. Nie znaczy to jeszcze, że model jest zupełnie bezużyteczny, bo znak stopy zwrotu jest jednym z prognozowanych elementów, a ponadto nie uwzględniliśmy w ogóle reestymacji. Za to najlepszy model dla DAC na próbie testowej to ARMA(12, 8)-Arch_In_Mean. Na 30 obserwacji ok. 23 jest poprawnie prognozowanych (76,67%).  Dalej wyciągamy parametry:

# parametry - HQC

wyciagnij_liczby <- function(wyr, x) {

  as.numeric(gsub("[^0-9]", "", substr(wyr, x, x+2), ))

}


rzadAR_hqc <- wyciagnij_liczby(colnames(kryterium_hqc), 6)

rzadMA_hqc <- wyciagnij_liczby(colnames(kryterium_hqc), 9)

czy_archm_hqc <- grepl("ARCH_In_Mean", colnames(kryterium_hqc))


# parametry - DAC out of sample


wyciagnij_liczby <- function(wyr, x) {

  as.numeric(gsub("[^0-9]", "", substr(wyr, x, x+2), ))

}


rzadAR_DAC <- wyciagnij_liczby(colnames(kryterium_dacPrognoza), 6)[1]

rzadMA_DAC <- wyciagnij_liczby(colnames(kryterium_dacPrognoza), 9)[1]

czy_archm_DAC <- grepl("ARCH_In_Mean", colnames(kryterium_dacPrognoza))[1]


Budujemy modele wg tych dwóch kryteriów.

1. wg HQIC

# budowa modelu wg HQC
mojGarch_hqc <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,
                                             garchOrder = c(1, 1)), 
                         mean.model = list(armaOrder = c(rzadAR_hqc, rzadMA_hqc), archm = czy_archm_hqc), 
                       distribution.model = rozklad)


dopasGarch_hqc <- ugarchfit(spec = mojGarch_hqc, data = stopaZwrotu, out.sample = wyrzuc, solver="hybrid")
dopasGarch_hqc    

> dopasGarch_hqc

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: eGARCH(1,1)
Mean Model	: ARFIMA(13,0,10)
Distribution	: snorm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.001778    0.000001    3187.4        0
ar1     0.613189    0.000060   10174.9        0
ar2    -0.230097    0.000038   -6082.9        0
ar3     0.222988    0.000043    5243.7        0
ar4     0.217702    0.000065    3368.1        0
ar5    -0.327908    0.000031  -10717.8        0
ar6    -0.045269    0.000000 -106355.2        0
ar7     0.204282    0.000025    8227.4        0
ar8    -0.295914    0.000029  -10116.2        0
ar9     0.914463    0.000003  327459.8        0
ar10   -0.282518    0.000051   -5548.4        0
ar11    0.049116    0.000005   10504.4        0
ar12   -0.271441    0.000008  -32509.3        0
ar13    0.137039    0.000020    6840.4        0
ma1    -0.231903    0.000043   -5375.9        0
ma2     0.200755    0.000037    5378.9        0
ma3    -0.006765    0.000002   -4222.9        0
ma4    -0.347592    0.000053   -6613.8        0
ma5     0.336234    0.000058    5793.3        0
ma6     0.136536    0.000055    2460.6        0
ma7    -0.346302    0.000069   -4984.0        0
ma8     0.448850    0.000068    6615.9        0
ma9    -1.371759    0.000196   -7006.0        0
ma10   -0.033945    0.000011   -3077.0        0
omega  -0.947041    0.000183   -5171.6        0
alpha1 -0.386418    0.000083   -4648.3        0
beta1   0.875887    0.000131    6680.6        0
gamma1 -0.303940    0.000138   -2201.8        0
skew    0.721121    0.000567    1270.8        0

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.001778    0.000003    572.79        0
ar1     0.613189    0.000158   3888.48        0
ar2    -0.230097    0.000052  -4414.02        0
ar3     0.222988    0.000133   1678.64        0
ar4     0.217702    0.000236    923.73        0
ar5    -0.327908    0.000005 -62812.59        0
ar6    -0.045269    0.000011  -4310.29        0
ar7     0.204282    0.000040   5126.05        0
ar8    -0.295914    0.000116  -2561.42        0
ar9     0.914463    0.000487   1878.31        0
ar10   -0.282518    0.000185  -1528.32        0
ar11    0.049116    0.000007   6918.33        0
ar12   -0.271441    0.000099  -2749.41        0
ar13    0.137039    0.000032   4308.67        0
ma1    -0.231903    0.000103  -2252.38        0
ma2     0.200755    0.000146   1378.90        0
ma3    -0.006765    0.000004  -1662.03        0
ma4    -0.347592    0.000118  -2938.63        0
ma5     0.336234    0.000067   5008.16        0
ma6     0.136536    0.000331    412.02        0
ma7    -0.346302    0.000170  -2038.25        0
ma8     0.448850    0.000028  15997.18        0
ma9    -1.371759    0.000816  -1681.73        0
ma10   -0.033945    0.000046   -735.19        0
omega  -0.947041    0.001079   -877.72        0
alpha1 -0.386418    0.000124  -3120.09        0
beta1   0.875887    0.000227   3858.24        0
gamma1 -0.303940    0.000381   -798.69        0
skew    0.721121    0.000122   5933.58        0

LogLikelihood : 486 

Information Criteria
------------------------------------
                    
Akaike       -4.5698
Bayes        -4.0915
Shibata      -4.6051
Hannan-Quinn -4.3762

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                          statistic p-value
Lag[1]                        2.234  0.1350
Lag[2*(p+q)+(p+q)-1][68]     22.335  1.0000
Lag[4*(p+q)+(p+q)-1][114]    38.545  0.9999
d.o.f=23
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      3.788 0.05163
Lag[2*(p+q)+(p+q)-1][5]     4.638 0.18454
Lag[4*(p+q)+(p+q)-1][9]     5.538 0.35487
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]   0.01662 0.500 2.000  0.8974
ARCH Lag[5]   1.46976 1.440 1.667  0.6003
ARCH Lag[7]   1.97098 2.315 1.543  0.7235

Nyblom stability test
------------------------------------
Joint Statistic:  no.parameters>20 (not available)
Individual Statistics:               
mu     0.026996
ar1    0.011957
ar2    0.012313
ar3    0.011328
ar4    0.013683
ar5    0.015958
ar6    0.013125
ar7    0.007785
ar8    0.008588
ar9    0.007616
ar10   0.006301
ar11   0.006253
ar12   0.008755
ar13   0.009341
ma1    0.013380
ma2    0.007712
ma3    0.016855
ma4    0.013099
ma5    0.007400
ma6    0.009297
ma7    0.013352
ma8    0.013728
ma9    0.013263
ma10   0.007913
omega  0.386518
alpha1 0.245545
beta1  0.405043
gamma1 0.825370
skew   0.776156

Asymptotic Critical Values (10% 5% 1%)
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias            1.027 0.3056    
Negative Sign Bias   1.329 0.1853    
Positive Sign Bias   1.295 0.1968    
Joint Effect         5.945 0.1143    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20      24.6      0.17413
2    30      33.7      0.25041
3    40      56.8      0.03258
4    50      63.0      0.08625


Parametry skośności i odpowiedzialne za efekt EGARCH są trochę niestabilne (wartość powyżej 0,47), ale poza tym jest ok - nie ma autokorelacji w resztach, więc uwzględniliśmy wszystkie autokorelacje, także w kwadratach reszt, a wszystkie parametry są istotne. Jedynie dopasowanie w grupie 3 w teście dobroci Pearsona jest nieco za słabe - prawdopodobnie rozkład reszt jest nieco inny niż wybrano "snorm", ale jest to duże odchylenie.

Dla upewnienia się czy nie ma gdzieś pomyłki, możemy użyć testu na DAC:

> # DAC w próbie
> dopas_hqc <- fitted(dopasGarch_hqc)
> stopaZwrotuZgodne <- dopas_hqc + residuals(dopasGarch_hqc)
> DACTest(forecast=dopas_hqc, actual=stopaZwrotuZgodne, test="PT")
$Test [1] "Pesaran and Timmermann" $Stat [1] 5.854 $p.value [1] 0.000000002401 $H0 [1] "Independently Distributed" $Decision [1] "Reject H0" $DirAcc [1] 0.71

 0,71 to tyle samo co wcześniej otrzymaliśmy (zmienna kryterium_hqc), więc wszystko się zgadza. Możemy dodatkowo sprawdzić jak wygląda DAC dla ostatnich 50 danych:

> Ost <- 50
> DACTest(forecast = tail(as.numeric(dopas_hqc), Ost), actual = tail(as.numeric(stopaZwrotuZgodne), Ost), test = "AG")
$Test
[1] "Anatolyev and Gerko"

$Stat
[1] 4.756

$p.value
[1] 0.0000009883

$H0
[1] "No Predictability"

$Decision
[1] "Reject  H0"

$DirAcc
[1] 0.78


Akurat tutaj użyłem drugiego rodzaju (AG), który daje najczęściej takie same wartości jak PT. W próbie uczącej model radzi sobie przez cały czas całkiem dobrze, ponieważ systematycznie hipoteza braku prognozowalności zostaje odrzucona.

Przejdziemy teraz stopniowo do próby testowej. Najpierw obliczamy DAC dla całej próby (ucząca + testowa). Estymujemy więc najpierw uzyskany dotychczas model dla całej próby za pomocą ugarchfit():

> dopasGarchCalaProba_hqc <- tryCatch(
+   {
+     ugarchfit(spec = mojGarch_hqc, data = stopaZwrotu, solver = "hybrid")
+   }, warning = function(w) {
+     message("Ostrzeżenie: ", conditionMessage(w))
+     ugarchfit(spec = mojGarch_hqc, data = stopaZwrotu, solver = "hybrid")
+   }
+ )
Ostrzeżenie: prawdopodobny problem zbieżności: 'optim' zwrócił kod= 1
Komunikat ostrzegawczy:
W poleceniu 'arima(data, order = c(modelinc[2], 0, modelinc[3]), include.mean = modelinc[1], ':
  prawdopodobny problem zbieżności: 'optim' zwrócił kod= 1
> dopasGarchCalaProba_hqc
*---------------------------------* * GARCH Model Fit * *---------------------------------* Conditional Variance Dynamics ----------------------------------- GARCH Model : eGARCH(1,1) Mean Model : ARFIMA(13,0,10) Distribution : snorm Optimal Parameters ------------------------------------ Estimate Std. Error t value Pr(>|t|) mu 0.000876 0.000000 3244.73 0 ar1 -0.527835 0.000196 -2694.35 0 ar2 0.684617 0.000062 11063.85 0 ar3 -0.142261 0.000398 -357.86 0 ar4 -0.820883 0.000061 -13496.80 0 ar5 0.293763 0.000021 14042.21 0 ar6 0.452103 0.000086 5227.05 0 ar7 -0.232687 0.000099 -2340.97 0 ar8 -0.070984 0.000108 -656.74 0 ar9 0.357832 0.000037 9762.27 0 ar10 0.020505 0.000048 427.04 0 ar11 -0.185835 0.000043 -4301.92 0 ar12 0.003256 0.000011 298.65 0 ar13 0.069660 0.000171 408.22 0 ma1 0.571280 0.000127 4507.56 0 ma2 -0.660489 0.000091 -7260.13 0 ma3 0.161063 0.000917 175.66 0 ma4 0.835642 0.000242 3456.81 0 ma5 -0.327134 0.000012 -26400.84 0 ma6 -0.538823 0.000154 -3503.44 0 ma7 0.192344 0.000114 1683.00 0 ma8 0.130609 0.000157 833.86 0 ma9 -0.234300 0.000108 -2175.97 0 ma10 -0.237432 0.000026 -9220.62 0 omega -0.765569 0.000111 -6870.67 0 alpha1 -0.254834 0.000028 -9192.27 0 beta1 0.893351 0.000117 7658.60 0 gamma1 -0.320261 0.000257 -1247.88 0 skew 0.815441 0.001631 500.05 0 Robust Standard Errors: Estimate Std. Error t value Pr(>|t|) mu 0.000876 0.000002 414.684 0 ar1 -0.527835 0.001263 -417.874 0 ar2 0.684617 0.000229 2984.727 0 ar3 -0.142261 0.001288 -110.408 0 ar4 -0.820883 0.000460 -1784.982 0 ar5 0.293763 0.000459 640.153 0 ar6 0.452103 0.000029 15504.063 0 ar7 -0.232687 0.000452 -514.996 0 ar8 -0.070984 0.000447 -158.903 0 ar9 0.357832 0.000662 540.440 0 ar10 0.020505 0.000189 108.288 0 ar11 -0.185835 0.000216 -861.202 0 ar12 0.003256 0.000051 64.180 0 ar13 0.069660 0.001136 61.315 0 ma1 0.571280 0.000962 593.674 0 ma2 -0.660489 0.000891 -741.613 0 ma3 0.161063 0.004921 32.733 0 ma4 0.835642 0.001461 571.974 0 ma5 -0.327134 0.000133 -2452.617 0 ma6 -0.538823 0.000077 -6953.614 0 ma7 0.192344 0.000347 554.317 0 ma8 0.130609 0.000598 218.284 0 ma9 -0.234300 0.000451 -519.405 0 ma10 -0.237432 0.000342 -694.913 0 omega -0.765569 0.000321 -2386.021 0 alpha1 -0.254834 0.000660 -386.185 0 beta1 0.893351 0.001459 612.425 0 gamma1 -0.320261 0.000021 -15025.018 0 skew 0.815441 0.004400 185.319 0 LogLikelihood : 512.7 Information Criteria ------------------------------------ Akaike -4.2057 Bayes -3.7722 Shibata -4.2330 Hannan-Quinn -4.0309 Weighted Ljung-Box Test on Standardized Residuals ------------------------------------ statistic p-value Lag[1] 1.434 0.2310 Lag[2*(p+q)+(p+q)-1][68] 25.530 1.0000 Lag[4*(p+q)+(p+q)-1][114] 42.082 0.9988 d.o.f=23 H0 : No serial correlation Weighted Ljung-Box Test on Standardized Squared Residuals ------------------------------------ statistic p-value Lag[1] 1.396 0.2374 Lag[2*(p+q)+(p+q)-1][5] 3.374 0.3430 Lag[4*(p+q)+(p+q)-1][9] 6.010 0.2971 d.o.f=2 Weighted ARCH LM Tests ------------------------------------ Statistic Shape Scale P-Value ARCH Lag[3] 2.544 0.500 2.000 0.1107 ARCH Lag[5] 3.226 1.440 1.667 0.2586 ARCH Lag[7] 5.209 2.315 1.543 0.2042 Nyblom stability test ------------------------------------ Joint Statistic: no.parameters>20 (not available) Individual Statistics: mu 0.02266 ar1 0.02024 ar2 0.02111 ar3 0.02888 ar4 0.01956 ar5 0.02272 ar6 0.02255 ar7 0.15612 ar8 0.02344 ar9 0.02388 ar10 0.02475 ar11 0.02319 ar12 0.02538 ar13 0.02162 ma1 0.02144 ma2 0.02181 ma3 0.03071 ma4 0.02657 ma5 0.02350 ma6 0.06516 ma7 0.02290 ma8 0.02397 ma9 0.02463 ma10 0.02744 omega 0.02269 alpha1 0.02531 beta1 0.18295 gamma1 0.06803 skew 0.03629 Asymptotic Critical Values (10% 5% 1%) Individual Statistic: 0.35 0.47 0.75 Sign Bias Test ------------------------------------ t-value prob sig Sign Bias 0.0987 0.9215 Negative Sign Bias 1.3704 0.1719 Positive Sign Bias 1.1653 0.2451 Joint Effect 3.2359 0.3567 Adjusted Pearson Goodness-of-Fit Test: ------------------------------------ group statistic p-value(g-1) 1 20 20.96 0.3392 2 30 37.65 0.1303 3 40 46.52 0.1903 4 50 46.96 0.5563

 

Kryteria informacyjne pogorszyły się. 

W końcu dochodzimy do prognozy kroczącej z reestymacją.

# Prognoza z reestymacją

klaster <- makeCluster(detectCores() - 1)

prognozaGarchReest_hqc <- ugarchroll(spec = mojGarch_hqc, data = stopaZwrotu, forecast.length = wyrzuc,

                                n.start = NULL, refit.every = 1, refit.window = "recursive",

                                window.size = NULL, solver = "hybrid", calculate.VaR = TRUE, 

                                VaR.alpha = c(0.01, 0.05),

                                cluster = klaster)

prognozaGarchReest_hqc <- resume(prognozaGarchReest_hqc)

stopCluster(klaster)


Jak widać znów stosujemy klaster. Funkcja resume() natomiast stosowana po to, by w razie braku konwergencji podczas obliczeń, jeszcze raz oszacować parametry przy pomocy innego algorytmu.

Po zakończeniu dostajemy obiekt, na którym możemy wykonać różne funkcje, jak plot(), jednakże lepiej samemu utworzyć taką funkcję, którą łatwiej sterować. Najpierw wyciągamy wartości oczekiwane:

prognozaReest_hqc <- prognozaGarchReest_hqc@forecast$density$Mu

i obliczymy jej DAC:

> DACTest(forecast=prognozaReest_hqc, 
+         actual=tail(as.numeric(stopaZwrotu), wyrzuc),
+         test="PT")
$Test
[1] "Pesaran and Timmermann"

$Stat
[1] 1.576

$p.value
[1] 0.05755

$H0
[1] "Independently Distributed"

$Decision
[1] "Fail to Reject H0"

$DirAcc
[1] 0.6333


Chociaż prognoza nie zdała testu prognozowalności kierunku, to otrzymane 63% poprawnych odpowiedzi jest to sporo więcej niż bez reestymacji (było 53%). Niemniej przy 30 danych (dla przypomnienia jest to prognoza próby testowej) musimy dostać min. 65%, żeby uznać wynik za coś więcej jak przypadek.

Wyciągnijmy daty utworzonej prognozy (są to daty próby testowej) i zamieńmy w xts:

datyTest <- as.Date(rownames(prognozaGarchReest_hqc@forecast$density))
 

Konwertujemy prognozę z reestymacją na obiekt xts:

prognozaReest_xts <- xts(prognozaReest_hqc, order.by = datyTest)


Mamy już prognozę z reestymacją, ale jeszcze potrzebujemy bez reestymacji, a więc ten sam kod, co stosowany poprzednio:

prognozaGarch_hqc <- ugarchforecast(dopasGarch_hqc, n.ahead = 1, n.roll = wyrzuc - 1)
prognoza_hqc <- as.numeric(fitted(prognozaGarch_hqc))
prognoza_xts <- xts(prognoza_hqc, order.by = datyTest)

Zauważmy tylko dwie rzeczy:
(1) n.roll zostało ustawione na 29 danych (wyrzuc = 30), dlatego że pierwsza prognoza nie dotyczy rolowania (n.ahead = 1). 
(2) obiekt xts dla prognozy bez reestymacji ma te same daty co z reestymacją, dlatego je tutaj znów wykorzystujemy. W tym miejscu jedna mała uwaga. Sama funkcja fitted() na obiekcie ugarchforecast:

> prognoza_hqc <- fitted(prognozaGarch_hqc)
> prognoza_hqc
    2024-03-15 2024-03-22 2024-03-28 2024-04-05 2024-04-12 2024-04-19 2024-04-26 2024-05-02
T+1   0.004522   -0.02022    0.02879   -0.01329     0.0068   0.001236   0.006209   -0.02933
    2024-05-10 2024-05-17 2024-05-24 2024-05-31 2024-06-07 2024-06-14 2024-06-21 2024-06-28
T+1   -0.02474    0.03675   -0.03833  -0.005701   -0.02924    0.02577   -0.01997    0.02392
    2024-07-05 2024-07-12 2024-07-19 2024-07-26 2024-08-02 2024-08-09 2024-08-16 2024-08-23
T+1   -0.01662   -0.06338    0.02571   -0.02336   -0.01967   -0.07568    0.09247   -0.05133
    2024-08-30 2024-09-06 2024-09-13 2024-09-20 2024-09-27 2024-10-04
T+1   0.009721   -0.02933   -0.04628    0.01524  -0.006271    0.05237


jest tablicą, której nazwy kolumn są datami, ale nie dotyczącymi okresu realizacji prognozy, tylko powstania prognozy. Nazwy wierszy dopiero pokazują okres, którego dotyczy prognoza. tak że "T+1"  oznacza datę z kolumny powiększoną o 1 okres, czyli tutaj jeden tydzień. Dlatego właśnie ostatnia data kończy się na 04.10.2024, a nie 11.10.2024. Stąd nie możemy bezpośrednio wyciągnąć stąd dat i potraktować jako daty prognozy. Inaczej sprawa wygląda w przypadku obiektu ugarchroll:

> prognozaGarchReest_hqc

*-------------------------------------*
*              GARCH Roll             *
*-------------------------------------*
No.Refits		: 30
Refit Horizon	: 1
No.Forecasts	: 30
GARCH Model		: eGARCH(1,1)
Distribution	: snorm 

Forecast Density:
                Mu  Sigma   Skew Shape Shape(GIG) Realized
2024-03-22  0.0045 0.0115 0.7211     0          0   0.0113
2024-03-28  0.0004 0.0195 0.9559     0          0   0.0242
2024-04-05  0.0055 0.0188 0.9982     0          0   0.0141
2024-04-12 -0.0166 0.0130 0.8980     0          0  -0.0100
2024-04-19  0.0207 0.0110 0.9934     0          0   0.0073
2024-04-26  0.0036 0.0229 0.7858     0          0   0.0061

..........................
                Mu  Sigma   Skew Shape Shape(GIG) Realized
2024-09-06  0.0229 0.0184 2.5994     0          0  -0.0397
2024-09-13 -0.0048 0.0302 0.8268     0          0   0.0079
2024-09-20 -0.0065 0.0379 0.7288     0          0  -0.0245
2024-09-27  0.0121 0.0353 0.8327     0          0   0.0442
2024-10-04  0.0075 0.0290 0.8249     0          0  -0.0391
2024-10-11  0.0048 0.0327 0.8784     0          0   0.0183

 
Tutaj ostatnia data faktycznie kończy się 11.10.2024, a prognozy są od razu porównywalne z prawdziwymi obserwacjami. Dlatego te daty mogliśmy od razu wyciągnąć i użyć.

Mamy zatem prognozę krocząca z reestymacją (prognozaReest_xts) oraz prognozę kroczącą bez reestymacji (prognoza_xts). Porównajmy je w końcu na wykresie. Aby zobaczyć je w szerszej perspektywie, dodamy też faktyczną stopę zwrotu:

plot(prognoza_xts, col = "blue", lwd = 2)
lines(prognozaReest_xts, col = "red", lwd = 2)
lines(stopaZwrotu, col="gray", lwd=2)
addLegend(legend.loc = "topleft", legend.names = c("Prognoza krocząca bez reestymacji", 
                                                   "Prognoza krocząca z reestymacją", 
                                                   "Tygodniowa stopa zwrotu WIG20"),
          lwd=3, bty="o", cex=1.05, col=c("blue", "red", "gray"))




Wygląda na to, że rzeczywiście prognoza reestymowana (czerwona) jest lepsza. 

Brakuje tu jeszcze prognozy na kolejny tydzień. Niestety ten artykuł długo piszę, więc dotyczyć będzie ona już kończącego się tygodnia. Chcemy stworzyć prognozę dla obu rodzajów. W przypadku niereestymowanym sytuacja jest prosta - wystarczy, że dodamy jeden okres do n.roll w ugarchforecast():

prognozaGarch_hqc <- ugarchforecast(dopasGarch_hqc, n.ahead = 1, n.roll = wyrzuc)
prognoza_hqc <- as.numeric(fitted(prognozaGarch_hqc))

Jednak zamiana na xts jest już troszkę trudniejsza, bo nie mamy nigdzie ustawionej daty z przyszłości:

czest <- "week"
nowaData <- seq(from = datyTest[length(datyTest)], by = czest, length.out = 2)[-1]
datyWykres <- c(datyTest, nowaData)
prognoza_xts <- xts(prognoza_hqc, order.by = datyWykres)

W przypadku prognozy reestymowanej nie ma możliwości bezpośredniego stworzenia prognozy poza próbę testową. Stąd należy po prostu wykonać zwykłą prognozę na całej próbie (tylko próba ucząca). Efekt będzie ten sam, bo to nic innego jak reestymacja po całej próbie:

prognozaGarchCalaProba_hqc <- ugarchforecast(dopasGarchCalaProba_hqc, n.ahead = 1)
prognozaCalaProba_hqc <- as.numeric(fitted(prognozaGarchCalaProba_hqc))
prognozaReest_xts <- xts(c(prognozaReest_hqc, prognozaCalaProba_hqc), 
                         order.by = datyWykres)
  
I wtedy powtarzamy kod na wykres:

plot(prognoza_xts, col = "blue", lwd = 2)
lines(prognozaReest_xts, col = "red", lwd = 2)
lines(stopaZwrotu, col="gray", lwd=2)
addLegend(legend.loc = "topleft", legend.names = c("Prognoza krocząca bez reestymacji", 
                                                   "Prognoza krocząca z reestymacją", 
                                                   "Tygodniowa stopa zwrotu WIG20"),
          lwd=3, bty="o", cex=1.05, col=c("blue", "red", "gray"))



Biorąc pod uwagę, że ta z reestymacją jest bardziej wiarygodna, w tym tygodniu bardziej prawdopodobny jest wzrost.


2. Wg DAC w próbie testowej

# budowa modelu wg DAC
mojGarch_DAC <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,
                                             garchOrder = c(1, 1)), 
                       mean.model = list(armaOrder = c(rzadAR_DAC, rzadMA_DAC), archm = czy_archm_DAC), 
                       distribution.model = rozklad)

dopasGarch_DAC <- ugarchfit(spec = mojGarch_DAC, data = stopaZwrotu, out.sample = wyrzuc, solver="hybrid")
dopasGarch_DAC

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: eGARCH(1,1)
Mean Model	: ARFIMA(12,0,8)
Distribution	: snorm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.001648    0.000000   4997.32        0
ar1    -0.663854    0.000077  -8578.89        0
ar2    -0.262185    0.000024 -10887.50        0
ar3    -0.051030    0.000004 -11599.32        0
ar4     0.309857    0.000032   9822.05        0
ar5     0.342351    0.000041   8407.70        0
ar6     0.359999    0.000043   8357.61        0
ar7     0.580322    0.000067   8602.21        0
ar8     0.010979    0.000001   8979.89        0
ar9     0.137562    0.000004  39254.41        0
ar10    0.026932    0.000003   7747.40        0
ar11    0.104187    0.000017   6313.71        0
ar12   -0.001168    0.000002   -514.54        0
ma1     0.921020    0.000120   7703.09        0
ma2     0.336957    0.000034  10038.87        0
ma3     0.097933    0.000013   7816.74        0
ma4    -0.370421    0.000018 -21002.82        0
ma5    -0.588972    0.000061  -9683.56        0
ma6    -0.596177    0.000061  -9768.79        0
ma7    -0.861720    0.000114  -7588.11        0
ma8    -0.046927    0.000005  -9285.27        0
archm  -0.029292    0.000004  -7966.21        0
omega  -0.808888    0.000100  -8114.56        0
alpha1 -0.311233    0.000853   -364.68        0
beta1   0.891072    0.000109   8198.37        0
gamma1 -0.364227    0.000340  -1070.40        0
skew    0.561526    0.002689    208.79        0

Robust Standard Errors:
        Estimate  Std. Error    t value Pr(>|t|)
mu      0.001648    0.000000  11903.419        0
ar1    -0.663854    0.000158  -4196.632        0
ar2    -0.262185    0.000105  -2494.059        0
ar3    -0.051030    0.000011  -4762.124        0
ar4     0.309857    0.000099   3118.566        0
ar5     0.342351    0.000054   6362.322        0
ar6     0.359999    0.000174   2072.837        0
ar7     0.580322    0.000191   3039.481        0
ar8     0.010979    0.000003   3309.169        0
ar9     0.137562    0.000056   2467.595        0
ar10    0.026932    0.000001  29139.146        0
ar11    0.104187    0.000020   5260.849        0
ar12   -0.001168    0.000003   -392.560        0
ma1     0.921020    0.000310   2975.374        0
ma2     0.336957    0.000039   8634.698        0
ma3     0.097933    0.000004  22486.872        0
ma4    -0.370421    0.000091  -4059.529        0
ma5    -0.588972    0.000149  -3952.367        0
ma6    -0.596177    0.000038 -15636.003        0
ma7    -0.861720    0.000496  -1736.513        0
ma8    -0.046927    0.000008  -6015.220        0
archm  -0.029292    0.000023  -1255.820        0
omega  -0.808888    0.000116  -6952.590        0
alpha1 -0.311233    0.002262   -137.605        0
beta1   0.891072    0.000602   1479.867        0
gamma1 -0.364227    0.000235  -1552.219        0
skew    0.561526    0.006557     85.641        0

LogLikelihood : 457 

Information Criteria
------------------------------------
                    
Akaike       -4.2996
Bayes        -3.8543
Shibata      -4.3305
Hannan-Quinn -4.1194

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                     0.02059  0.8859
Lag[2*(p+q)+(p+q)-1][59]  21.65058  1.0000
Lag[4*(p+q)+(p+q)-1][99]  35.98965  0.9985
d.o.f=20
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      0.565  0.4523
Lag[2*(p+q)+(p+q)-1][5]     5.352  0.1272
Lag[4*(p+q)+(p+q)-1][9]     6.283  0.2670
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]     2.216 0.500 2.000  0.1366
ARCH Lag[5]     2.468 1.440 1.667  0.3768
ARCH Lag[7]     2.517 2.315 1.543  0.6090

Nyblom stability test
------------------------------------
Joint Statistic:  no.parameters>20 (not available)
Individual Statistics:              
mu     0.02688
ar1    0.03818
ar2    0.03886
ar3    0.03903
ar4    0.03875
ar5    0.03790
ar6    0.03681
ar7    0.18916
ar8    0.03516
ar9    0.03328
ar10   0.03355
ar11   0.02755
ar12   0.03261
ma1    0.03305
ma2    0.03679
ma3    0.03647
ma4    0.03476
ma5    0.01836
ma6    0.03795
ma7    0.35048
ma8    0.03726
archm  0.02762
omega  0.02795
alpha1 0.01723
beta1  0.04913
gamma1 0.07816
skew   0.20351

Asymptotic Critical Values (10% 5% 1%)
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.4007 0.6890    
Negative Sign Bias  1.0956 0.2746    
Positive Sign Bias  0.4800 0.6317    
Joint Effect        1.4322 0.6980    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20      12.2       0.8769
2    30      37.0       0.1462
3    40      34.0       0.6970
4    50      34.0       0.9491


Model jest poprawny.

Cały kod jest taki sam, a jedynie tam gdzie trzeba, należy zastąpić końcówkę "_hqc" przez "_DAC". 
W samej próbie prognoza kierunku jest stosunkowo słaba, bo jedynie 63,5% poprawnych odpowiedzi, ale test odrzucił losowość. Ostatnie 50 danych prognozowanych jest już lepiej (68%), jednak w przypadku HQC było 78%. Dopiero w próbie testowej dostaniemy poprawę (prawie 77%). Pytanie, czy jest to prawdziwa poprawa? Aby odpowiedzieć na to pytanie, tworzymy prognozę z reestymacją.


# Prognoza z reestymacją
klaster <- makeCluster(detectCores() - 1)
prognozaGarchReest_DAC <- ugarchroll(spec = mojGarch_DAC, data = stopaZwrotu, forecast.length = wyrzuc,
                                     n.start = NULL, refit.every = 1, refit.window = "recursive",
                                     window.size = NULL, solver = "hybrid", calculate.VaR = TRUE, 
                                     VaR.alpha = c(0.01, 0.05),
                                     cluster = klaster)
prognozaGarchReest_DAC <- resume(prognozaGarchReest_DAC)
stopCluster(klaster)

prognozaReest_DAC <- prognozaGarchReest_DAC@forecast$density$Mu

DACTest(forecast=prognozaReest_DAC, 
        actual=tail(as.numeric(stopaZwrotu), wyrzuc),
        test="PT")

Niestety DAC spada do 60% tracąc istotność. Czyli w żadnym przypadku nie uzyskałem modelu prognozującego znak stopy zwrotu w sposób wiarygodny. Zanim to przeanalizuję, robię ten wykres dla modelu wg DAC (out of sample):



Oba typy prognoz wskazują, że ten tydzień będzie rosnący. Predykcja zostaje wzmocniona przez HQC, którego prognoza reestymująca również sugeruje wzrost.  Między prognozami prognozaReest_hqc i prognozaReest_DAC występuje dość silna korelacja na poziomie 0,67. Załóżmy więc, że strategia będzie polegać na tym, że sygnałem będzie sytuacja, gdy obie prognozy mają ten sam znak. Sprawdzamy DAC dla tej sytuacji:

te_same_znaki <- which(sign(prognozaReest_hqc) ==  sign(prognozaReest_DAC))
DACTest(forecast = prognozaReest_hqc[te_same_znaki], actual = as.numeric(tail(stopaZwrotu, wyrzuc))[te_same_znaki],
        test="PT")


Otrzymałem DAC = 65,2% (p-value = 0.028), co przemawia za zdolnością strategii do prognozowania kierunku lepiej niż rzucanie monetą. 

W tym miejscu warto nadmienić, że użycie drugiego z dostępnych testów DAC - "AG" (Anatolyev-Gerko) przynosi wynik na granicy istotności statystycznej:

> DACTest(forecast = prognozaReest_hqc[te_same_znaki], actual = as.numeric(tail(stopaZwrotu, wyrzuc))[te_same_znaki],
+         test="AG")
$Test
[1] "Anatolyev and Gerko"

$Stat
[1] 1.618

$p.value
[1] 0.05288

$H0
[1] "No Predictability"

$Decision
[1] "Fail to Reject H0"

$DirAcc
[1] 0.6522


Różnica między PT a AG polega na tym, że ten pierwszy bada zdolność prognostyczną znaku zmiennej, a ten drugi średniej warunkowej. Tak więc o ile PT skupia się na badaniu prawdopodobieństwa odgadnięcia znaku stopy zwrotu, o tyle AG szacuje siłę zdolności do prognozowania oczekiwanej stopy zwrotu [Anatolyev, S., Gerko, A. 2005, "A trading approach to testing for predictability, Journal of Business and Economic Statistics, 23(4), 455–461]. Jest więc oczywiste, że zawsze ten drugi będzie miał wyższe p-value, bo po prostu trudniej szacować średnią niż sam jej znak. Test AG ma interpretację bardziej ekonomiczną niż statystyczną, bo nawet jeśli jesteśmy w stanie prognozować znak stopy zwrotu, to nie znaczy jeszcze, że nasze zyski są ponadprzeciętne.

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.