środa, 3 lipca 2024

Modelowanie ryzyka z rugarch (język R)

Ryzyko można szacować i modelować na różne sposoby. Standardowo obliczamy wariancję stopy zwrotu i będzie to nawet poprawne, jeśli się ona nie zmienia. Pytanie tylko dlaczego uważamy, że powinna być stała? Co więcej, nawet gdy średnia wariancja się nie zmienia, to nie znaczy, że nie można znaleźć dokładniejszej miary ryzyka - taką jest wariancja warunkowa, która stanowi istotę procesu GARCH

Pakietu rugarch w R już używałem w tym wpisie , ale tym razem skupię się tylko na nim i zrobię przykład na WIG20. Ładujemy pakiet:

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

  install.packages("rugarch")

}

  library("rugarch")

Dodatkowo załaduję "xts", który przyda się dla prawdziwych danych oraz "forecast" do symulacji ARMA.

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

  install.packages("xts")

}

  library("xts")

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

  install.packages("forecast")

}

library("forecast")


Specyfikacja modelu

Stwórzmy sztuczny proces GARCH. Pierwsza sprawa to budowa teoretycznego modelu - specyfikacja. To się robi przy pomocy funkcji ugarchspec(). Na specyfikację składają się głównie 4 elementy:

1) model (funkcja) wariancji warunkowej: variance.model

2) model (funkcja) średniej warunkowej: mean.model

3) rozkład warunkowy zmian (np. stóp zwrotu, nie warunkowej wariancji): distribution.model

4) wartości parametrów (przy symulacji procesu należy podać wszystkie parametry modelu): fixed.pars


Ad 1) variance.model jako lista składa się głównie z dwóch części: rodzaju modelu (model) oraz rzędów (garchOrder). Z rodzajów modelu do wyboru mamy:

“sGARCH” - standardowy GARCH,

“eGARCH” - wykładniczy GARCH (The exponential GARCH)

“gjrGARCH” -  Glosten-Jagannathan-Runkle GARCH

“apARCH” - The Asymmetric Power ARCH

“iGARCH” - zintegrowany GARCH (The integrated GARCH)

 “csGARCH” - The Component sGARCH

“fGARCH” - full GARCH - rodzina modeli podzielona na podmodele: 

a) “GARCH” - znów sGARCH

b) “TGARCH” - progowy GARCH (The Threshold GARCH)

c) “AVGARCH” - The Absolute Value GARCH

 d) “NGARCH” - nieliniowy GARCH

e) “NAGARCH” - nieliniowy asymetryczny GARCH

f) “APARCH” -  znowu APARCH (nie wiem dlaczego, ale to nie jest dokładnie to samo)

g) “GJRGARCH”- znowu gjrGARCH (jak wyżej)


Drugim elementem tej części są rzędy, czyli garchOrder. Należy zwrócić uwagę, że jest to dokładna analogia do ARMA(p, q), z tym że kolejność jest odwrotna: najpierw było ARCH(q), czyli średnia krocząca wariancji warunkowej, a potem dołożono część autoregresyjną tej wariancji, w sumie GARCH(q, p). Standardowo i na początek zawsze używamy q = 1, p = 1. Rzędy będą wektorem, którego elementami będzie liczba rzędów przyporządkowana danemu parametrowi. Każdy model ma swoje parametry, tak że wybierając odpowiedni typ modelu należy mieć to na względzie przy zapisywaniu tego wektora. Czyli ARCH(1) zapiszemy jako garchOrder = c(1, 0), a GARCH(1, 1) jako  garchOrder = c(1, 1). 

Przykład 1: ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 2)))

Przykład 2, dla podmodeli: ugarchspec(variance.model = list(model = "fGARCH", submodel = "NGARCH", garchOrder = c(1, 1)))


Ad 2) mean.model  jako lista zawiera dwie podstawowe części: armaOrder oraz external.regressors. Jest to po prostu ARMAX dla zmiennej, jak stopa zwrotu. Pierwszą z wymienionych części stosujemy analogicznie jak dla garchOrder. Druga część to zmienne egzogeniczne, czyli zewnętrzne regresory ARMAX. Specjalnym przypadkiem ARMAX jest GARCH-in-Mean, w którym za zmienną egzogeniczną wstawiamy samą wariancję warunkową (lub jej pierwiastek). Jest to ciekawy przykład sprzężenia zwrotnego między średnią warunkową a jej wariancją, np. między oczekiwaną (warunkową) stopą zwrotu a jej ryzykiem. Nie musimy jednak sami wyciągać tej wariancji i przyłączać do modelu, bo taki model jest do dyspozycji tutaj. Po prostu zamiast external.regressors stosujemy dodatkowy argument archm = TRUE. 

Przykład1: ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 1)))

Przykład2: ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 0), archm = TRUE))

Warto dodać, że w tej części także można ustalić czy średnia ma się zachowywać zgodnie z procesem ARFIMA.


Ad 3) distribution.model, tj. rozkład prawdopodobieństwa naszej zmiennej może to być:

“norm” - rozkład normalny

“snorm” - skośno-normalny (skew-normal distribution)

“std” - t-Studenta

“sstd” skośny t-studenta (skew-student), 

“ged” - uogólniony rozkład błędów (the generalized error distribution), 

“sged” - skośny uogólniony rozkład błędów (the skew-generalized error distribution), 

“nig” - normalny odwrotny rozkład gaussowski (the normal inverse gaussian distribution), 

“ghyp” - uogólniony hiperboliczny (the Generalized Hyperbolic),

“jsu” - Johnson’s SU distribution (SU to nie skrót, ale po prostu nazwa zmiennej w tym rozkładzie, podobnie jak litera "t" w t-Studenta).

