Giełdowy Racjonalista: maja 2025

Jeśli powiem, że każda prognoza gospodarcza czy giełdowa jest okraszona niepewnością, to będzie to banał. Ale jeśli powiem, że celem prognozy powinno być nie tyle odgadnięcie przyszłej wartości, co minimalizacja jej niepewności, to już banałem nie jest. Dla uproszczenia założymy, że niepewność to przedział ufności dla prognozy naszego modelu. Im więcej informacji zawiera model, tym przedział ufności staje się mniejszy, zbliżając się tym samym do prognozy punktowej. Widzimy więc, że idea nie jest wcale banalna: należy zebrać jak najwięcej istotnych informacji, aby zawężać przedział prognozy, zmniejszając tym samym ryzyko.

Analogicznie do przedziału ufności, Hansen et al.* zaproponowali zbiór ufności dla wielu modeli: z całej puli testowanych modeli, poszukiwane są te najlepsze, które wpadają do zbioru ufności. Im więcej posiadamy informacji, tym mniejszy będzie ten zbiór, zbliżając się tym samym do jednego najlepszego modelu.

Standardową metodą wyboru najlepszego modelu jest użycie kryterium informacyjnego, np. AIC, BIC czy HQIC. Niestety nie uwzględniają one niepewności czy przypadkowości parametrów. Z tym się ciągle zmagamy. Przykładowo, mam zawsze wątpliwości co do używania modeli parametrycznych, ponieważ wydają się one sztuczne, realizując uśrednianie na siłę. Dlatego tak ważne są testy stabilności parametrów. Ale stabilność to dopiero początek, bo:

(1) jeśli spośród wielu modeli różnica w kryterium informacyjnym jest niewielka, to może ona wynikać po prostu z przypadku,

(2) dodanie tylko jednej obserwacji bardzo często zmienia wybór modelu.

W rugarch jest funkcja mcsTest (Model Confidence Set - MCS), która testuje, czy dany model należy zachować w zbiorze ufności czy go odrzucić. Ma ona następujące argumenty:

- losses: macierz lub ramka danych z wartościami strat dla każdego modelu (każda kolumna odpowiada jednemu modelowi). To w jaki sposób definiujemy straty, zależy od nas. W przypadku rozkładów symetrycznych, najczęściej stosuje się kwadraty odchyleń od prognozy lub wartości absolutne tych odchyleń. Dla rozkładów niesymetrycznych stosować można semi-odchylenia albo miary oparte na wartości zagrożonej.

- alpha: poziom istotności (domyślnie 0.05).

- nboot: liczba replikacji bootstrapowych (domyślnie 100). Bootstrap to procedura wielokrotnego losowania próbki ze zwracaniem. nboot to właśnie liczba tych losowań.

- nblock: długość bloku dla bootstrappu blokowego (domyślnie 1). Jeżeli obserwacje mogą być zależne od siebie, to bootstrap musi losować blok występujących po sobie x obserwacji. Jeżeli uważamy, że każda obserwacja jest niezależna od innej, to x = 1. Im rząd autokorelacji większy, tym blok powinien być naturalnie też większy.

- boot: rodzaj bootstrapu ("stationary" lub "block"). Chociaż w obydwu przypadkach mamy do czynienia z blokami, to metoda blokowa, "block", oznacza, że całą próbę dzielimy na stałą liczbę przedziałów, tak że każdy przedział składa się bloku obserwacji o długości wyznaczonej dokładnie przez nblock. Potem dopiero losujemy blok, powtarzamy to losowanie dokładnie nboot razy. Natomiast metoda stacjonarna, "stationary", tworzy bloki o losowej długości ze średnią równą nblock. Jeśli występuje autokorelacja w szeregu czasowym, to blok nie może być za krótki. Najczęściej "stationary" jest lepszy od "block".

Co właściwie robi ten test? W skrócie, kiedy mamy macierz strat dla każdego modelu, liczymy różnice między każdym modelem (prognozą każdego modelu) i są one odpowiednio uśredniane. Tworzymy bloki bootstrapowe i dla każdego bloku obliczamy średnie różnice między średnią stratą w bloku, średnią stratą dla całego modelu i średnią stratą ze wszystkich modeli. Dla każdego modelu będą różne wartości tworzące rozkład różnic. Wartości zbyt odchylone na plus (tzn. prawostronne ogony rozkładu) zostają odrzucone i tak uzyskujemy zbiór ufności.

