Tytuł brzmi koślawo, bo słowo "krocząca" jest synonimem do "ruchomej" i np. średnia krocząca i średnia ruchoma oznaczają to samo. Z tym że ta średnia oznacza średnią z ostatnich x okresów, a prognoza jest zawsze do przodu. W ogóle określenie prognoza krocząca nie jest najlepsze i powinno się raczej mówić o prognozie nawijanej (rolującej) lub rekurencyjnej. W języku angielskim mówi się o recursive forecast. Natomiast rolling forecast zwykło się nazywać prognozę kroczącą z ruchomym oknem będącą analogią do średniej ruchomej. Tym przypadkiem zajmę się obecnie.
Krótko na temat samej funkcji w rugarch. Interesuje nas postać:
ugarchroll(spec, data, forecast.length = 500, n.start = NULL, refit.every = 25, refit.window = c("recursive", "moving"), window.size = NULL, solver = "hybrid", cluster = NULL)
gdzie:
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"), tzn. wybieramy tu prognozę nawijaną na kolejny okres, 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)
W tym badaniu ustawiam forecast.length = 30, refit.every = 3, refit.window = "moving", window.size = 300.
Temat ogólnie stanowi kontynuację prognozy za pomocą pakietu rugarch. Stanęło na tym, że modelowanie na dużej próbie jest utrudnione z uwagi na niestabilność parametrów wariancji warunkowej - przynajmniej dla tygodniowych stóp zwrotu. Wychwycił to specjalny test na proces TV-GARCH. Od teraz zawsze będę stosował ten test w pierwszej kolejności.
W końcu udało mi się znaleźć w miarę optymalną próbę - nie za małą i nie za dużą dla GARCH - jest to próba między 300 a 350 tygodniowych stóp zwrotu. Wybrałem 330, z czego ostatnie 30 będą out-of-sample.
Tym razem będę modelował WIG zamiast WIG20, a potem S&P 500 (źródło stooq.pl). Wszystkie potrzebne pakiety:
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")
if (require("tvgarch")==FALSE) {
install.packages("tvgarch")
}
library("tvgarch")
WIG
#Przygotowanie danych:
nazwaPliku = "wig_d.csv"
Obs <- 331
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
Tak jak powiedziałem, zaczynam od testowania TV-GARCH. Testuję ostatnie 330 stóp zwrotu WIG (zmienna stopaZwrotu):
# Sprawdzamy TV-GARCH
tvGarchTest <- tvgarchTest(y = as.numeric(stopaZwrotu))
tvGarchTest
Testing GARCH(1,1), i.e., the model under H0, against
TV-GARCH(1,1), i.e., the model under H1:
Estimation results for model under H0:
Model: GARCH(1,1)
intercept.h arch1 garch1
Estimate: 0.0002362 0.2062 0.51256
Std. Error: 0.0003463 0.1080 0.03947
Log-likelihood: 723.2
Transition variable in TV-GARCH(1,1): time
Results from the Non-Robust TR^2 Test:
NonRobTR2 p-value
H0:B3=B2=B1=0 2.4250 0.4890
H03:B3=0 1.4251 0.2326
H02:B2=0|B3=0 0.7273 0.3938
H01:B1=0|B3=B2=0 0.2775 0.5984
Results from the Robust TR^2 Test:
RobTR2 p-value
H0:B3=B2=B1=0 3.749 0.2898
H03:B3=0 1.122 0.2895
H02:B2=0|B3=0 2.245 0.1340
H01:B1=0|B3=B2=0 1.370 0.2419
Single
No. of locations (alpha = 0.05) 0
Żaden test nie sugeruje nawet, że może wystąpić niestabilność.
Jeśli chodzi o ogólną postać modelu wariancji, to przyjmuję EGARCH(1, 1), która okazywała się dotychczas najlepszą pod względem dopasowania. Kłopot zaczyna się dla części ARMA(p, q). Gdybym wiedział od początku jaką postać powinienem przyjąć, wtedy stosuję dwie funkcje: ugarchspec() oraz ugarchroll(). Nie jest już potrzebna ugarchfit(), którą częściowo zastępuje ugarchroll(). A nie jest potrzebna, bo normalnie służy nie tylko do utworzenia modelu empirycznego, ale i też do diagnostyki modelu. A skoro uznajemy, że prawidłową postać modelu znam od początku, to nie jest konieczna diagnostyka. Ale nawet jeżeli nie znam modelu, a interesuje mnie jedynie prognoza out-of-sample, to diagnostyka także jest zbędna. Może inaczej - nie jest niezbędna. Można jej użyć, ale tylko żeby sie dowiedzieć, czy model jest dobrze dopasowany do próby, przede wszystkim in-sample.
Funkcja ugarchfit() przyda się jednak do prognozowania naprawdę przyszłych, czyli nieistniejących jeszcze danych. To znaczy do samej prognozy użyję ugarchforecast(), której argumentem jest ugarchfit().
W sumie, tak jak poprzednio, ugarchroll() posłuży do testowania prognozy out-of-sample, a ugarchfit() + ugarchforecast() do prognozy na kolejne okresy.
Kryterium wyboru najlepszego modelu będzie się tym razem opierać wyłącznie na out-of-sample: kierunkowej dokładności (ang. directional accuracy criterion - DAC) oraz na sumie odchyleń (SO) od średniej prognozy. DAC musi być jak największe, SO jak najmniejsze. Zwrócę uwagę, że nie będę korygował tych miar o liczbę współczynników, jak to robią AIC, BIC i HQIC, dlatego należy traktować to badanie ad hoc.
Ostatnia sprawa to rozkład reszt. Tym razem dla przyspieszenia obliczeń przyjmę rozkład normalny, chociaż zdaję sobie sprawę, że mogę dostać przez to niedopasowaną wariancję warunkową. Jednak będzie to punkt wyjścia do dalszych analiz.
rozklad <- "norm"
hqc <- list()
model_Garch <- "eGARCH"
submodel_Garch <- NULL
wyrzuc <- 30
q <- 1
p <- 1
# 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)
clusterExport(klaster, varlist = c("stopaZwrotu", "rozklad", "model_Garch", "submodel_Garch", "q", "p", "wyrzuc"))
clusterEvalQ(klaster, library(rugarch))
# Utworzenie siatki parametrów dla pętli
siatka_parametrow <- expand.grid(k = 0:1, i = 0:20, j = 0:20)
# Iteracja po siatce parametrów bez foreach
wyniki <- lapply(1:nrow(siatka_parametrow), function(param) {
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({
prognoza <- NULL
# Define the GARCH specification
mojGarch <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,
garchOrder = c(q, p)),
mean.model = list(armaOrder = c(i, j), archm = as.logical(k)),
distribution.model = rozklad)
# Try ugarchroll and handle potential errors by attempting to use resume()
prognoza <- try(
ugarchroll(spec = mojGarch, data = stopaZwrotu, forecast.length = wyrzuc,
n.start = NULL, refit.every = 3, refit.window = "moving",
window.size = 300, solver = "hybrid", calculate.VaR = TRUE,
VaR.alpha = c(0.01, 0.05),
cluster = klaster),
silent = TRUE
)
# Check if ugarchroll returned an error
if (inherits(prognoza, "try-error")) {
message("Błąd w ugarchroll, przechodzę do resume(): ", attr(prognoza, "condition")$message)
return(list(error = attr(prognoza, "condition")$message, rzad_modelu = rzad_modelu))
}
# Use resume() if no error occurred
prognoza <- resume(prognoza)
# Extract forecasts and calculate metrics
prognozaSrednia <- prognoza@forecast$density$Mu # średnia prognoza
prognozaSigma <- prognoza@forecast$density$Sigma # warunkowe odchylenie st
stopaZwrotuZgodna <- prognoza@forecast$density$Realized # faktyczne stopy zwrotu
odchylenie <- sum((prognozaSrednia - stopaZwrotuZgodna)^2)
# Initialize DAC value
dac_wartosc <- NA
# Skip DACTest for ARMA(0, 0) with k = 0
if (i + j + k > 0) {
dac_wartosc <- DACTest(forecast = as.numeric(prognozaSrednia),
actual = as.numeric(stopaZwrotuZgodna),
test = "PT")$DirAcc # Liczenie DAC dla prognozy
}
return(list(
rzad_modelu = rzad_modelu,
dac = dac_wartosc,
odchylenie = odchylenie
))
}, error = function(e) {
# Catch block for unexpected errors
message("Błąd w ARMA(", i, ", ", j, ") z k = ", k, ": ", conditionMessage(e))
return(list(error = conditionMessage(e), rzad_modelu = rzad_modelu))
})
})
# Zatrzymaj klaster po przetwarzaniu
stopCluster(klaster)
dac <- list()
odchylenie <- list()
# Przetwórz wyniki i przechowaj je w oryginalnych listach
for (wynik in wyniki) {
if (!is.null(wynik)) {
for (nazwa_modelu in names(wynik)) {
dac[[wynik$rzad_modelu]] <- wynik$dac
odchylenie[[wynik$rzad_modelu]] <- wynik$odchylenie
#var5[[wynik$rzad_modelu]] <- var5
}
}
}
# DAC w próbie testowej
dacDf <- as.data.frame(do.call(cbind, dac))
rownames(dacDf) <- "DAC (out of sample)"
# odchylenie w próbie testowej
odchylenieDf <- as.data.frame(do.call(cbind, odchylenie))
rownames(odchylenieDf) <- "odchylenie (out of sample)"
porownanie <- bind_rows(dacDf, odchylenieDf)
# Znajdź wszystkie kolumny z maksymalną wartością w trzecim wierszu (DAC w próbie testowej)
max_dac <- max(round(porownanie[1, ], 4), na.rm = TRUE)
kryterium_dac <- porownanie[, which(round(porownanie[1, ], 4) == max_dac), drop = FALSE]
# Znajdź wszystkie kolumny z minimalną wartością w pierwszym wierszu (HQC)
min_odchylenie <- min(round(porownanie[2, ], 4), na.rm = TRUE)
kryterium_odchylenie <- porownanie[, which(round(porownanie[2, ], 4) == min_odchylenie), drop = FALSE]
kryterium_dac
kryterium_odchylenie
Z punktu widzenia DAC najlepszy okazał się ARMA(1, 3)-ARCH_In_Mean (DAC = 0,667 , SO = 0,196). Z punktu widzenia sumy odchyleń najlepszy model to ARMA(7, 8)-ARCH_In_Mean (DAC = 0,633, SO = 0,178).
Te dwa modele sprawdzimy.
Kryterium DAC
Wyciągamy parametry i ustalamy je w formie zmiennych:
# parametry - DAC
wyciagnij_liczby <- function(wyr, x) {
as.numeric(gsub("[^0-9]", "", substr(wyr, x, x+2), ))
}
rzadAR_DAC <- wyciagnij_liczby(colnames(kryterium_dac), 6)[1]
rzadMA_DAC <- wyciagnij_liczby(colnames(kryterium_dac), 9)[1]
czy_archm_DAC <- grepl("ARCH_In_Mean", colnames(kryterium_dac))[1]
Specyfikacja otrzymanego modelu:
mojGarch_DAC <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,
garchOrder = c(p, q)),
mean.model = list(armaOrder = c(rzadAR_DAC, rzadMA_DAC), archm = czy_archm_DAC),
distribution.model = rozklad)
Następnie stosujemy właśnie prognozę kroczącą - rolującą (robię to drugi raz, pierwszy raz w pętli tylko w celu wyszukania optimum, teraz pełna analiza - chociaż można by to wszystko zapisać w pętli).
# Prognoza z reestymacją - prognoza krocząca
klaster <- makeCluster(detectCores() - 1)
prognozaGarch_DAC <- ugarchroll(spec = mojGarch_DAC, data = stopaZwrotu, forecast.length = wyrzuc,
n.start = NULL, refit.every = 3, refit.window = "moving",
window.size = 300, solver = "hybrid", calculate.VaR = TRUE,
VaR.alpha = c(0.01, 0.05),
cluster = klaster)
prognozaGarch_DAC <- resume(prognozaGarch_DAC)
stopCluster(klaster)
prognozaMu_DAC <- prognozaGarch_DAC@forecast$density$Mu
prognozaSigma_DAC <- prognozaGarch_DAC@forecast$density$Sigma
Aby sprawdzić czy rzeczywiście dostałem ten model co chciałem, robię test DAC:
DACTest(forecast=prognozaReest_DAC,
actual=tail(as.numeric(stopaZwrotu), wyrzuc),
test="PT")
Stat
[1] 1.975
$p.value
[1] 0.02411
$H0
[1] "Independently Distributed"
$Decision
[1] "Reject H0"
$DirAcc
[1] 0.6667
Czyli rzeczywiście to co powinno być.
Ostatnia część analityczna to utworzenie prognozy na przyszłe okresy, których jeszcze nie ma. Wykorzystamy całą próbę i funkcje ugarchfit() i ugarchforecast():
# Prognoza na podstawie całej próby na 1 przyszły okres
dopasGarchCalaProba_DAC <- tryCatch(
{
ugarchfit(spec = mojGarch_DAC, data = stopaZwrotu, solver = "hybrid")
}, warning = function(w) {
message("Ostrzeżenie: ", conditionMessage(w))
ugarchfit(spec = mojGarch_DAC, data = stopaZwrotu, solver = "hybrid")
}
)
Normalnie w tym miejscu przeprowadziłbym diagnostykę czyli wpisał dopasGarchCalaProba_DAC, ale tego nie robię, bo nie jest to teraz najważniejsze.
Tworzę prognozę na 1 okres do przodu zarówno dla średniej jak i wariancji warunkowej (ryzyka):
prognozaGarchNowa_DAC <- ugarchforecast(dopasGarchCalaProba_DAC, n.ahead = 1)
# średnia
prognozaMuNowa_DAC <- as.numeric(fitted(prognozaGarchNowa_DAC))
# ryzyko
prognozaSigmaNowa_DAC <- as.numeric(sigma(prognozaGarchNowa_DAC))
Dalsza część to już dostosowanie do wykresu:
# Dla wykresu - tworzenie daty, której jeszcze nie ma - 1 okres na przyszłość
datyTest <- tail(index(stopaZwrotu), wyrzuc)
czest <- "week"
nowaData <- seq(from = datyTest[length(datyTest)], by = czest, length.out = 2)[-1]
# Dla wykresu - cały zakres dat
datyWykres <- c(datyTest, nowaData)
# Dla wykresu - średnia prognoza i ryzyko, cały zakres (z reestymacją + na przyszłość)
prognozaMu_xts <- xts(c(prognozaMu_DAC, prognozaMuNowa_DAC), order.by = datyWykres)
prognozaSigma_xts <- xts(c(prognozaSigma_DAC, prognozaSigmaNowa_DAC), order.by = datyWykres)
# wykres
plot(prognozaMu_xts, col = "blue", lwd = 2)
lines(prognozaSigma_xts, col = "red", lwd = 2)
lines(-prognozaSigma_xts, col = "red", lwd = 2)
lines(tail(stopaZwrotu, wyrzuc), col="gray", lwd=2)
addLegend(legend.loc = "topleft", legend.names = c("Średnia prognoza krocząca",
"Średnie ryzyko prognozy kroczącej",
"Tygodniowa stopa zwrotu WIG"),
lwd=3, bty="o", cex=0.85, col=c("blue", "red", "gray"))
Sporo błędów, model DAC = 67%, ale nie widać tego specjalnie, szczególnie w ostatnich tygodniach.
Kryterium sumy odchyleń
Kod jest analogiczny, jedynie zamiast kryterium_DAC wstawiamy kryterium_odchylenie i następnie tak samo z pozostałymi zmiennymi, nie będę więc już powtarzał kodu.
Z pozoru dziwne się może wydać, że tak dużo odchyleń wychodzi poza obszar oczekiwanego ryzyka. Niestety może być to cena przyjęcia rozkładu normalnego dla przyspieszenia obliczeń.
Uśrednione prognozy
Oba kryteria można uśrednić. Wtedy dostaniemy:
Widać, że to uśrednienie ma sens - tam gdzie kryterium DAC było gorsze, tam lepiej wypadało kryterium odchyleń. Można powiedzieć, że prognoza na następny tydzień jest delikatnie wzrostowa, może wokół zera.
Natomiast ryzyka nie udało się tu dobrze zamodelować, przypuszczalnie z powodu tego nieszczęsnego rozkładu normalnego.
S&P 500
Jedyną zmianą w kodzie jest nazwaPliku = "^spx_d.csv".
Kryterium DAC
Otrzymałem 2 modele: ARMA(12, 20) (DAC = 0,8 i SO = 0,008) oraz ARMA(9, 14)-ARCH_In_Mean (DAC = 0,8 i SO = 0,0081).
ARMA(12, 20)
ARMA(9, 14)-ARCH_In_Mean
Kryterium sumy odchyleń (SO)
W tym wypadku najlepsze okazały się ARMA(11, 7)-ARCH_In_Mean (DAC = 0,67 i SO = 0,0079) oraz ARMA(7, 9)-ARCH_In_Mean (DAC = 0,67 i SO = 0,00791)
ARMA(11, 7)-ARCH_In_Mean
ARMA(7, 9)-ARCH_In_Mean
Uśrednione prognozy
Wszystkie 4 prognozy możemy uśrednić i dostajemy:
W sumie oczekiwać można spadków amerykańskiego indeksu w następnym tygodniu przy średnim ryzyku. Niemniej trzeba mieć na względzie, że badanie zawiera wiele uproszczeń, jak nieuwzględnienie poprawnego rozkładu reszt, brak korekty o liczbę współczynników modelu oraz wnioskowanie na podstawie zaledwie 30 ostatnich obserwacji, tzn. out-of-sample, wyłączając dopasowanie in-sample.