Przykład: ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 2)), distribution.model = "norm")


Ad 4) fixed.pars używamy w zasadzie tylko do symulacji. Ustalamy wartości parametrów dla wybranego modelu. Nazewnictwo parametrów:

alpha1, alpha2, ... dla ARCH(q)

beta1, beta2, ..., dla GARCH(p)

omega - wyraz wolny funkcji wariancji

gamma1, gamma2, ...  dla asymetrii, tj. APARCH, EGARCH i gjrGARCH

delta - siła asymetrii (APARCH)

vxreg1, vxreg2, ... dla zmiennych egzogenicznych (zewnętrznych regresorów) modelu wariancji

ar1, ar2, ... dla części AR

ma1, ma2, ... dla części MA

mu dla wyrazu wolnego ARMAX

archm - wspomniany ARCH-in-Mean

arfima jako ułamkowy rząd d

mxreg1, mxreg2, ... dla zmiennych egzogenicznych (zewnętrznych regresorów) w ARMAX

eta11, eta12, ... dla parametrów asymetrii obrotu (FGARCH)
eta22, eta21, ... dla parametrów asymetrii przesunięcia (FGARCH)

lambda - siła warunkowej sigmy (wariancji)

skew - dla skośności rozkładu badanej zmiennej

shape - kształt rozkładu 

ghlambda - lambda dla rozkładu "ghyp"

Przykład 1, dla rozkładu normalnego: 

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

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

                                                                      alpha1 = 0.3, 

                                                                      beta1 = 0.2, 

                                                                      mu = 0.000001,

                                                                      ar1 = 0.000001,

                                                                      ma1 = 0.000001))


Przykład 2, dla skośnego rozkładu t-Studenta: 

ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 2)), 

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

                                                                      alpha1 = 0.3,

                                                                      beta1 = 0.2,

                                                                      beta2 = 0.1, 

                                                                      mu = 0.000001,

                                                                      ar1 = 0.000001,

                                                                      ma1 = 0.000001,

                                                                      skew = 1.5,

                                                                      shape = 5))

Symulacja procesu GARCH

Symulację wykonamy za pomocą funkcji ugarchpath(). Głównymi argumentami są wybrany model (ze specyfikacji) oraz liczba danych. Funkcja zwraca nie tyle samą symulację, ale cały zestaw jej obiektów. Aby uzyskać same wartości, najłatwiej użyć funkcji fitted(). 

Przykład:

set.seed(3)

mojGarch <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 0)), 

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

                                                                      alpha1 = 0.9, 

                                                                      mu = 0.000001,

                                                                      ar1 = 0.000001,

                                                                      ma1 = 0.000001))


symGarch <- ugarchpath(mojGarch, n.sim = 500)

symStopaZwrotu <- fitted(symGarch) 

# ewentualnie

#symStopaZwrotu <- symGarch@path$seriesSim

plot(symStopaZwrotu, type="l")

To co może zaskakiwać to to, że ustawiliśmy rozkład normalny, a wykres przypomina raczej rozkład z grubymi ogonami. Jednak jest to warunkowy rozkład normalny w tym sensie, że zależy od przeszłości. Sama zmienna (stopa zwrotu) jest tu niezależna od poprzednich wartości (ar1 i ma1  ustawione na zero), ale jej składnik losowy już tak. Składnik losowy jest po prostu procesem średniej kroczącej (normalnie nazywany MA, ale aby nie mylić z procesem dla samej zmiennej nazywany ARCH). 

Z symulacji możemy wydobyć też wariancję warunkową albo warunkowe odchylenie standardowe, tj. nasze ryzyko:

symGarch <- ugarchpath(mojGarch, n.sim = 500)
symStopaZwrotu <- fitted(symGarch) 
# ewentualnie
#symStopaZwrotu <- symGarch@path$seriesSim
plot(symStopaZwrotu, type="l")

symRyzyko <- sigma(symGarch)
# ewentualnie
#symRyzyko <- symGarch@path$sigmaSim
plot(symRyzyko, type="l")



Najlepiej połączyć obydwa wykresy w jeden: stopa zwrotu + jego zmienne odchylenie:

zakresY <- c(min(symStopaZwrotu, -symRyzyko), max(symStopaZwrotu, symRyzyko))

plot.ts(symStopaZwrotu, ylim = zakresY)

lines(symRyzyko, col = "blue")

lines(-symRyzyko, col = "blue")


Model empiryczny

Powiedzmy, że dostajemy samą taką symulację bez znajomości modelu. Chcemy dopasować otrzymane dwa szeregi do odpowiedniego modelu GARCH. Stosujemy funkcję ugarchfit(). Najważniejsze jej argumenty to dane (data), wybrany model (spec) oraz liczba prognoz poza próbą (out.sample). Zakładamy, że nic nie wiemy o próbie, więc na początek tworzymy standardowy sGARCH  - oczywiście tym razem bez ustalania parametrów, bo mamy zadanie je dopiero oszacować. Za dane podstawiamy wcześniej zdefiniowaną zmienną symStopaZwrotu:

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

dopasGarch <- ugarchfit(spec = mojGarch, data = simArch, out.sample = 0)

print(dopasGarch)



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

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,1)
Mean Model	: ARFIMA(1,0,1)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.039416    0.040145  0.98183  0.32618
ar1    -0.427181    0.391769 -1.09039  0.27554
ma1     0.480952    0.368062  1.30671  0.19131
omega   0.598849    0.071227  8.40760  0.00000
alpha1  0.804083    0.103489  7.76977  0.00000
beta1   0.005350    0.016619  0.32191  0.74752

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.039416    0.042395  0.92975  0.35250
ar1    -0.427181    0.358426 -1.19183  0.23333
ma1     0.480952    0.327567  1.46826  0.14203
omega   0.598849    0.065481  9.14534  0.00000
alpha1  0.804083    0.103111  7.79820  0.00000
beta1   0.005350    0.006337  0.84419  0.39856