Oczywiście najpierw sami budujemy funkcję strat. Ja użyłem definicji najbardziej dla mnie intuicyjnej: jeżeli znak prognozy różni się od znaku wartości rzeczywistej, to bierzemy moduł z różnicy między nimi. Jeżeli znak jest ten sam, to strata wynosi zero.

Przetestujmy tę funkcję. Będziemy używać mcsTest na zbiorze uczącym. Reguły używania poniższych funkcji w rugarch omówiłem szczegółowo tutaj.

Stwierdziłem wcześniej, że próba między 300 a 400 tygodniowych danych giełdowych dla GARCH wydaje się już w miarę odpowiednia. Zrobię symulację EGARCH o próbie równej 330. Powiedzmy, że prawdziwy proces to EGARCH(1, 1)-ARMA(3, 1). Obok niego mamy też fałszywe - wszystkie EGARCH: ARMA(0, 0), ARMA(1, 0), ARMA(2, 0), ARMA(3, 0). Zobaczmy które z nich zostaną wykluczone przez mcsTest, a które przyjęte do zbioru ufności. Zrobimy test w dwóch wariantach: blokowy i stacjonarny. Ustawienia: nblock = 20 oraz nboot = 1000.

library("rugarch")

library("xts")

# 1. Definicje stałych parametrów

stale_3_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

ma1 = 0.4,

gamma1 = 0.05

)

stale_0_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

gamma1 = 0.05

)

stale_1_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

gamma1 = 0.05

)

stale_2_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

gamma1 = 0.05

)

stale_3_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

gamma1 = 0.05

)

stale_0_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ma1 = 0.4,

gamma1 = 0.05

)

# 2. Specyfikacje modeli EGARCH(1,1) z różnymi specyfikacjami średniej:

spec_3_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,1), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_3_1

)

spec_0_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_0_0

)

spec_1_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(1,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_1_0

)

spec_2_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(2,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_2_0

)

spec_3_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_3_0

)

spec_0_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,1), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_0_1

)

# 3. Symulacja 330 obserwacji według poprawnego modelu (ARMA(3,1))

set.seed(765)

symulacja <- ugarchpath(spec_3_1, n.sim = 330)

dane_symulowane <- as.numeric(fitted(symulacja))

# 4. Dopasowanie modeli do całego szeregu (in-sample)

dopasowanie_3_1 <- ugarchfit(spec = spec_3_1, data = dane_symulowane)

dopasowanie_0_0 <- ugarchfit(spec = spec_0_0, data = dane_symulowane)

dopasowanie_1_0 <- ugarchfit(spec = spec_1_0, data = dane_symulowane)

dopasowanie_2_0 <- ugarchfit(spec = spec_2_0, data = dane_symulowane)

dopasowanie_3_0 <- ugarchfit(spec = spec_3_0, data = dane_symulowane)

dopasowanie_0_1 <- ugarchfit(spec = spec_0_1, data = dane_symulowane)

# 5. Zapisanie wyników dopasowania w liście (dla przejrzystości)

dopasowania <- list(

model_3_1 = dopasowanie_3_1,

model_0_0 = dopasowanie_0_0,

model_1_0 = dopasowanie_1_0,

model_2_0 = dopasowanie_2_0,

model_3_0 = dopasowanie_3_0,

model_3_0 = dopasowanie_0_1

)

print("Dopasowane modele:")

print(names(dopasowania))

# 6. Definicja funkcji strat:

Funkcja_strat <- function(rzeczywiste, prognoza) {

roznica <- rzeczywiste - prognoza

strata <- ifelse(sign(rzeczywiste) != sign(prognoza), abs(roznica), 0)

return(strata)

}

# 7. Obliczenie funkcji strat dla każdego modelu (in-sample)

strata_3_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_1))

strata_0_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_0))

strata_1_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_1_0))

strata_2_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_2_0))

strata_3_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_0))

strata_0_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_1))

# Łączenie strat w macierz – wiersze odpowiadają obserwacjom, kolumny modelom

