Prognoza krocząca z ruchomym oknem w rugarch (język R)

 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 do 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 się dowiedzieć, czy model jest dobrze dopasowany do próby 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 modelu, 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.

Klin zwyżkujący na NASDAQ

 Energia z jaką NASDAQ pnie się w góre jest porażająca i jednocześnie groźna. To wygląda na klin zwyżkujący, który kończy się zazwyczaj zmianą trendu:

STS i RSI też znajdują się w miejscach wykupienia, ten pierwszy na skrajnie wysokim poziomie:

Punkt kulminacyjny na S&P500

 Obecnie amerykański indeks osiąga poziom, który będzie psychologicznie wpływał na wycofywanie się traderów. Oczekuję w końcu tej bessy.

Trading Economics od dłuższego już czasu prognozuje podobnie:

Straszenie Trumpem już nabiera absurdalnych rozmiarów

Dzień po wyborach dziennikarze już wiedzą, co Trump zrobi z Ukrainą, ale to dopiero początek straszenia. Nie minął tydzień, a poziom straszenia osiąga taki poziom, że dziennikarze nawet nie zauważają własnego absurdu. Nie absurdalnego toku rozumowania - tylko absurdu - bo rozumowania tu nie ma. Business Insider, czyli Onet, napisał artykuł To może być pierwsza bitwa Trumpa. Jeszcze nigdy żaden prezydent USA się nie odważył, w którym przytacza słowa eksperta o tym, że Trump "nie będzie miał żadnych problemów, żeby ewentualnie przegłosować jakąś nowelizację ustawy o banku centralnym", "na podstawie obecnie obowiązującego prawa on mógłby wymienić właściwie całą Radę Gubernatorów, łącznie z przewodniczącym Rady Gubernatorów (...), żeby wszystkie osoby, które tam będą zasiadały, były osobami z jego nominacji i z tego grona być może chciałby wskazać nowego przewodniczącego Rady ".

A w następnym akapicie: 

To rodzi bardzo poważne konsekwencje i to może być nawet zagrożeniem dla całego systemu finansowego Stanów Zjednoczonych, ponieważ bank centralny USA cieszy się nimbem autonomii. Jako jedyna instytucja federalna ma swój własny budżet i nie musi się kongresu prosić o pieniądze.  

Dotychczas też żaden prezydent nie ośmielił się ingerować w pracę instytucji. — Odkąd przyjęto ustawę o rezerwie federalnej, z jednym wyjątkiem podczas II Wojny Światowej, wszyscy prezydenci szanowali autonomię banku centralnego. A Trump ma takie pomysły i wcale się z tym nie kryje. Szedł do tych wyborów właśnie z hasłami, że on po prostu będzie w stanie lepiej przewidzieć, jakie stopy procentowe powinny być.

Czyli tak, najpierw ekspert twierdzi, że wszystko byłoby zgodne z prawem, a następnie dodaje, że to zagrożenie dla całego systemu finansowego USA. Dobre. Nie wiedziałem, że mają tam dyktaturę i jeden człowiek może rozwalić cały system.

Ekspert Onetu śmieje się, że Trump wie lepiej od całego FED-u, jakie powinny być stopy proc. A przecież wiadomo, że FED utrzymywał za długo zerowe stopy, doprowadzając w 2022 r. do 9% inflacji. To ja w 2021 r. (zob. tu) pokazałem, że stopa referencyjna (effective federal funds rate) została zaniżona o ok. 2 pkt proc. Już wtedy inflacja wzrosła do ponad 5%. A skoro ja to potrafiłem, to tym bardziej Trump ma analityków, którzy to potrafili. Jeżeli jakaś instytucja się skompromitowała, to ją należy krytykować, a nie miliardera, który krytykuje instytucję.

Mało tego, mam wątpliwości czy w ogóle FED jest taki potrzebny. W artykule  FED do likwidacji?! dowiodłem, że wpływ FEDu na gospodarkę USA był marniutki, właściwie negatywny.

W ostatnich latach mówi się często o niebezpieczeństwie podważania wszelkich autorytetów, szczególnie naukowych, i Trumpa wrzuca się w ten wszechobecny trend. Fakty jednak są takie, że w niektórych przypadkach to podważanie jest niezbędne, aby wywołać jakąś dyskusję na temat jakości tych instytucji. Jeżeli nie spełniają swoich funkcji (celów), to albo należy je poprawić, albo zlikwidować.

Czy indeksy GPW są źle liczone?

Poprawne porównywanie stóp zwrotu z różnych indeksów, szczególnie w dłuższych okresach, wymaga odpowiedniej wiedzy. Np. WIG jest indeksem dochodowym w tym sensie, że uwzględnia nie tylko zmiany cen akcji, ale też dywidendy. Z kolei WIG20 i S&P 500 są indeksami cenowymi w tym sensie, że jedynie pokazują zmiany cen akcji, stąd nie są korygowane o dywidendy. Stanowią po prostu ważone kapitalizacje odpowiednich spółek. Aby móc porównać WIG z WIG20, stworzono WIG20TR - czyli właśnie odpowiednik WIGu dla największych 20 spółek. Również S&P 500 ma swój odpowiednik indeksu dochodowego - S&P 500 Total Return.

Tak więc porównanie WIG z S&P 500 nie jest poprawne - należy porównać go z S&P 500 TR. Można powiedzieć, że Total Return (TR) oznacza zwrot z indeksu dochodowego. Wydaje się więc, że sposób obliczania powinien być co do zasady ten sam - przynajmniej w kwestii dodawania dywidend.

I tu dochodzimy do meritum. Przeanalizowałem wzory, jakie stosują GPW i NYSE. Wg mnie sposób obliczania S&P 500 TR jest poprawny, a WIG20TR (a więc i WIG) niepoprawny.

Total Return na GPW

Metodologię obliczania indeksów, w tym WIG20TR, znajdziemy w regulaminach na tej stronie GPW -> W Regulaminie Rodziny Indeksów Giełdowych. 

Ogólny wzór na każdy indeks GPW jest następujący (pkt 4.2.1. Regulaminu):

gdzie:

Indeks(t) ‒ wartość Indeksu na Sesji „t”

M(t) ‒ kapitalizacja Portfela Indeksu na Sesji „t”