LogLikelihood : -797.842 

Information Criteria
------------------------------------
                   
Akaike       3.2154
Bayes        3.2659
Shibata      3.2151
Hannan-Quinn 3.2352

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.067  0.3016
Lag[2*(p+q)+(p+q)-1][5]     1.357  0.9994
Lag[4*(p+q)+(p+q)-1][9]     3.017  0.8880
d.o.f=2
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                    0.02626  0.8713
Lag[2*(p+q)+(p+q)-1][5]   0.21540  0.9914
Lag[4*(p+q)+(p+q)-1][9]   1.84882  0.9223
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]  0.001223 0.500 2.000  0.9721
ARCH Lag[5]  0.274535 1.440 1.667  0.9470
ARCH Lag[7]  1.653454 2.315 1.543  0.7901

Nyblom stability test
------------------------------------
Joint Statistic:  1.0574
Individual Statistics:              
mu     0.03849
ar1    0.33674
ma1    0.31603
omega  0.38146
alpha1 0.21883
beta1  0.27590

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.49 1.68 2.12
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.9554 0.3399    
Negative Sign Bias  0.1454 0.8845    
Positive Sign Bias  0.2359 0.8136    
Joint Effect        1.2101 0.7506    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     16.00       0.6573
2    30     28.96       0.4672
3    40     34.40       0.6796
4    50     53.60       0.3023


Wydruk dostarcza nie tylko estymacji współczynników, ale też wszystkich potrzebnych informacji o poprawności modelu. Wykonujemy jego diagnostykę.

 Są dwie wersje testowania istotności parametrów. Pierwsza to konwencjonalna metoda największej wiarygodności (Optimal parameters), a druga odporna na nieprawidłowe założenie rozkładu (Robust Standard Errors). Ta druga jest więc mocniejsza i powinniśmy na nią patrzeć jako pierwszą. Parametry się w niej nie zmieniają, a jedynie ich przedziały istotności. Widzimy, że jedynie istotne są omega i alpha1, czyli poprawnie, natomiast same ich wartości też dość dobrze oszacowane (alpha1 jest równy 0.9, a oszacowano go na 0.8; omega równa się 0.5, a w estymacji = 0.6). Oczywiście widząc, że beta1 oraz ma1 są nieistotne, moglibyśmy poprawić teraz model. Dokończmy jednak najpierw diagnostykę.

Dalej mamy kryteria informacyjne (Information Criteria), które w tym wypadku nic nam nie powiedzą, bo nie porównujemy modelu z innym.

Ważony test Ljunga-Boxa sprawdzamy na dwóch poziomach: resztach (Weighted Ljung-Box Test on Standardized Residuals) oraz kwadratach reszt (Weighted Ljung-Box Test on Standardized Squared Residuals). Jeśli p-value > 0.05, to możemy przyjąć, że nie występują autokorelacje w pierwszym przypadku  w resztach, w drugim - w kwadratach reszt. Ponieważ p-value dla wielu rzędów q i p jest mocno powyżej 0.05, przyjęty model jest poprawny.

Podobnym testem jest Weighted ARCH LM, ale skupia się na samym efekcie ARCH, a dokładniej czy przyjęty model nie pominął tego efektu, np. dla pewnego rzędu. Jeżeli p-value > 0.05, to znaczy, że nie pominął, w przeciwnym wypadku prawdopodobnie pominął. Widzimy, że dla wielu rzędów q i p jest mocno powyżej 0.05, przyjęty model jest poprawny.

Test Nybloma (Nyblom stability test) bada stabilność modelu, tzn. czy parametry zmieniają się w czasie. Jest to bardzo ważny test z punktu widzenia właśnie analizy finansowej. Joint Statistic: 1.0574 dotyczy stabilności całego modelu. Porównujemy to z teoretycznymi wartościami dla kolejnych p-value, szczególnie patrzymy na wartość korespondującą z 5% - tutaj jest to 1.68. Ponieważ 1.68 > 1.0574 , nie ma podstaw do odrzucenia hipotezy zerowej o stabilności modelu. Dodatkowa analiza dotyczy każdego współczynnika osobno. Patrzymy na Individual Statistic. Mamy tu współczynniki: 

mu  =  0.03849, ar1 = 0.33674, ma1 = 0.31603, omega  0.38146, alpha1 0.21883, beta1  0.27590

Wartości te porównujemy z wartością teoretyczną odpowiadającą 5% - wynosi 0.47. Ponieważ 0.47 jest większa od każdego wspołczynnika, to znaczy, że wszystkie parametry są stabilne, zarówno średniej, wariancji, jak i rozkładu.

Test błędu znaku (Sign Bias Test) pozwala z kolei ocenić czy w szeregu nie występują szoki dodatnie lub ujemne. Niekonsekwentne oznaczenie  prob sig to nic innego jak p-value. Zarówno dla dodatniego, ujemnego, jak i łącznego efektu, jest ono większe od 0.05, czyli szoki nie występują. Gdyby występowały, trzeba by zmienić model na EGARCH, APARCH, TGARCH lub gjrGARCH.

W końcu test dobroci dopasowania modelu (Adjusted Pearson Goodness-of-Fit Test) sprawdza czy rozkład reszt jest odpowiedni. Mamy tu podział na kolejne kwantyle, tak że możemy zobaczyć jak reszty się zachowują w całym rozkładzie. Jeśli p-value > 0.05 dla któregoś kwantyla, znaczyłoby, że przyjęty rozkład jest niewłaściwy. W tym wypadku wszystkie wskazują poprawność rozkładu - przyjęliśmy normalny i oceniony jako też normalny.

