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 więc "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:

plot(prognozaReest_hqc, type = "l", col = "blue", lwd = 2, xaxt = "n")
lines(prognozaReest_DAC, col = "darkred", lwd = 2)
axis(side = 1, at = index(prognozaReest_hqc), labels = datyTest, cex.axis = 0.8)
abline(h=0, col="grey")
grid(nx=NULL, ny=NULL, lwd=3)
legend(x="top", legend=c("Prognoza z reestymacją", "Prognoza bez reestymacji"), lty=c(1,1),
       lwd=c(2,2), col=c("blue", "darkred"), cex = 0.7, bg = "white")
korelacja <- round(cor(prognozaReest_hqc, prognozaReest_DAC), 2)
usr <- par("usr")
text(x = usr[2] - 4, y = usr[3]+0.005, paste("Korelacja =", korelacja))



Załóżmy, ż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.

Brak komentarzy:

Prześlij komentarz