M(0) ‒ kapitalizacja Portfela Indeksu w Dniu Bazowym

Indeks(0) ‒ wartość Indeksu w Dniu Bazowym

K(t) ‒ Współczynnik Korygujący Indeksu na Sesji „t”

Dla WIG20TR (pkt 5.3.4.): 

gdzie:

M(t’) ‒ kapitalizacja Portfela Indeksu po rewizji, korekcie lub korekcie nadzwyczajnej obliczana na podstawie ostatniego kursu

M(t) ‒ kapitalizacja Portfela Indeksu zamknięcia przed rewizją, korektą lub korektą nadzwyczajną 

Dla uproszczenia załóżmy, że prawa poboru i nowe emisje akcji nie występują. Wtedy:

D(t) ‒ wartość dywidendy z akcji lub wartość teoretyczna dywidendy z akcji (obliczana przez Giełdę zgodnie ze Szczegółowymi Zasadami Obrotu Giełdowego), które na Sesji „t+1” po raz pierwszy będą notowane „bez dywidendy”; w przypadku ustalenia dywidendy w walucie obcej, kwota dywidendy zostaje przeliczona na złote.

Wyjaśnienie: w mianowniku występuje korekta K(t), która jest ułamkiem (indeksu bez dywidendy do indeksu z dywidendą), a więc podnosi cały indeks o odwrotność tego ułamka. Stąd WIG20TR nie będzie spadał o dywidendę. 

Ta konstrukcja jest identyczna dla WIG. Z kolei WIG20 nie zawiera owej korekty (zob. pkt 5.2.5. Regulaminu) i dlatego będzie spadał o średnioważoną dywidendę.

 

Total Return na NYSE i NASDAQ

Metodologia obliczenia TR jest zawarta w dokumencie Index Mathematics Methodology -> dok. pdf -> str. 21. 

Wtedy dzienny całkowity zwrot (Daily Total Return - DTR) zapiszemy:

gdzie IndexLevel(t) to typowy indeks giełdowy, jak S&P 500.

Wówczas indeks całkowitego zwrotu będzie określony wzorem:

Wyjaśnienie: Indeks z poprzedniego dnia jest powiększony o zwrot z tytułu zmiany cen i otrzymanej dywidendy nowego dnia.

Czy obydwa indeksy są równe?

Aby odpowiedzieć na to pytanie, podstawmy korektę K(t) do wzoru na Indeks:

Za K(t-1) analogicznie:

Powtarzając ten proces aż do K(0) dostaniemy:

W okresie 0 nie ma korekty albo przyjmujemy, że wszystkie wartości poprzednie M(t<0) = M(t<0') = M(0), stąd K(0) = 1. Dalej, w uproszczeniu możemy powiedzieć, że M(t'), M(t-1'), ... są to wartości na otwarcie kolejnego dnia, podczas gdy M(t), M(t-1), ... to wartości na zamknięcie tego samego dnia. Przy takim wyobrażeniu możemy zapisać, że:

dzięki czemu wzór na indeks się skraca:

Wyprowadźmy wzór na r'. Wiemy, że

M(t) = M(t') + D(t)

Czyli dla t-1:

M(t-1) = M(t-1') + D(t-1)

Teraz przekształcamy:

M(t) = M(t-1) + [M(t) - M(t-1)] = M(t-1') + D(t-1) + M(t) - M(t-1') - D(t-1) =

= M(t-1')*[M(t-1') / M(t-1') + [D(t-1) + M(t) - M(t-1') - D(t-1)]/M(t-1')] = 

= M(t-1')*[1 + [M(t) - M(t-1')]/M(t-1')] = M(t-1')*(1 + r')

Mamy więc

M(t) = M(t-1')*(1 + r'), gdzie

r(t') = [M(t) - M(t-1')]/M(t-1') = M(t) / M(t-1') - 1

Czyli zauważmy, że 

Jest to stopa zwrotu z WIG i WIG20TR. Problem polega na tym, że to wygląda tak, jakby ceną na zamknięcie była cena po odcięciu dywidendy - chociaż jeszcze tego odcięcia nie było, ponieważ nastąpi dopiero w kolejnym dniu. I wtedy w tym kolejnym dniu cena M(t) jest odnoszona do tej sztucznej ceny niby po odcięciu. Jak dla mnie to nie ma sensu. 

M(t) możemy też tak przekształcić, by uzyskać prostą stopę zwrotu:

M(t) = M(t-1) + [M(t) - M(t-1)] = M(t-1)/M(t-1)*M(t-1) + [M(t) - M(t-1)]/M(t-1)*M(t-1)   =

= M(t-1)*(1+r),

r(t) = [M(t) - M(t-1)]/M(t-1) = M(t) / M(t-1) - 1

r(t) jest zwykłą stopą zwrotu, typową dla WIG20 lub S&P 500, czyli bez dywidend. W takim razie całkowita dzienna stopa zwrotu z S&P 500 musi być równa:

Jest to oczywiście ten sam wzór, który wcześniej był podany w metodologii obliczania S&P.

Wynika z tego, że w ogólnym przypadku TR dla WIG nie będzie równa TR dla S&P 500, a dokładniej dla WIG będzie większa! Dla WIG dywidendę odejmuje się w mianowniku i od tego punktu liczy się wzrost. Dla S&P 500 TR w liczniku dodaje się dywidendę w stosunku do zamknięcia poprzedniego dnia. Z powodu tej różnicy powstaje przesunięcie czasowe w ustaleniu okresu dywidendy. W przypadku WIG i WIG20TR dywidenda nie dotyczy dnia, w którym nastąpiło jej odcięcie od cen, ale od poprzedniego dnia. W przypadku S&P 500 TR dywidenda dotyczy tego dnia, w którym następuje odcięcie. 

Przykład

Żeby zobaczyć różnicę zrobimy sztuczny przykład. W dniu 0 kupujemy 1 akcję po 10 zł z prawem do dywidendy po 1 zł. Jutro, w dniu 1, dostajemy dywidendę. Kurs akcji spadnie o wielkość dywidendy (bez odcięcia wyszedłby na zero). W dniu 2 kurs rośnie o 15%, a dodatkowo ustalona zostaje znów dywidenda 1 zł. W dniu 3 cena na otwarcie spada o wielkość tej dywidendy, ale od otwarcia cena rośnie o 10%:

(1) WIG / WIG20TR:

dzień 0: 

D(0) = 1

M(0) = 10

M(0') = M(0) - D(0) = 10 - 1 = 9

K(0) = 1

indeks(0) = M(0) = 10

dzień 1: 

D(1) = 0

M(1) = M(0') = 9

M(1') = M(1) - D(1) = 9 - 0 = 9

K(1) = M(0') / M(0) * K(0) = 9 / 10 * 1 = 0,9

indeks(1) = M(1) / (M(0) * K(1)) * 10 = 9 / (10 * 0,9) * 10 = 10

DTR(1) = indeks(1) / indeks(0) - 1 = 10 / 10 - 1 = 0

Sprawdźmy czy DTR(1) jest równe r':

r' = M(1) / M(0') - 1 = 9 / 9 - 1 = 0

Zgadza się: r' = DTR(1)

dzień 2: 

D(2) = 1

M(2) = M(1')*1,15 = 9*1,15 = 10.35

M(2') = M(2) - D(2) = 10,35 - 1 = 9,35

K(2) = M(1') / M(1) * K(1) = 9 / 9 * 0,9 = 0,9

indeks(2) = M(2) / (M(0) * K(2)) * 10 = 10,35 / (10 * 0,9) * 10 = 11,5

DTR(2) = indeks(2) / indeks(1) - 1 = 11,5 / 10 = 0,15

Sprawdźmy czy DTR(2) jest równe r':

r' = M(2) / M(1') - 1 = 10.35 / 9 - 1 = 0,15

Zgadza się: r' = DTR(2)

dzień 3: 

D(3) = 0

M(3) = M(2')*1,1 = 9,35*1,1 = 10,29

M(3') = M(3) - D(3) = 10,29 - 0 = 10,29

K(3) = M(2') / M(2) * K(2) = 9,35 / 10,35 * 0,9 = 0,813

indeks(3) = M(3) / (M(0) * K(3)) * 10 = 10,29 / (10 * 0,813) * 10 = 12,66

Policzmy zmianę w trzecim dniu:

DTR(3) = indeks(3) / indeks(2) = 12,66 / 11,5 - 1 = 0,10

Sprawdźmy czy DTR(3) jest równe r':

r' = M(3) / M(2') - 1 = 10,29 / 9,35 - 1 = 0,10

Zgadza się: r' = DTR(3)

W okresie od 0 do 3-go dnia indeks wzrósł o:

TR = (12,66 / 10 - 1)*100% = 26,6%

(2) S&P 500:

dzień 0: 

D(0) = 0

M(0) = 10

indeks(0) = M(0) = 10

dzień 1:  

D(1) = 1

M(1) = M(0) - D(1) = 10 - 1 = 9

DTR(1) = (M(1) + D(1)) / M(0) - 1 = (9 + 1) / 10 - 1 = 0

indeks(1) = indeks(0)*(1 + DTR(1)) = 10*(1 + 0) = 10

Zmiana indeksu:

indeks(1) / indeks(0) - 1 = 10/10 - 1 = 0 = DTR(1) 

dzień 2: 

D(2) = 0

M(2) = (M(1) - D(2))*1,15 = 9*1,15 = 10,35

DTR(2) = (M(2) + D(2)) / M(1) - 1 = (10,35 + 0) / 9 - 1 = 0,15

indeks(2) = indeks(1)*(1 + DTR(2)) = 10*(1 + 0,15) = 11,5

dzień 3: 

D(3) = 1

M(3) = (M(2) - D(3))*1,1 = (10,35 - 1)*1,1 = 9,35*1,1 = 10,29

DTR(3) = (M(3) + D(3)) / M(2) - 1 = (10,29 + 1) / 10,35 - 1 = 0,091

indeks(3) = indeks(2)*(1 + DTR(3)) = 11,5*(1 + 0,091) = 12,55

W okresie od 0 do 3-go dnia indeks wzrósł o:

TR = (12,55 / 10 - 1)*100% = 25,5%

Zatem widzimy, że - mimo identycznych założeń: poziomów kapitalizacji M(t), dziennych zmian ceny i identycznych dywidend - WIG rośnie mocniej od S&P 500 TR o 1,1 pkt proc.

Podsumowanie

No więc pytanie czy to możliwe, żeby GPW źle obliczała WIG i WIG20TR oraz inne indeksy TR? Są tylko 3 możliwości. Albo czegoś nie rozumiem, albo regulamin nie jest wystarczająco precyzyjny i GPW dokonuje jeszcze dodatkowej korekty indeksu, albo naprawdę źle liczą ten indeks. W USA jest natomiast prawidłowo - indeks TR jest konstruowany na podstawie powszechnie stosowanego wzoru na całkowitą stopę zwrotu.  

Mimo wszystko nawet jeżeli mam rację, to nie ma to w rzeczywistości wielkiego wpływu, bo mało prawdopodobne, żeby kurs zmieniał się o 10% i więcej dziennie, jak w przykładzie. W dodatku sam S&P 500 - bez dywidend - radzi sobie dużo, dużo lepiej niż nasz WIG z dywidendami. 

Lepiej żeby Trump wygrał?

Już 5 listopada odbędą się wybory w USA. Prof. A. Lichtman, którego metoda zawiodła tylko raz w 9-ciu na 10 przypadków (w ciągu ostatnich 40 lat) - przewiduje wygraną Kamali Harris (zob. np. tu , tu i tu). Przy tym trzeba zaznaczyć, że model Lichtmana jest krytykowany za zbytni subiektywizm. Mnie też zastanawia mocna krytyka Trumpa z jego strony (ostatnio razem z synem nazwał go faszystą, zob. tu). Oczywiście można chwalić za szczerość, ale wolałbym, żeby mówił to już po wyborach - byłby dla mnie bardziej wiarygodny. 

W artykule Dlaczego Demokraci dają lepsze stopy zwrotu niż Republikanie? pokazałem, że podczas rządów demokratów stopy zwrotu z S&P 500 są istotnie wyższe. Od tego czasu jednak indeks urósł 15%. W dodatku za nami 2 lata potężnej hossy (wzrost o ponad 40%). Gdyby więc Trump wygrał, wszystko stałoby się jasne: w Polsce indeksy by runęły (zresztą już chyba dyskontowanie Trumpa następuje), a w USA bessa mogłaby się rozpędzić. Dzięki temu nie tylko można by zarobić na bardziej pewnych spadkach, ale też łatwiej byłoby stopniowo budować portfel na przyszłość.

Spójrzmy na taki Microsoft:

Zapowiada się głębsze obsunięcie. Obecna wycena jest naprawdę okej (jak wyceniałem z rok temu, to dostałem ok. 400$), ale skoro poszybowało do 460, to i "musi" spaść do 340 dla równowagi.

Trochę przypomina to moment rozpoczęcia wojny w 2022 r., kiedy ludzie bali się wojny z Rosją. Wtedy bessa wydawała mi się zbyt oczywista - sądziłem, że spadki będą głębokie, ale krótkie. Teraz byłoby to dużo bardziej naturalne: dojście Trumpa do władzy wywoła strach nie tylko w Polsce, ale nie będzie tak gwałtowny jak wtedy w 2022 r. Wprawdzie wtedy USA przygotowywało publikę przez wiele miesięcy na wojnę, ale wielu (może nawet większość) w to nie wierzyło. Teraz jest niemalże na odwrót, media straszą, że jak Trump dojdzie do władzy, to za chwilę zobaczymy rosyjskie czołgi przy granicy. Więc będzie to znowu dobra okazja do kupowania co niektórych akcji. Portfel 60% największych spółek S&P 500, 20% kilka spółek z WIG i 20% coś innego - ETF, futures (np. na WIG20), fundusz, zagraniczne akcje, może surowce. To pierwsze przybliżenie portfela, jaki bym dziś budował.

Lichtman twierdzi, że Harris wygra, ale ekonomiczne argumenty idą na przekór. W USA jest niskie bezrobocie i niska inflacja - a to skłania ludzi do wyboru republikanów, jak było wyjaśnione we wspomnianym artykule. Poza tym jest jeszcze coś, o czym nikt nie mówi. Harris jest kobietą, a więc byłaby pierwszym prezydentem USA płci żeńskiej. Przynajmniej nasze media jakoś w ogóle nie poruszają tej kwestii, która przecież jest istotna. Bo skoro nigdy kobieta nie została prezydentem (w lewicowej nowomowie będzie się na pewno mówić "prezydentka"), to znaczy, że kobietom jest trudniej wygrać. Biorąc to wszystko pod uwagę, nastawiam się raczej na powrót Trumpa. To będzie prawdziwy test dla metody Lichtmana - jeśli znowu będzie miał rację, to zamknie usta krytykom. Jeśli się pomyli, to nikt już nie będzie go traktował poważnie i stanie się rodzajem anegdoty.

Za mała próba czy zły model?

 I znowu model zawiódł. Prawdopodobieństwo wzrostu WIG20 w poprzednim tygodniu szacowałem na ponad 0,6. Jednak indeks spadł. Diagnostyka nie wykryła jakichś nieprawidłowości, szczególnie model wybrany zgodnie z kryterium DAC w próbie testowej 30 obserwacji wydawał się całkiem poprawny. Doszedłem do wniosku, że problem może leżeć w zbyt niskim DAC po reestymacji. Wówczas dostałem 63,5% poprawnych odpowiedzi, co było niewystarczające dla 30 danych (nieistotne), tak że ciągle był to raczej przypadek. O tym żeby uznać to za coś więcej, przekonała mnie korelacja 0,67 między reestymowanymi prognozami kroczącymi na podstawie różnych kryteriów (HQC i DAC-of-sample). Mógł to być oczywiście pech, ale trzeba wziąć pod uwagę inne powody.

Pierwsza sprawa to wielkość próby. Użyte 200 obserwacji w próbie uczącej to nie jest za dużo jak na taki model. Wybrałem 200 tygodni ad hoc na takiej zasadzie, żeby z jednej strony obliczenia nie trwały za długo, a z drugiej, żeby uwzględnić przynajmniej jeden cykl giełdowy, który mniej więcej tyle trwa. Może być to jednak za mało. Z drugiej strony rodzi się problem reestymacji, który szczegółowo omawiałem poprzednio. Jeżeli uznajemy, że prawdziwy model nie istnieje i należy przez cały czas go reestymować, to lepiej aby próba testowa nie była za duża. Nie chodzi nawet o czasochłonność obliczeniową, ale o to, że po wielu obserwacjach testowych, sama specyfikacja (model teoretyczny) się zmienia. Rzecz w tym, że nawet przy ponad 500 danych dodanie jednej tylko obserwacji zmienia optymalny model. Tak więc, gdybym użył 501 danych już bym miał inne optimum, 502 jeszcze inne itd. Teoretycznie przy kolejnej reestymacji parametry dla zbędnych rzędów powinny zostać przez algorytm ustawione na zero, ale to z kolei rodzi pytanie o ustawienie maksymalnego rzędu. Przykładowo przy założeniu max 15 i przy 400 obserwacjach próby uczącej dostałem  ARMA(12, 15)-EGARCH-ARCH_In_Mean dla HQC, co znaczy, że albo max rzędów jest za mały, albo próba jest za mała. 

Aby zminimalizować ten problem, zwiększam próbę uczącą do 600 obserwacji (tygodniowe log-stopy zwrotu WIG20), a maksimum rzędów ARMA do 20. Natomiast wielkość próby testowej wyniesie tylko 15.  Poza tym kod pozostaje ten sam. W ten sposób wg min. HQC otrzymałem:

                    ARMA(9, 9)-ARCH_In_Mean
HQC                                 -4.4735
DAC (in sample)                      0.6133 

DAC (out of sample) 0.6000 

Ciekawe, że DAC in-sample i out-of-sample są w tym modelu takie same. Może przypadek, może nie. Wykonajmy diagnostykę:

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

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

Optimal Parameters
------------------------------------
        Estimate  Std. Error      t value Pr(>|t|)
mu     -0.000611    0.000000   -72428.506        0
ar1    -1.017510    0.000786    -1295.050        0
ar2     0.510817    0.000964      529.650        0
ar3     0.688893    0.000020    34654.649        0
ar4     0.042736    0.000037     1150.044        0
ar5     0.031574    0.000027     1166.599        0
ar6     0.802410    0.000001  1581829.476        0
ar7     0.981174    0.000001   839176.987        0
ar8    -0.423505    0.000006   -74666.301        0
ar9    -0.723909    0.000000 -1840101.437        0
ma1     1.014618    0.000022    46381.828        0
ma2    -0.521428    0.000221    -2362.691        0
ma3    -0.732535    0.000111    -6598.699        0
ma4    -0.121873    0.000055    -2210.317        0
ma5    -0.089846    0.000052    -1728.580        0
ma6    -0.894943    0.000205    -4362.136        0
ma7    -1.067543    0.000164    -6508.449        0
ma8     0.517891    0.000070     7402.304        0
ma9     0.849641    0.000082    10319.607        0
archm   0.000132    0.000000      315.194        0
omega  -1.188795    0.005020     -236.822        0
alpha1 -0.261897    0.000051    -5159.294        0
beta1   0.842014    0.000210     4006.494        0
gamma1  0.009565    0.000146       65.525        0
skew    0.741782    0.000413     1796.442        0

Robust Standard Errors:
        Estimate  Std. Error     t value Pr(>|t|)
mu     -0.000611    0.000000  -14098.715        0
ar1    -1.017510    0.000898   -1132.916        0
ar2     0.510817    0.001865     273.969        0
ar3     0.688893    0.000038   18278.881        0
ar4     0.042736    0.000151     282.800        0
ar5     0.031574    0.000111     283.606        0
ar6     0.802410    0.000003  307956.247        0
ar7     0.981174    0.000004  279684.711        0
ar8    -0.423505    0.000019  -22631.465        0
ar9    -0.723909    0.000001 -628047.744        0
ma1     1.014618    0.000123    8251.585        0
ma2    -0.521428    0.000587    -887.549        0
ma3    -0.732535    0.000303   -2414.904        0
ma4    -0.121873    0.000083   -1475.789        0
ma5    -0.089846    0.000103    -870.758        0
ma6    -0.894943    0.000507   -1766.114        0
ma7    -1.067543    0.000361   -2955.822        0
ma8     0.517891    0.000449    1154.102        0
ma9     0.849641    0.000402    2111.852        0
archm   0.000132    0.000002      80.308        0
omega  -1.188795    0.020893     -56.900        0
alpha1 -0.261897    0.000597    -438.494        0
beta1   0.842014    0.000153    5495.866        0
gamma1  0.009565    0.000228      41.865        0
skew    0.741782    0.000726    1021.332        0

LogLikelihood : 1388 

Information Criteria
------------------------------------
                    
Akaike       -4.5448
Bayes        -4.3616
Shibata      -4.5481
Hannan-Quinn -4.4735

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                      0.6302  0.4273
Lag[2*(p+q)+(p+q)-1][53]   21.5958  1.0000
Lag[4*(p+q)+(p+q)-1][89]   43.0403  0.6430
d.o.f=18
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic     p-value
Lag[1]                     0.6644 0.415008325
Lag[2*(p+q)+(p+q)-1][5]   22.9348 0.000002998
Lag[4*(p+q)+(p+q)-1][9]   25.8744 0.000005721
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.2817 0.500 2.000  0.5956
ARCH Lag[5]    0.3440 1.440 1.667  0.9283
ARCH Lag[7]    0.7047 2.315 1.543  0.9563

Nyblom stability test
------------------------------------
Joint Statistic:  no.parameters>20 (not available)
Individual Statistics:              
mu     0.01828
ar1    0.04389
ar2    0.04270
ar3    0.02659
ar4    0.01544
ar5    0.01972
ar6    0.02274
ar7    0.01703
ar8    0.01997
ar9    0.01801
ma1    0.03099
ma2    0.03226
ma3    0.03189
ma4    0.02808
ma5    0.02623
ma6    0.01821
ma7    0.06530
ma8    0.03776
ma9    0.01859
archm  0.01796
omega  0.15436
alpha1 0.19570
beta1  0.14572
gamma1 0.05599
skew   0.14442

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

Sign Bias Test
------------------------------------
                   t-value    prob sig
Sign Bias          1.72190 0.08561   *
Negative Sign Bias 1.37158 0.17071    
Positive Sign Bias 0.06065 0.95166    
Joint Effect       4.79282 0.18761    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     26.60      0.11432
2    30     38.30      0.11581
3    40     46.13      0.20113
4    50     62.50      0.09319

Słabo to wygląda. Test znaków i kwadratów reszt sugeruje, że model jest niewłaściwy.

I tak dochodzimy do drugiej kwestii - samej specyfikacji modelu. Jeszcze dla 300 obserwacji EGARCH(1, 1) był prawidłowy, ale 600 zmienia ten pogląd. Jest jeszcze wprawdzie możliwość, że zły rozkład, ale raczej nie o to chodzi.

Wszystko od początku

Całą analizę zacznijmy od początku. Na tych 600 danych treningowych i 15 testowych trzeba najpierw znaleźć najlepszy model ARMA (bez GARCH). Kod będzie bardzo podobny do omawianego poprzednio (wszystkie pakiety i zmienne te same), ale bez ARCH-In-Mean i użyć funkcji arfimaspec(), arfimafit(), arfimaforecast() zamiast ugarchspec(), ugarchfit(), ugarchforecast():

rozklad <- "snorm"

hqc <- list()

dacDopas <- list()

dacPrognoza <- list()

wyrzuc <- 15

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 bez k (ARCH-in-mean)

siatka_parametrow <- expand.grid(i = 0:20, j = 0:20)

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

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

  i <- siatka_parametrow$i[param]

  j <- siatka_parametrow$j[param]

  rzad_modelu <- sprintf("ARMA(%d, %d)", i, j)

  tryCatch({

    dopasArma <- NULL

    mojArma <- arfimaspec(mean.model = list(armaOrder = c(i, j)), 

                          distribution.model = rozklad)

    dopasArma <- arfimafit(spec = mojArma, data = stopaZwrotu, out.sample = wyrzuc,

                           solver = "hybrid")

    stopaZwrotuDopas <- fitted(dopasArma)

    stopaZwrotuZgodna <- stopaZwrotuDopas + residuals(dopasArma)

    prognozaArma <- arfimaforecast(dopasArma, 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)) {

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

                                  actual = as.numeric(stopaZwrotuZgodna), 

                                  test = "PT")$DirAcc

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

    }

    return(list(

      rzad_modelu = rzad_modelu,

      hqc = infocriteria(dopasArma)[4,],

      dacDopas = dacDopas_wartosc,

      dacPrognoza = dacPrognoza_wartosc

    ))

  }, error = function(e) {

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

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

  })

}

# 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

    }

  }

}

# 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

Dostałem:

> kryterium_hqc
                    ARMA(7, 11)
HQC                     -4.3222
DAC (in sample)          0.6017
DAC (out of sample)      0.8667
> kryterium_dacDopas
                    ARMA(16, 11)
HQC                      -4.2633
DAC (in sample)           0.6283
DAC (out of sample)       0.6000
> kryterium_dacPrognoza
                    ARMA(6, 4) ARMA(6, 18)
HQC                    -4.2574      -4.264
DAC (in sample)         0.5417       0.570
DAC (out of sample)     1.0000       1.000

Dalej wyciągamy optymalne rzędy:

# 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)

# parametry - DAC out of sample

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

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

Wykonajmy diagnostykę, analogicznie jak poprzednio wg HQC:

dopasArma_hqc <- arfimafit(spec = mojArma_hqc, data = stopaZwrotu, out.sample = wyrzuc, solver="hybrid")

dopasArma_hqc

> dopasArma_hqc

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(7,0,11)
Distribution	: snorm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error     t value Pr(>|t|)
mu    -0.000163    0.000001     -118.88        0
ar1    1.708939    0.000055    31171.71        0
ar2   -1.589454    0.000048   -33434.75        0
ar3    0.941900    0.000042    22383.28        0
ar4   -1.292676    0.000001 -1403760.91        0
ar5    1.688498    0.000060    28287.27        0
ar6   -1.389024    0.000042   -32791.53        0
ar7    0.337973    0.000043     7801.90        0
ma1   -1.803387    0.000060   -29855.67        0
ma2    1.876807    0.000059    31845.44        0
ma3   -1.393470    0.000044   -31363.61        0
ma4    1.874833    0.000068    27413.36        0
ma5   -2.382959    0.000069   -34491.56        0
ma6    2.152583    0.000074    29113.30        0
ma7   -1.013392    0.000043   -23773.12        0
ma8    0.566727    0.000027    20790.58        0
ma9   -0.357758    0.000017   -20796.19        0
ma10   0.187237    0.000010    18174.45        0
ma11  -0.066118    0.000014    -4891.35        0
sigma  0.026299    0.000060      439.52        0
skew   0.803218    0.001613      498.00        0

Robust Standard Errors:
       Estimate  Std. Error    t value Pr(>|t|)
mu    -0.000163    0.000001    -245.88        0
ar1    1.708939    0.000100   17094.56        0
ar2   -1.589454    0.000060  -26331.39        0
ar3    0.941900    0.000053   17769.70        0
ar4   -1.292676    0.000004 -321135.66        0
ar5    1.688498    0.000068   24959.04        0
ar6   -1.389024    0.000018  -75620.19        0
ar7    0.337973    0.000035    9540.15        0
ma1   -1.803387    0.000209   -8649.06        0
ma2    1.876807    0.000157   11965.81        0
ma3   -1.393470    0.000031  -45389.79        0
ma4    1.874833    0.000165   11333.48        0
ma5   -2.382959    0.000120  -19852.77        0
ma6    2.152583    0.000209   10290.45        0
ma7   -1.013392    0.000129   -7881.81        0
ma8    0.566727    0.000037   15484.79        0
ma9   -0.357758    0.000041   -8716.28        0
ma10   0.187237    0.000009   19894.66        0
ma11  -0.066118    0.000021   -3178.62        0
sigma  0.026299    0.000071     371.02        0
skew   0.803218    0.000424    1892.55        0

LogLikelihood : 1336 

Information Criteria
------------------------------------
                    
Akaike       -4.3821
Bayes        -4.2282
Shibata      -4.3844
Hannan-Quinn -4.3222

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                       1.659 0.19768
Lag[2*(p+q)+(p+q)-1][53]    35.169 0.00000
Lag[4*(p+q)+(p+q)-1][89]    54.654 0.03086

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic               p-value
Lag[1]                      1.762 0.1844147994362065335
Lag[2*(p+q)+(p+q)-1][2]    35.427 0.0000000004027803646
Lag[4*(p+q)+(p+q)-1][5]    57.765 0.0000000000000003331


ARCH LM Tests
------------------------------------
             Statistic DoF             P-Value
ARCH Lag[2]      67.45   2 0.00000000000000222
ARCH Lag[5]      67.76   5 0.00000000000029932
ARCH Lag[10]     70.05  10 0.00000000004330802

Nyblom stability test
------------------------------------
Joint Statistic:  no.parameters>20 (not available)
Individual Statistics:             
mu    0.05184
ar1   0.08997
ar2   0.06947
ar3   0.05230
ar4   0.02741
ar5   0.04652
ar6   0.07092
ar7   0.08179
ma1   0.09581
ma2   0.07342
ma3   0.05758
ma4   0.03403
ma5   0.03789
ma6   0.04300
ma7   0.07302
ma8   0.10824
ma9   0.07579
ma10  0.04847
ma11  0.03839
sigma 0.66621
skew  0.02743

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

I co się okazuje? Że występuje jednak efekt ARCH. Sama ARMA jest niedopuszczalna. 

Sprawdźmy jeszcze model wg kryterium DAC w próbie testowej:

> dopasArma_DAC

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(6,0,18)
Distribution	: snorm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error   t value Pr(>|t|)
mu    -0.000484    0.000003   -181.01        0
ar1    1.652690    0.000052  31715.16        0
ar2   -0.569904    0.000021 -26571.17        0
ar3   -0.691440    0.000024 -29358.26        0
ar4   -0.370370    0.000016 -23705.84        0
ar5    1.335052    0.000028  46954.51        0
ar6   -0.810241    0.000022 -36509.08        0
ma1   -1.807455    0.000065 -27683.56        0
ma2    0.896344    0.000043  20628.63        0
ma3    0.479002    0.000017  28772.19        0
ma4    0.307185    0.000013  23760.96        0
ma5   -1.300847    0.000056 -23158.98        0
ma6    0.999150    0.000045  22099.33        0
ma7   -0.255216    0.000013 -19665.03        0
ma8    0.146197    0.000009  15602.83        0
ma9    0.010041    0.000005   1948.42        0
ma10  -0.192279    0.000009 -20638.50        0
ma11   0.181860    0.000009  19169.60        0
ma12  -0.126310    0.000007 -19332.60        0
ma13   0.129111    0.000008  15919.81        0
ma14  -0.135228    0.000010 -13197.50        0
ma15   0.080054    0.000005  16270.07        0
ma16  -0.257084    0.000014 -19028.90        0
ma17   0.334883    0.000013  24885.88        0
ma18  -0.170361    0.000002 -78823.90        0
sigma  0.026828    0.000039    684.76        0
skew   0.736525    0.000896    822.24        0

Robust Standard Errors:
       Estimate  Std. Error    t value Pr(>|t|)
mu    -0.000484    0.000005    -88.585        0
ar1    1.652690    0.000225   7358.377        0
ar2   -0.569904    0.000030 -18787.065        0
ar3   -0.691440    0.000020 -33985.665        0
ar4   -0.370370    0.000007 -53955.284        0
ar5    1.335052    0.000061  22010.081        0
ar6   -0.810241    0.000078 -10327.488        0
ma1   -1.807455    0.000418  -4325.866        0
ma2    0.896344    0.000228   3928.523        0
ma3    0.479002    0.000019  25779.508        0
ma4    0.307185    0.000031  10005.522        0
ma5   -1.300847    0.000312  -4166.365        0
ma6    0.999150    0.000079  12603.512        0
ma7   -0.255216    0.000037  -6950.382        0
ma8    0.146197    0.000019   7556.457        0
ma9    0.010041    0.000026    388.235        0
ma10  -0.192279    0.000006 -33015.132        0
ma11   0.181860    0.000014  12830.990        0
ma12  -0.126310    0.000003 -39589.998        0
ma13   0.129111    0.000016   8046.554        0
ma14  -0.135228    0.000022  -6183.826        0
ma15   0.080054    0.000032   2482.382        0
ma16  -0.257084    0.000028  -9035.050        0
ma17   0.334883    0.000027  12423.452        0
ma18  -0.170361    0.000007 -26122.951        0
sigma  0.026828    0.000057    474.242        0
skew   0.736525    0.003314    222.255        0

LogLikelihood : 1329 

Information Criteria
------------------------------------
                    
Akaike       -4.3406
Bayes        -4.1427
Shibata      -4.3444
Hannan-Quinn -4.2636

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                          statistic p-value
Lag[1]                        5.469 0.01936
Lag[2*(p+q)+(p+q)-1][71]     35.802 0.63160
Lag[4*(p+q)+(p+q)-1][119]    62.706 0.30838

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic           p-value
Lag[1]                     0.6252 0.429113495341529
Lag[2*(p+q)+(p+q)-1][2]   28.6335 0.000000024213054
Lag[4*(p+q)+(p+q)-1][5]   45.7761 0.000000000001031


ARCH LM Tests
------------------------------------
             Statistic DoF           P-Value
ARCH Lag[2]      55.77   2 0.000000000000775
ARCH Lag[5]      57.73   5 0.000000000035680
ARCH Lag[10]     58.90  10 0.000000005842875

Nyblom stability test
------------------------------------
Joint Statistic:  no.parameters>20 (not available)
Individual Statistics:             
mu    0.07987
ar1   0.12280
ar2   0.10719
ar3   0.05115
ar4   0.04223
ar5   0.07941
ar6   0.11593
ma1   0.11203
ma2   0.06554
ma3   0.03642
ma4   0.07160
ma5   0.11188
ma6   0.12608
ma7   0.10324
ma8   0.05420
ma9   0.03809
ma10  0.07410
ma11  0.11379
ma12  0.12568
ma13  0.09919
ma14  0.05052
ma15  0.03972
ma16  0.07790
ma17  0.11582
ma18  0.12465
sigma 0.63364
skew  0.06743

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

To samo - mamy tu do czynienia z efektami ARCH. Aby model poprawnie działał, musimy je uwzględnić.

Wiemy już na pewno, że sama ARMA jest błędna. Dochodzimy do wniosku, że musi być to jakiś model GARCH, ale nie EGARCH. Sprawdziłem więc standardowy GARCH(1,1) oznaczany "sGARCH" - ale także okazał się niewłaściwy. Trzeba szukać dalej. Wiadomo tylko, że rzędy części GARCH (1, 1) są prawidłowe, bo gdyby nie były, to w resztach pojawiłyby się efekty ARCH nie uwzględnione przez model - a ich nie było w EGARCH (ani w sGARCH jak sprawdzałem). 

Przychodzą mi dwie możliwości, które trzeba po kolei sprawdzić. Po pierwsze może być to inny model, jak TGARCH, APARCH, GJR-GARCH, NGARCH, NAGARCH, iGARCH i csGARCH. Po drugie EGARCH może być prawidłowy, ale może występować jakaś niestabilność czy zmiana struktury wariancji, której test Nybloma nie uwzględnił. Stąd może być potrzebny zupełnie inny model albo bardziej zaawansowany GARCH, jak np. TV-GARCH (time-varying GARCH), którego nie ma obecnie w rugarch.

Jak dotąd sprawdziłem we wszystkich konfiguracjach sGARCH, TGARCH, GJR-GARCH, IGARCH i APARCH. Mimo użycia przetwarzania równoległego poszukiwanie optimum dla każdego zajmuje wiele godzin. Okazało się, że żaden model wcale nie był lepszy od EGARCH. Sprawdziłem dodatkowo - ale już bez ARMA - wiele modeli z GARCH innych rzędów, czyli GARCH(1, 2), GARCH(2, 1), GARCH(2, 2). Sprawdzałem pod tym kątem też NGARCH i NAGARCH. Nic to nie dało.

Jedną z potencjalnych przyczyn trudności modelowania warunkowej wariancji, jest zmienność jej parametrów w czasie. Test Nybloma w rugarch nie wykryje tej zmienności, jeśli zmiana nie będzie systematyczna, ale za to będzie systematycznie zanikać w czasie. Moglibyśmy użyć innych testów. Dotychczas sprawdzałem zmienność wariancji przy pomocy testu Goldfelda-Quandta, Breuscha-Pagana, Harrisona-McCabe (zob. tu) oraz CUSUM (tu i tu). Niektóre z tych testów mogą jednak pomylić warunkową heteroskedastyczność z niestabilnością wariancji niewarunkowej. 

TV-GARCH

Rozwiązaniem problemu jest pakiet tvgarch, służący do symulacji, estymacji i analizy procesu TV-GARCH (Time Varying GARCH). Pakiet zawiera funkcję tvgarchTest(), która testuje proces TV-GARCH(1,1). Instalujemy i ładujemy pakiet:

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

  install.packages("tvgarch")

}

library("tvgarch")

Najpierw przeprowadzimy kilka symulacji sGARCH i eGarch z rozkładem skośno-normalnym i wykonamy test:

set.seed(30)

sym <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 1)), 

                   distribution.model = "snorm", fixed.pars = list(omega = 0.5, 

                                                                   alpha1 = 0.6,

                                                                   beta1 = 0.3,

                                                                   mu = 0.01,

                                                                   ar1 = 0.6,

                                                                   ma1 = 0.2,

                                                                   skew = 0.7))