Diagnostyka zakończona. Tutaj warto wspomnieć ciekawą opcję, którą można by nazwać graficzną diagnostyką. Wystarczy wpisać:

plot(dopasGarch)

i dostaniemy listę:

Wybierzmy 8:


Dzięki tej opcji możemy naocznie ocenić poprawność przyjętego rozkładu.

Również zamiast testu dobroci możemy użyć wykresu kwantylowego (QQ), tzn. wybieramy z listy 9:


Tak jak wspomniałem, skoro niektóre parametry okazały się nieistotne, to moglibyśmy teraz model poprawić, tak że kod wyglądałby następująco:

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

dopasGarch <- ugarchfit(spec = mojGarch, data = symStopaZwrotu, out.sample = 0)

print(dopasGarch)

Poniżej wydruk tylko kilku pierwszych statystyk:


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

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,0)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.039045    0.039600  0.98599  0.32414
omega   0.617144    0.067416  9.15423  0.00000
alpha1  0.800232    0.102546  7.80364  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.039045    0.042187  0.92552  0.35469
omega   0.617144    0.063213  9.76291  0.00000
alpha1  0.800232    0.102620  7.79800  0.00000

LogLikelihood : -799.3085 

Information Criteria
------------------------------------
                   
Akaike       3.2092
Bayes        3.2345
Shibata      3.2092
Hannan-Quinn 3.2192


Chociaż wiemy już, że drugi model jest lepszy, to aby ocenić to obiektywnie, porównamy ich kryteria informacyjne, które właśnie tutaj się przydadzą. Jeżeli wszystkie kryteria: AIC, BIC, SIC i HQIC wskazują mniejsze wartości, to znaczy, że ten model jest lepszy. Tak rzeczywiście się dzieje

Akaike       3.2092 < 3.2154
Bayes        3.2345 < 3.2659
Shibata      3.2092 < 3.2151
Hannan-Quinn 3.2192 < 3.2352

Otrzymujemy jednoznaczny wniosek, że ten drugi model jest lepszy od pierwszego.

Prognoza

Mając dobry model możemy przejść do prognozy. Prognozę można wykonać zarówno dla stopy zwrotu, jak i jej ryzyka. Ale ponieważ prognoza stopy zwrotu = oczekiwana stopa zwrotu, to znaczy, że w naszym przypadku jest to po prostu wyraz wolny bliski zera (nie mamy tu procesu ARMA), więc tak naprawdę interesuje nas tylko prognoza ryzyka. Prognozę wykonujemy przy pomocy ugarchforecast(). Jej główne argumenty to

1. model empiryczny z ugarchfit, ewentualnie model teoretyczny z podaniem danych (fitORspec)

2. ile wykonać prognoz do przodu (n.ahead)

3. ile prognoz ma się nakładać, aby stworzyć pierwszą poza próbą (n.roll)

4. ile ostatnich danych ma być zachowanych poza testowaną próbą (out.sample)

Punkty 3 i 4 są powiązane i na razie je zostawiam.

Stworzymy prognozę dla pięciu nowych punktów:

noweOkresy <- 5

prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = noweOkresy)

prognozaStopaZwrotu <- fitted(prognozaGarch)

prognozaRyzyko <- sigma(prognozaGarch)

Aby zobaczyć to wszystko graficznie, powinniśmy połączyć próbę z prognozą. Jednak łączenie różnych typów danych i pokazywanie na wykresie może być w R wyzwaniem. Ponieważ działamy na szeregach czasowych, to najpierw zamieńmy wszystkie dane na ts i połączymy z prognozą:

# zamiana na ts i łączenie z prognozą

stopaZwrotu <- as.ts(symStopaZwrotu)

ryzyko <- as.ts(symRyzyko)

ostatniOkres <- end(stopaZwrotu)

dlugoscProby <- length(stopaZwrotu)

stopaZwrotu <- ts(c(as.numeric(stopaZwrotu), rep(NA, noweOkresy)), start = start(stopaZwrotu))

prognozaStopaZwrotu <- ts(c(stopaZwrotu[dlugoscProby], as.numeric(prognozaStopaZwrotu)), start = ostatniOkres)

ryzyko <- ts(c(symRyzyko, rep(NA, noweOkresy)), start = start(stopaZwrotu))

prognozaRyzyko <- ts(c(ryzyko[dlugoscProby], as.numeric(prognozaRyzyko)), start = ostatniOkres)


Pokażmy 30 ostatnich obserwacji i 5 nowych

# wykres

zakresX <- c(dlugoscProby - 30, length(stopaZwrotu))

zakresY <- c(min(na.omit(stopaZwrotu), -na.omit(ryzyko)), max(na.omit(stopaZwrotu), -na.omit(ryzyko)) + 2)

plot(stopaZwrotu, xlim=zakresX, ylim=zakresY, lwd=2)

lines(prognozaStopaZwrotu, lwd=2, col="green")

lines(ryzyko, lwd=2, col="red")

lines(-ryzyko, lwd=2, col="red")

lines(prognozaRyzyko, lwd=2, col="brown")

lines(-prognozaRyzyko, lwd=2, col="brown")

legend(x="topleft", legend=c("stopa zwrotu", "prognoza stopy zwrotu", "modelowane ryzyko", 

                             "prognoza ryzyka"), lwd=c(2,2,2,2), col=c("black",

                                                                       "green",

                                                                       "red",

                                                                       "brown"))



Ostatnia symulacja testowa

Na koniec przeprowadzę ostatnią symulację i jej test. Tym razem na odwrót - stworzę zwykły proces AR(1) o 500 obserwacjach i sprawdzę jak rugarch radzi sobie z tym, gdy ma zły model.

