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)
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.1. Wariant stacjonarny
mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "stationary")
# 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] 5 1
$pvalsR
[,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.071
[6,] 1.000
$excludedR
[1] 6 2 4 3
$includedSQ
[1] 5 1
$pvalsSQ
[,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.071
[6,] 1.000
$excludedSQ
[1] 6 2 4 3Jak widać, ponownie można by podwyższyć p-value do 10%, żeby dostać tylko model 1. Tak czy inaczej test działa bardzo dobrze także dla innych rozkładów.
Oczywiście to bardzo wzorcowe przykłady i w praktyce będzie to trudniejsze, ale to jest dobry punkt wyjścia.
* Hansen, P., A. Lunde, and J. Nason (2011). The model confidence set. Econometrica 79,
453–497.