straty <- cbind(

model_3_1 = strata_3_1,

model_0_0 = strata_0_0,

model_1_0 = strata_1_0,

model_2_0 = strata_2_0,

model_3_0 = strata_3_0,

model_3_0 = strata_0_1

)

# 8. Usunięcie wierszy składających się tylko z zer

straty <- straty[rowSums(straty != 0) > 0, ]

straty <- coredata(straty)

# 9. Test MCS

# 9.1. Wariant 1: stacjonarny

mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "stationary")

# 9.2. Wariant 1: blokowy

mcs_block <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "block")

# 10. Wyświetlenie wyników:

cat("\nWyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):\n")

print(mcs_stationary)

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 1 5

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.004
[5,] 0.697
[6,] 1.000

$excludedR
[1] 6 2 4 3

$includedSQ
[1] 1 5

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.003
[5,] 0.697
[6,] 1.000

$excludedSQ
[1] 6 2 4 3

Należy zauważyć, że mamy tu 2 typy testów: statystyka R (Range statistic) oraz SQ (Semi-Quadratic statistic). Pierwszy z nich skupia się na różnicach w maksymalnych stratach między modelami, a drugi na różnicach w kwadratach strat. Najlepiej, jeśli obydwa wskazują ten sam zbiór, bo wtedy dostajemy potwierdzenie, że zbiór ufności zawiera odpowiednie modele.

Jak widać obydwa testy rzeczywiście wskazują to samo: included 1 oraz 5, czyli zbiór ufności zawiera tylko dwa modele - pierwszy oraz piąty. Pozostałe modele zostają odrzucone (excluded), zarówno przez R jak SQ. Możemy powiedzieć, że test zadziałał poprawnie, mimo przyjęcia błędnego modelu ARMA(3, 0). Model ten jest podobny do prawidłowego ARMA(3, 1), a co więcej poprawny model ma największe p-value = 1 (piąty 0.697). O tym może jeszcze powiem. Otóż wydrukowana lista może być trochę myląca, bo kolejność nie pokazuje p-value w kolejności od modelu 1 do 6, tylko w kolejności od najbardziej istotnego do najmniej istotnego wyniku dotyczącego odrzucenia modelu ze zbioru ufności. Im mniejsze p-value, tym większe prawdopodobieństwo, że dany model należy wyrzucić ze zbioru.

Żeby zawęzić zbiór ufności, należy zwiększyć próbę. Po zwiększeniu jej do 1600, dostałem:

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 1

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.047
[6,] 1.000

$excludedR
[1] 6 2 4 3 5

$includedSQ
[1] 1

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.047
[6,] 1.000

$excludedSQ
[1] 6 2 4 3 5

Przy tak dużej próbie test nie ma już wątpliwości, że jedyny prawidłowy model to ten pierwszy.

Bardzo podobne wyniki dostałem dla boot = "block" i nie ma sensu pokazywać ich. Podobnie jak dla "stationary" musiałem ustawić próbę liczną 1600, aby uzyskać wykluczenie wszystkich fałszywych modeli. Co jednak zauważyłem, to to, że w przypadku "block" wykluczenie modelu 5 było ledwo, ledwo, na styk, bo p-value wyniosło 0.049 (dla stationary też, ale niższe, bo 0.047). Stacjonarny test generalnei jest lepszą alternatywą.

Powyższe przykłady zakładały rozkład normalny reszt modelu. Dla rynków może być lepszy inny, np. skośny rozkład t-Studenta. W takiej sytuacji potrzebne są dwa dodatkowe parametry do stałych: skew oraz shape. Zamiast "norm" będzie "sstd". Przykładowo kod będzie taki:

# 1. Definicje stałych parametrów

stale_3_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

ma1 = 0.4,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_0_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_1_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_2_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_3_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_0_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ma1 = 0.4,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

# 2. Specyfikacje modeli EGARCH(1,1) z różnymi specyfikacjami średniej:

spec_3_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,1), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_3_1

)

spec_0_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_0_0

)

spec_1_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(1,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_1_0

)

spec_2_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(2,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_2_0

)

spec_3_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_3_0

)

spec_0_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,1), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_0_1

)

# 3. Symulacja 330 obserwacji według poprawnego modelu (ARMA(3,1))

set.seed(55)