symGarch <- ugarchpath(sym, n.sim = 600)

symGarch <- symGarch@path$seriesSim

ts.plot(symGarch)

I włączamy funkcję:

tvGarchTest <- tvgarchTest(y = as.numeric(symGarch))

tvGarchTest

Efekt:

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.5632 0.8630 0.28048
Std. Error:      0.1232 0.1134 0.04657
                     
Log-likelihood: -1309


Transition variable in TV-GARCH(1,1): time 

Results from the Non-Robust TR^2 Test: 

                 NonRobTR2 p-value
H0:B3=B2=B1=0       3.7317  0.2919
H03:B3=0            0.7273  0.3937
H02:B2=0|B3=0       2.2520  0.1334
H01:B1=0|B3=B2=0    0.7589  0.3837

Results from the Robust TR^2 Test: 

                 RobTR2 p-value
H0:B3=B2=B1=0     6.268  0.0993
H03:B3=0          1.393  0.2378
H02:B2=0|B3=0     2.195  0.1385
H01:B1=0|B3=B2=0  2.694  0.1007

                                Single
No. of locations (alpha = 0.05)      0

Są tu dwie wersje testu: nieodporny (non-robust) i odporny (robust), podobnie jak to było w rugarch. Odporność należy tu rozumieć jako odporność na odchylenia od rozkładu normalnego. Ponieważ nasze reszty modelu nie mają rozkładu normalnego, to ten odporny będzie ważniejszy - więc tylko na niego patrzymy. Dzięki temu odróżnieniu możemy ocenić czy błędy odchyleń zachowują się gaussowsko. Parametr B1 można interpretować jako średnią warunkową wariancji niewarunkowej. Parametr B2 - jako wariancję warunkową wariancji niewarunkowej. Parametr B3 - jako skośność warunkową wariancji niewarunkowej. Najważniejsza jest pierwsza hipoteza. Ponieważ wskazuje p-value > 0,05 , to znaczy, że nie ma podstaw do odrzucenia hipotezy stałego GARCH. Tak też wskazuje ostatnia statystyka  No. of locations = 0 liczbę punktów przejść z jednego modelu do drugiego.