set.seed(3)

symStopaZwrotu = arima.sim(list(order=c(1,0,0), ar=0.6), n=500, mean=1)

mojGarch <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(1, 1)), 

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

dopasGarch <- ugarchfit(spec = mojGarch, data = symStopaZwrotu, out.sample = 0)

print(dopasGarch) 

Częściowy wydruk:

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
mu      2.607355    0.115327    22.608  0.00000
ar1     0.605802    0.033249    18.220  0.00000
omega   0.001275    0.003960     0.322  0.74745
alpha1  0.000000    0.003644     0.000  1.00000
beta1   0.999000    0.000065 15463.568  0.00000

LogLikelihood : -724.7 

Information Criteria
------------------------------------
                   
Akaike       2.9188
Bayes        2.9610
Shibata      2.9186
Hannan-Quinn 2.9354

O ile algorytm prawidłowo oszacował ar1 na 0.6 i alpha1 na zero, o tyle błędnie przypisał becie 1. Ta pełna jedynka pokazuje, że nie możemy kierować się bezmyślnie samą istotnością: skoro beta równa się 1, to może znaczyć, że  model został źle wyspecyfikowany. Aby to ocenić, usuńmy betę ze specyfikacji:

mojGarch <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(1, 0)), 

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

dopasGarch <- ugarchfit(spec = mojGarch, data = symStopaZwrotu, out.sample = 0)

print(dopasGarch)

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu       2.60768    0.120767   21.593        0
ar1      0.60566    0.035427   17.096        0
omega    1.06510    0.081482   13.072        0
alpha1   0.00000    0.048245    0.000        1

LogLikelihood : -725 

Information Criteria
------------------------------------
                   
Akaike       2.9160
Bayes        2.9497
Shibata      2.9158
Hannan-Quinn 2.9292


Teraz widzimy, że model rzeczywiście był źle wyspecyfikowany, ponieważ wszystkie kryteria się poprawiły (mają mniejszą wartość). Jedynie omega zmieniła wartość, a więc to ją beta1 wcześniej psuła. Oczywiście w symulacji nie ustawialiśmy odchylenia standardowego, a więc zostało ono automatycznie ustalone na 1, a wobec tego jego kwadrat - omega - też musi się równać 1 i tyle też wyniosła estymacja (i również w końcu stała się istotna). Jedynym błędem jaki się nam ostał jest zbyt duży wyraz wolny (miał się równać 1). Ten błąd jednak nie wynika z błędnego algorytmu, bo gdy będziemy estymować zwykłą arimą:

arima(x=symStopaZwrotu, order=c(1,0,0))

dostaniemy to samo. Dopiero użycie MNK:

summary(lm(symStopaZwrotu[-1] ~ symStopaZwrotu[-length(symStopaZwrotu)]))

przyniesie prawidłowy wyraz wolny równy ok. 1. Na marginesie, w tej próbce akurat średnia arytmetyczna rzeczywiście równa się 2.6, co wynika po prostu z tego, że składa się ze średniej niewarunkowej i warunkowej. MNK lepiej oddziela te dwie części. Algorytmy dla ARIMA i GARCH są przeznaczone dla bardziej złożonych modeli, a nie zwykłego AR, stąd kosztem są takie błędy. 


Model na prawdziwych danych

Dzienne stopy zwrotu WIG20

Zakres czasu: 2023-04-21 - 2024-07-02, czyli ostatnie 300 sesji. Dane ściągnięte ze stooq.pl. Przygotowujemy analizę:

# ustawiamy folder roboczy

sciezka = r"(C:\testy)"

#zamieniam "\" na "/"

sciezka = gsub("\\", "/", sciezka, fixed=T)

setwd(sciezka_testy)

nazwaPliku = "wig20_d.csv"

plik <- read.csv(nazwaPliku)

cena <- plik$Zamkniecie

Obs <- 300

stopaZwrotu <- tail(diff(log(cena)), Obs)

okres <- as.Date(tail(plik$Data[-1], Obs))  

Zwracam uwagę, że używam logarytmicznej stopy zwrotu, co generalnie ułatwia nie tylko kodowanie, ale testowanie (zmienna staje się "bardziej" normalna). Ostatnie 300 stóp zwrotu WIG20 okazały się trudniejsze do modelowania niż sądziłem. Przykładowo, ten sam model co dla symulacji, okazał się niedobry, ponieważ parametr  ARCH(1) był nieistotny:

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

dopasGarch <- ugarchfit(spec = mojGarch, data = simArch, out.sample = 0)

print(dopasGarch)

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

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,0)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000933    0.000713   1.3079  0.19092
omega   0.000153    0.000014  11.0540  0.00000
alpha1  0.000000    0.039697   0.0000  1.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000933    0.000661   1.4103  0.15844
omega   0.000153    0.000016   9.3305  0.00000
alpha1  0.000000    0.033403   0.0000  1.00000

LogLikelihood : 892.5 

Information Criteria
------------------------------------
                    
Akaike       -5.9298
Bayes        -5.8928
Shibata      -5.9300
Hannan-Quinn -5.9150

Zobaczmy z jak ogromną kombinacją parametrów do określenia mamy tu do czynienia: rodzaj modelu zarówno GARCH, jak i ARFIMA (w tym określenie czy GARCH umieścić w ARMA czy nie) , rzędy GARCH(q, p), rząd AR, rząd MA, typ rozkładu.  Wykonałem wiele testów i doszedłem (przynajmniej na dziś) do wniosku, że najlepszym modelem jest TGARCH(1, 1)-in-Mean z rozkładem skośno-normalnym. TGARCH uwzględnia inne zachowanie się warunkowej wariancji, gdy są spadki i inne, gdy są wzrosty. Dostajemy więc dużo dokładniejsze zmiany ryzyka w kolejnych okresach niż w standardowym GARCH. Nazwa TGARCH pochodzi od Threshold GARCH, czyli progowy GARCH. Progiem jest tu... zero, poniżej którego następuje przeskoczenie na parametr dla ujemnych zmian (nie jest to nic skomplikowanego, po prostu rozdziela się składnik losowy na część dodatnią i ujemną).  

 mojGarch <- ugarchspec(variance.model = list(model = "fGARCH", submodel = "TGARCH",

                                             garchOrder = c(1, 1)), 

                       mean.model = list(armaOrder = c(0, 0), archm = T), distribution.model = "snorm")


dopasGarch <- ugarchfit(spec = mojGarch, data = stopaZwrotu, out.sample = 0)

print(dopasGarch)


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

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: fGARCH(1,1)
fGARCH Sub-Model	: TGARCH
Mean Model	: ARFIMA(0,0,0)
Distribution	: snorm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error   t value Pr(>|t|)
mu     -0.011509    0.000678   -16.983        0
archm   1.029024    0.057292    17.961        0
omega   0.000539    0.000038    14.082        0
alpha1  0.041037    0.003865    10.618        0
beta1   0.922511    0.000036 25297.879        0
eta11   1.000000    0.182817     5.470        0
skew    1.161197    0.095820    12.118        0

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
mu     -0.011509    0.000354   -32.488        0
archm   1.029024    0.023033    44.677        0
omega   0.000539    0.000018    30.293        0
alpha1  0.041037    0.001523    26.947        0
beta1   0.922511    0.000050 18517.286        0
eta11   1.000000    0.075490    13.247        0
skew    1.161197    0.089539    12.969        0

LogLikelihood : 902.5 

Information Criteria
------------------------------------
                    
Akaike       -5.9702
Bayes        -5.8838
Shibata      -5.9713
Hannan-Quinn -5.9356

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.1093  0.7409
Lag[2*(p+q)+(p+q)-1][2]    0.5386  0.6758
Lag[4*(p+q)+(p+q)-1][5]    1.9115  0.6395
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.715  0.1903
Lag[2*(p+q)+(p+q)-1][5]     2.352  0.5379
Lag[4*(p+q)+(p+q)-1][9]     3.638  0.6506
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.7304 0.500 2.000  0.3928
ARCH Lag[5]    1.2072 1.440 1.667  0.6726
ARCH Lag[7]    2.2430 2.315 1.543  0.6659

Nyblom stability test
------------------------------------
Joint Statistic:  1.287
Individual Statistics:              
mu     0.11514
archm  0.10348
omega  0.08008
alpha1 0.06191
beta1  0.06585
eta11  0.08665
skew   0.19526

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.69 1.9 2.35
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.3681 0.7131    
Negative Sign Bias  1.3763 0.1698    
Positive Sign Bias  1.1213 0.2631    
Joint Effect        3.1698 0.3662    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     14.67       0.7435
2    30     26.80       0.5825
3    40     31.47       0.7990
4    50     44.00       0.6756


Na pierwszy rzut oka może się wydać dziwne, że niewarunkowa średnia stopa zwrotu jest ujemna (mu = -0.0115), ale trzeba pamiętać, że model zawiera dodatkowy element - średnią warunkową zależną od warunkowego odchylenia standardowego i jej parametr jest dodatni (archm = 1,029). Powoduje to, że średnia stopa zwrotu zmienia się w czasie. Warto odnotować, że taki efekt może wcale nie zostać wykryty przez testy stabilności. 

Aby zobaczyć czy rozkład się zgadza, można użyć komendy

plot(dopasGarch, which=8)


i dodatkowo

plot(dopasGarch, which=9)


Wykonujemy prognozę. Zauważmy, że w tym przypadku możemy nieco poprawić precyzję prognozowanej stopy zwrotu: nie jest to po prostu wyraz wolny (średnia niewarunkowa), ale funkcja ryzyka (wariancji warunkowej). Z tego powodu wyłuskamy z empirycznego modelu nie tylko estymację ryzyka z próby, ale też "prognozy" (oczekiwanej stopy zwrotu) z próby. Prognozy poza próbą na 1 dzień do przodu także zrobimy.

# prognoza w próbie (dopasowanie do próby) 

stopaZwrotuDopas <- fitted(dopasGarch)

ryzyko <- sigma(dopasGarch)

# prognoza poza próbą

ileDoPrzodu <- 1

prognozaGarch <- ugarchforecast(dopasGarch, n.ahead = ileDoPrzodu)

prognozaStopaZwrotu <- fitted(prognozaGarch)

prognozaRyzyko <- sigma(prognozaGarch)

# zamiana na ts i łączenie z prognozą dla wykresu

stopaZwrotuWykres <- as.ts(stopaZwrotu)

stopaZwrotuDopasWykres <- as.ts(stopaZwrotuDopas)

ryzykoWykres <- as.ts(ryzyko)

ostatniOkres <- end(stopaZwrotuWykres)

dlugoscProby <- length(stopaZwrotuWykres)

stopaZwrotuWykres <- ts(c(as.numeric(stopaZwrotuWykres), rep(NA, ileDoPrzodu)), start = start(stopaZwrotuWykres))

stopaZwrotuDopasWykres <- ts(c(stopaZwrotuDopasWykres, rep(NA, ileDoPrzodu)), start = start(stopaZwrotuWykres))

prognozaStopaZwrotu <- ts(c(rep(NA, dlugoscProby-1), stopaZwrotuDopasWykres[dlugoscProby], 

                            as.numeric(prognozaStopaZwrotu)), start = stopaZwrotuWykres)