Dla set.seed(1100)

Dostałem podobne statystyki. 

Teraz symulacja EGARCH(1, 1):

sym <- ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 1)), 

                   distribution.model = "snorm", fixed.pars = list(omega = 0.5, 

                                                                   alpha1 = 0.6,

                                                                   beta1 = 0.3,

                                                                   mu = 0.01,

                                                                   ar1 = 0.6,

                                                                   ma1 = 0.2,

                                                                   gamma1 = 0.8,

                                                                   skew = 0.7))

set.seed(1100)

symGarch <- ugarchpath(sym, n.sim = 600)

symGarch <- symGarch@path$seriesSim

ts.plot(symGarch)

Tak samo test odrzucił tu TV-GARCH.

To teraz nasze (logarytmiczne) stopy zwrotu WIG20

plot(stopaZwrotu)

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.0001704 0.15324 0.63531
Std. Error:   0.0003021 0.09059 0.05965
                    
Log-likelihood: 1346


Transition variable in TV-GARCH(1,1): time 

Results from the Non-Robust TR^2 Test: 

                 NonRobTR2 p-value
H0:B3=B2=B1=0       5.7505  0.1244
H03:B3=0            3.7567  0.0526
H02:B2=0|B3=0       0.0233  0.8787
H01:B1=0|B3=B2=0    1.9829  0.1591

Results from the Robust TR^2 Test: 

                  RobTR2 p-value
H0:B3=B2=B1=0    10.2591  0.0165
H03:B3=0          3.5790  0.0585
H02:B2=0|B3=0     0.0203  0.8867
H01:B1=0|B3=B2=0  6.4987  0.0108

                                Single
No. of locations (alpha = 0.05)      1

Test rzeczywiście odrzuca stabilność wariancji niewarunkowej. Towarzyszy temu niskie p-value. Jednocześnie wyłapał 1 punkt zmiany struktury modelu wariancji. Dostajemy dowód, że albo należy posługiwać się mniejszą próbą, albo właśnie TV-GARCH. 

Chociaż jest jeszcze jedna możliwość - nawet bardziej naturalna. Można użyć nieparametrycznego GARCH, analogicznie do regresji nieparametrycznej. Skoro to takie naturalne, to po co w ogóle parametryczne modele? Głównie z powodu dostępności i tak jak w przypadku rugarch z całą paletą diagnostyczną. Drugim powodem jest to, że powinny być one zawsze punktem wyjścia - bo przecież mogą się okazać lepsze od wersji nieparametrycznej. Poza tym używanie skomplikowanych metod bez zrozumienia podstaw jest bez sensu.