ryzykoWykres <- ts(c(ryzykoWykres, rep(NA, ileDoPrzodu)), start = start(stopaZwrotuWykres))

prognozaRyzyko <- ts(c(rep(NA, dlugoscProby-1), ryzykoWykres[dlugoscProby], as.numeric(prognozaRyzyko)), start = stopaZwrotuWykres)

okres <- tail(plik$Data[-1], Obs)

okresPrognoza <- max(okres) + days(1)

okresWykres <- c(okres, okresPrognoza)

# wykres

Ost <- 15

zakresX <- c(okresWykres[dlugoscProby - Ost + 1], okresWykres[length(okresWykres)])

zakresY <- c(min(na.omit(stopaZwrotuWykres), -na.omit(ryzykoWykres)), max(na.omit(stopaZwrotuWykres), -na.omit(ryzykoWykres)) + max(na.omit(stopaZwrotuWykres), -na.omit(ryzykoWykres)))

plot(x=okresWykres, y=stopaZwrotuWykres, main="Dzienna stopa zwrotu WIG20", xaxt="n", xlim=zakresX, ylim=zakresY, lwd=2, type="l", xlab="")

axis(side=1, at=okresWykres, labels=okresWykres, las=2, cex.axis=0.8)

lines(x=okresWykres, y=stopaZwrotuDopasWykres, col="lightblue", lwd=2)

lines(x=okresWykres, y=prognozaStopaZwrotu, lwd=2, col="green")

lines(x=okresWykres, y=ryzykoWykres, lwd=2, col="red")

lines(x=okresWykres, y=-ryzykoWykres, lwd=2, col="red")

lines(x=okresWykres, y=prognozaRyzyko, lwd=2, col="brown")

lines(x=okresWykres, y=-prognozaRyzyko, lwd=2, col="brown")

abline(h=0, col="grey")

legend(x="topleft", legend=c("stopa zwrotu", "modelowana stopa zwrotu (GARCH-In-Mean)", "prognoza stopy zwrotu", "modelowane ryzyko (T-GARCH)", 

                             "prognoza ryzyka (T-GARCH)"), lwd=c(2,2,2,2,2), col=c("black",

                                                                       "lightblue",

                                                                       "green",

                                                                       "red",

                                                                       "brown"))



Miesięczne stopy zwrotu

Aby zagregować dane dzienne do miesięcznych użyłem pakietu xts. Obserwacji jest tak samo 300, więc cofamy się znacznie dalej w przeszłość. Zakres to lipiec 1999 - czerwiec 2024.

tabela <- as.data.frame(to.monthly(as.xts(plik), indexAt="endof", OHLC=FALSE))

cena <- tabela$Zamkniecie

Obs <- 300

stopaZwrotu <- tail(diff(log(cena)), Obs)

okres <- as.Date(tail(rownames(tabela)[-1], Obs))


Tym razem lepsze rezultaty daje EGARCH:

mojGarch <- ugarchspec(variance.model = list(model = "eGARCH", submodel = NULL,

                                             garchOrder = c(1, 1)), 

                       mean.model = list(armaOrder = c(0, 0), archm = T), distribution.model = "snorm")


I dalej już analogicznie do dziennych. Jedyną różnicą w kodzie będzie zastąpienie funkcji days(1) przez months(1).

okresPrognoza <- max(okres) + months(1)

Częściowy wydruk:

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu     -0.046871    0.002481 -18.8899 0.000000
archm   0.740374    0.038708  19.1271 0.000000
omega  -0.340842    0.042694  -7.9834 0.000000
alpha1 -0.081079    0.021606  -3.7526 0.000175
beta1   0.938421    0.007192 130.4857 0.000000
gamma1  0.129712    0.038211   3.3946 0.000687
skew    0.840676    0.063438  13.2519 0.000000


I wykresy





W pełnej skali analizowanego okresu dostaniemy następującą grafikę:


Obecnie mamy do czynienia ze spadającą zmiennością. Nie jest to przypadek, a zjawisko związane z grupowaniem wariancji. Po okresie dużej zmienności przychodzi mała (w podsumowaniu wyjaśnienie możliwej przyczyny).  Z tego powodu model prognozuje spadek stopy zwrotu w lipcu, a nawet indeksu. Nie podaję wartości, bo to nie ma większego sensu - widzimy, że nie jest to prawdziwa prognoza, tylko wartość oczekiwana, średnia. 

Podsumowanie

W ramach podsumowania oto wnioski i przemyślenia:

1. dla dziennych i miesięcznych stóp zwrotu WIG20 występuje efekt GARCH(1, 1). Dla ostatnich 300 dziennych lepszy jest TGARCH, dla ostatnich 300 miesięcznych EGARCH. Nie jest to przypadek. TGARCH lepiej radzi sobie z asymetrią wpływu pozytywnych i negatywnych informacji. Jak łatwo się domyślić, powinno (i jest) to bardziej widoczne dla zmian krótkoterminowych, bo wiadomości zmieniają się z dnia na dzień. Z kolei EGARCH będzie lepszy tam gdzie ujawnia się zjawisko grupowania wariancji. To bardziej tajemnicza sprawa, jak sądzę jest związana z cyklicznością przepływu kapitału. Kasa przecież generalnie spływa co miesiąc do gospodarki, a kalendarz wszyscy mają mniej więcej taki sam. Więcej kapitału to większa liczba inwestorów, to większa liczba wycen, a co za tym idzie większa rozbieżność opinii, stąd większa zmienność. Stąd oczywiste, że zjawisko grupowania wariancji prędzej będzie się odbywać w interwałach miesięcznych niż dziennych.  

2. GARCH stwierdza, że zarówno dzienne, jak i miesięczne stopy zwrotu WIG20 nie zawierają autokorelacji, natomiast korelują dodatnio z ryzykiem, tj. mamy do czynienia z GARCH-in-Mean. To dobra wiadomość dla teoretyków finansów, bo po pierwsze w końcu mamy dowód, że CAPM na polskim rynku działa oraz zachowana zostaje jego efektywność (bo właśnie nie ma autokorelacji, a jest dodatni związek między ryzykiem a oczekiwaną stopą zwrotu).

3. Chociaż zdecydowanie należy odrzucić rozkład normalny, to dla log-stóp zwrotu można spokojnie przyjąć skośno-normalny.

4. Biorąc pod uwagę poprzednie punkty, należy stwierdzić że mamy do czynienia z warunkowym CAPM ze skośnością oraz efektami GARCH.

P. S. Jeżeli test Nybloma wskazuje , że model jest niestabilny i poprawianie założeń tego nie zmienia, należy skorzystać z tych testów. Alternatywą jest jeszcze CUSUMQ lub ogólniejsza rodzina testów fluktuacji (dla kwadratów stopy zwrotu), które w pewnym sensie łączą warunkowe i niewarunkowe miary.

poniedziałek, 20 maja 2024

Ile powinien kosztować Orlen?

Orlen wygląda jakby wybudził się z marazmu. Przyjrzyjmy mu się.  

Z technicznego punktu widzenia jest tak sobie. Kurs znajduje się w obszarach różnych oporów, a RSI sugeruje wykupienie. Ostatnie kilka szczytów RSI w okolicy 70 dawało faktycznie sygnał sprzedaży i teraz jest w tej okolicy: 


Co do wyceny różni analitycy próbują nam wmówić, że Orlen jest za tani. Łatwo obalić to twierdzenie. W 2023 r. spółka ogłosiła nową strategię dywidendową polegającą na wypłacie 40% skorygowanych wolnych przepływów pieniężnych (tj. przepływy z działalności operacyjnej minus przepływy z działalności inwestycyjnej) z poprzedniego roku. Jednak nie będzie to mniej niż dywidenda bazowa (gwarantowana), której poziom został ustalony na 4,00 PLN na jedną akcję dla roku 2022 i będzie sukcesywnie rósł o 15 groszy co roku, aż do poziomu 5,20 PLN na jedną akcję w roku 2030 (źródło). Wynika z tego, że do wzoru na model zdyskontowanych dywidend można wstawić minimalne dywidendy, ale wtedy zamiast koszt kapitału własnego musimy ustawić rentowność zbliżoną do obligacji. Nie można jednak brać obligacji samego Orlenu, którego bieżąca rentowność wynosi poniżej 7%, bo jego zapadalność to grudzień 2025. Lepszym porównaniem jest KGHM, którego obligacje będą wykupywane w czerwcu 2029, dla których rentowność wynosi obecnie 7,5% (tyle co oprocentowanie). 

Te 7,5% to ciągle za mało. Po pierwsze np. Millenium też do 2029 płaci ponad 8%. A przecież, to po drugie, mówimy o rynku dłużnym. Dług zawsze będzie miał pierwszeństwo przed dywidendą. Orlen mówi o "dywidendzie gwarantowanej", ale tak naprawdę co to znaczy? Nie mamy przecież gwarancji w sensie prawnym. Jest to tylko trochę bardziej pewne niż "na gębę". 

Jest też taka spółka jak PCC Rokita, która ma wypracowaną markę w branży chemicznej i która płaci ok. 9%. 

Oczywiście tak wysoki procent wynika z wysokich stóp proc., które mogą spaść w ciągu najbliższych 5 lat. Ocena tej przyszłości ciągle się zmienia przez rynek (stąd m.in. wahania kursu), więc pozostaje jedynie sprawdzić różne scenariusze.

Dla dochodu rezydualnego (po 2031) utrzymam te same parametry co w latach poprzednich. 

W sumie założyłem, że różnica między minimalną stopą dyskontową a stopą wolną od ryzyka wynosi 2,2%. 

I tak, przy dzisiejszej stopie wolnej od ryzyka dla obligacji skarbowych 7% ustaliłem stopę dyskontową (dla minimalnej dywidendy) 9,2%. Otrzymałem wartość ponad 70 zł dla inwestora krótkoterminowego lub po prostu płacącego koszty transakcyjne. Jednak koszty te nie muszą być rozumiane dosłownie, bo może to być koszt analizy akcji (czy sprzedać, dokupić itp.). 


Jak widać uwzględniłem dywidendę na rok 2024, która nie jest dyskontowana.

Dla inwestora średnioterminowego i długoterminowego  ponoszącego koszty transakcyjne (prowizja maklerska lub koszty analizy momentu sprzedaży akcji) dostałem niecałe 73 zł. 




Dla inwestora dywidendowego (bez kosztów transakcji) - ok. 75 zł.



W scenariuszu zakładającym spadek stopy procentowej w kolejnych latach o 0,5 pkt. proc. wartość akcji rośnie nawet do 81 zł, a przy założeniu stopniowego spadku stopy o 1,5 pkt - do 98. Takie wyceny będą możliwe przy większym wolumenie, ponieważ pojawia się wtedy więcej graczy, a tym samym zwiększa się rozrzut wycen.

Ryzykanci mogą więc próbować jeszcze kupować. Natomiast inwestorzy konserwatywni powinni raczej szukać gdzie indziej. Pozornie niska wycena wynika po prostu ze słabych (w porównaniu z zyskami) dywidend. Dzisiejszy kurs 72,5 zł stanowi rozsądną wycenę w porównaniu do 65 zł jeszcze kilka tygodni temu.

Link do analizy