Analiza porównawcza ryzyka zrównoważonego portfela w wybranych modelach klasy GARCH

1 Streszczenie

Niniejsza praca ma na celu analizę porównawczą ryzyka rozumianego jako oszacowanie funkcji warunkowej wariancji w modelach klasy GARCH. Do przeprowadzenia badania wybrałem dwa rozszerzenia modeli GARCH: TGARCH i EGARCH, które zaimplementowane są w pakiecie rugarch. Zbudowałem zrównoważony portfel (ang. equally-weighted portfolio) składający się z 5 kryptowalut o największej kapitalizacji. Dla tak zbudowanego portfela przeprowadziłem analizę porównawczą oszacowań warunkowej wariancji oraz oszacowań wartości narażonej na ryzyko (ang. Value-at-Risk, VaR) zarówno w wybranym okresie in-sample jak i w okresie out-of-sample, bazując na jednookresowych prognozach funkcji warunkowej wariancji. Przeprowadziłem również analizę wrażliwości, która pokazuje jak przyjęcie początkowych założeń analizy wpływa na końcowy wynik i wnioski.

2 Wstęp

Celem poniższego badania była analiza porównawcza ryzka w portfelu zrównoważonym. Do utworzenia portfela wybrałem 5 kryptowalut o największej kapitalizacji. Postanowiłem skupić swoją uwagę na kryptowalutach, ze względu na olbrzymią zmienność jaką charakteryzuje się ten rynek oraz jego rosnącą popularność. Aby odpowiednio określić stopień ryzyka z utworzonego portfela inwestycyjnego, posłużyłem się oszacowaniem funkcji warunkowej wariancji oraz wartości narażonej na ryzyko dla modeli klasy GARCH. Następnie porównałem ze sobą otrzymane dla wybranych modeli wyniki.

2.1 Kryptowaluty

Kryptowaluty to waluty wirtualne, umożliwiające transfer pieniężny pomiędzy użytkownikami sieci bez angażowania pośredników, takich jak banki czy kantory. Ich główną zaletą jest zdecentralizowany system transakcyjny oparty o technologię blockchain. Do przechowywania transakcji używana jest rozproszona baza danych, rozprowadzona pomiędzy węzłami sieci peer-to-peer. Dzięki temu każdy użytkownik ma możliwość pobrania jej, a wszystkie zawierane transakcje są jawne, lecz identyfikatory użytkowników mogą pozostać anonimowe. Taki system sprzyja bezpieczeństwu transakcji oraz przejrzystości przepływów finansowych.
Najpopularniejszą, a zarazem pierwszą w obiegu, kryptowalutą jest Bitcoin. Kolejnymi kryptowalutami o największej rynkowej kapitalizacji są: Ethereum (Ether), Ripple, Litecoin oraz NEM. W tej chwili (stan na 22.07.2017) liczba kryptowalut przekroczyła już 1000 i w dalszym ciągu rośnie. Na rynku tym w niedalekiej przyszłości może dojść do znaczących zmian popularności poszczególnych kryptowalut, w związku z planowym podziałem Bitcoin’a na dwie oddzielne kryptowaluty: Bitcoin i Bitcoin Cash.

2.2 Modele klasy GARCH

Modele klasy GARCH to modele uogólnionej warunkowej autoregresyjnej heteroskedastyczności. Są one stosowane do przewidywania zmian wariancji błędu losowego.
Z uwagi na fakt, że szeregi finansowe cechują się występowaniem różnych efektów:

• efekt leptokurtozy - występowanie “grubych ogonów” i wyższego szczytu funkcji gęstości, niż to ma miejsce w rozkładzie normalnym;
• efekt grupowania wariancji - występowanie po sobie okresów nasilonej zmienności i okresów względnie stabilnych;
• efekt dźwigni - asymetryczny wpływ informacji pozytywnych i negatywnych na poziom przyszłej wariancji;
• efekt skośności - brak symetrii rozkładu stóp zwrotu względem średniej;
• efekt autokorelacji stóp zwrotu;
• efekt długiej pamięci - po znaczących wzrostach występują kolejne wzrosty, następnie nadchodzą znaczące spadki, a po nich dalsze spadki;

przyjęło się wykorzystywać różne rozszerzenia standardowego modelu GARCH, w celu jak najdokładniejszego opisu badanego zjawiska. W poniższym opracowaniu skupiłem się na analizie rozszerzeń uwzględniających efekt dźwigni (TGARCH, EGARCH) i porównaniu otrzymanych wyników ze standardowym modelem GARCH.

3 Budowa portfela

3.1 Wstępna analiza i uporządkowanie danych

Pierwszym punktem wykonywanego badania było zbudowanie zrównoważonego portfela inwestycyjnego. Aby przeprowadzić analizę skorzystałem z następujących pakietów:

library(dplyr)
library(ggplot2)
library(stargazer)
library(tidyr)
library(knitr)
library(pander)
library(DT)
library(plyr)
library(stats)
library(PerformanceAnalytics)
library(xts)
library(rugarch)
library(fGarch)
library(timeSeries)
library(quantmod)
library(tseries)
library(FinTS)
library(TSA)
library(reshape)
library(dygraphs)
library(parallel)
library(xtable)
library(grid)

Następnie wczytałem dane dotyczące notowań kryptowalut, które pochodzą ze strony coinmarketcap.com.

bitcoin.url <- "https://docs.google.com/spreadsheets/d/1nxJm_h7jTH4p97DxZySdPGrlD96fhVwRGBsnd27d8W4/pub?output=csv"
ethereum.url <- "https://docs.google.com/spreadsheets/d/1qYN1gkMUh5o3CkaOjFEJp58CLbF6mQ7qPEhxOYpDYtQ/pub?output=csv"
ripple.url <- "https://docs.google.com/spreadsheets/d/1Da0PL2nV48OTNWM-MU9zIWTo8jVz_To9QOkDuInHyQU/pub?output=csv"
litecoin.url <- "https://docs.google.com/spreadsheets/d/1brWFRaCK9z7CPN-fzbibYz00GYRbMWyZqdh0gHA2q1o/pub?output=csv"
nem.url <- "https://docs.google.com/spreadsheets/d/1qJwh3MTIe9wPHuoabyOVTTuVXF302ZHe7g0UhxwDNAk/pub?output=csv"

bitcoin <- read.table(bitcoin.url, head = TRUE, sep = ";", stringsAsFactors = FALSE)
ethereum <- read.table(ethereum.url, head = TRUE, sep = ";", stringsAsFactors = FALSE)
ripple <- read.table(ripple.url, head = TRUE, sep = ";", stringsAsFactors = FALSE)
litecoin <- read.table(litecoin.url, head = TRUE, sep = ";", stringsAsFactors = FALSE)
nem <- read.table(nem.url, head = TRUE, sep = ";", stringsAsFactors = FALSE)

Kolumny zawierające daty należało przekształić do formatu daty w R.

lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C")
bitcoin$Date <- as.Date(bitcoin$Date, format="%b %d, %Y")
ethereum$Date <- as.Date(ethereum$Date, format="%b %d, %Y")
ripple$Date <- as.Date(ripple$Date, format="%b %d, %Y")
litecoin$Date <- as.Date(litecoin$Date, format="%b %d, %Y")
nem$Date <- as.Date(nem$Date, format="%b %d, %Y")
Sys.setlocale("LC_TIME", lct)

Z uwagi na fakt, że dane dla poszczególnych kryptowalut różniły się między sobą przedziałami czasowymi:

min(bitcoin$Date) 

## [1] "2015-05-22"

min(ethereum$Date)

## [1] "2015-08-08"

min(ripple$Date)

## [1] "2015-05-22"

min(litecoin$Date)

## [1] "2015-05-22"

min(nem$Date)

## [1] "2015-06-22"

ustaliłem arbitralnie początek próby na dzień 10.08.2015 r.

bitcoin <- bitcoin[bitcoin$Date >= as.Date("2015-08-10"),]
ethereum <- ethereum[ethereum$Date >= as.Date("2015-08-10"),]
ripple <- ripple[ripple$Date >= as.Date("2015-08-10"),]
litecoin <- litecoin[litecoin$Date >= as.Date("2015-08-10"),]
nem <- nem[nem$Date >= as.Date("2015-08-10"),]

Za pomocą funkcji xts z pakietu xts przekształciłem powyższe ramki danych w szeregi czasowe.

bitcoinTS <- xts(bitcoin$Close, bitcoin$Date)
ethereumTS <- xts(ethereum$Close, ethereum$Date)
rippleTS <- xts(ripple$Close, ripple$Date)
litecoinTS <- xts(litecoin$Close, litecoin$Date)
nemTS <- xts(nem$Close, nem$Date)

Następnie korzystając z pakietu dygraphs stworzyłem dynamiczne wykresy pokazujące jak zmieniała się wycena poszczególnych kryptowalut w całym okresie trwania próby. Przykładowy kod do stworzenia wykresu:

dygraph(bitcoinTS, main="Kurs bitcoina w dolarach") %>% 
  dyRangeSelector() %>%
  dyOptions(stackedGraph = TRUE) %>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 1.5)) %>%
  dyAxis("y", label = "Price ($)") %>%
  dySeries(color = "#CCCC00") 

3.2 Zrównoważony portfel inwestycyjny

Wartość danego aktywa w portfelu definiowana jest jako: \[V_{i} = \lambda_{i}P_{i},\]
gdzie \(\lambda_{i}\) - liczba instrumentów danego typu w portfelu, \(P_{i}\) - wartość poszczególnego instrumentu. Całkowita wartość portfela:
\[V_{P} = \sum_{i=1}^{N}V_{i}.\]
Waga i-tego aktywa w portfelu definiowana jest jako: \[w_{i} = \frac{V_{i}}{V_{P}},\]
natomiast stopa zwrotu z portfela w czasie \(t\):
\[R_{t} = \frac{V_{P_{t}} - V_{P_{t-1}}}{V_{P_{t-1}}}.\] Posiłkując się powyższymi definicjami, stworzyłem zrównoważony portfel przy założeniu, że jego początkowa wartość wynosiła $1000. Ponieważ instrumenty mają posiadać jednakową wagę to \(w=0.2\). Jako \(k\) oznaczyłem ilość instrumentów zakupionych w dniu pierwszym.

ceny <- cbind(bitcoinTS, ethereumTS, rippleTS, litecoinTS, nemTS)
V_0 <- 1000
N <- ncol(ceny) 
w <- rep(1/N, N)
k <- w * V_0 / ceny[1,]
k_num <- as.numeric(k)

3.3 Wartość portfela w kolejnych okresach

V_t <- matrix(0, nrow(ceny), ncol(ceny), dimnames=dimnames(ceny))
for(i in 1:nrow(ceny)){
  V_t[i,] = ceny[i,] * k_num
} 
V_tTS <- data.frame(V_t)
V_tTS <- xts(V_tTS, rev(bitcoin$Date))

V_portfel <- rowSums(V_t) 
V_portfelTS <- data.frame(V_portfel)
V_portfelTS <- xts(V_portfelTS, rev(bitcoin$Date)) 

w_t <- V_t / V_portfel
w_tTS <- data.frame(w_t)
w_tTS <- xts(w_tTS, rev(bitcoin$Date))

3.4 Zwroty z portfela

R_t <- diff(V_portfel) / V_portfel[1:713]
R_tTS <- data.frame(R_t)
R_tTS <- xts(R_tTS, rev(bitcoin$Date[1:713]))

Zwroty logarytmiczne różnią się tym od zwrotów “prostych”, że są zawsze mniejsze. Najmniejszy możliwy zwrot prosty to -100% podczas gdy najmniejszy możliwy zwrot logarytmiczny to -infinity. Dla dziennych zwrotów oraz zwrotów o krótszych okresach niż dzienne (np. godzinowe) różnica pomiędzy oboma typami zwrotów nie będzie znacząca, ale w analizach przyjęło się powszechnie używać zwrotów logarytmicznych. Aby wyznaczyć logarytmiczne zwroty użyłem pakietu quantmod.

dR <- dailyReturn(V_portfelTS, subset=NULL, type="log")

W celu zobrazowania poszczególnych właściwości stóp zwrotów w analizowanym okresie, wywołałem tabelę z najpopularniejszymi statystykami.

table.Stats(dR)

##                 daily.returns
## Observations         713.0000
## NAs                    0.0000
## Minimum               -0.2184
## Quartile 1            -0.0224
## Median                 0.0020
## Arithmetic Mean        0.0082
## Geometric Mean         0.0063
## Quartile 3             0.0329
## Maximum                0.5054
## SE Mean                0.0023
## LCL Mean (0.95)        0.0036
## UCL Mean (0.95)        0.0127
## Variance               0.0038
## Stdev                  0.0618
## Skewness               1.5473
## Kurtosis               9.9424

Jak można zaobserwować maksymalna dzienna strata osiąga wielkość 0.2184, natomiast maksymalny dzienny zysk 0.5054, a średni dzienny zysk wynosi 0.0082 przy wariancji równej 0.0038. Ważną kwestią jest wysoka wartość kurtozy, która sugeruje iż zwroty nie przyjmują rozkładu normalnego. Aby to sprawdzić wykonałem test Jarque-Bera na normalność rozkładu.

jarque.bera.test(dR)

## 
##  Jarque Bera Test
## 
## data:  dR
## X-squared = 3221.2, df = 2, p-value < 2.2e-16

Ponieważ p-value osiągnęło wartość dużo niższą od przyjętego poziomu istotności rzędu 0.05, należy odrzucić hipotezę zerową o normalności rozkładu logarytmicznych stóp zwrotu. W celu lepszego zobrazowania rozkładu stóp zwrotu, wykonałem poniższe wykresy.

z <- as.numeric(dR)
y <- quantile(z, c(0.25, 0.75))
x <- qnorm(c(0.25, 0.75),mean(z),sd(z))
slope <- diff(y)/diff(x)
int <- y[1L] - slope * x[1L]
D <- with(density(z), data.frame(x,y)) 

QQ <- ggplot(NULL, aes(sample=z)) + 
  stat_qq(distribution=qnorm, dparams=list(mean=mean(z),sd=sd(z)), colour="red3") +
  geom_abline(slope=slope, intercept=int, size=1)

D_plot <- ggplot(data=D, aes(x=x, y=y)) +
  geom_line(lwd=0.2,col="red3") +
  geom_area(aes(fill="red"), alpha=0.5) +
  stat_function(fun=dnorm, args=list(mean=mean(z),sd=sd(z)), colour= "black") +
  theme(legend.position="none")

pushViewport(viewport(layout=grid.layout(1, 2)))
print(QQ, vp=viewport(layout.pos.row=1, layout.pos.col=1))
print(D_plot, vp=viewport(layout.pos.row=1, layout.pos.col=2))

Widać wyraźnie, że rozkład zwrotów jest silnie leptokurtyczny. Cechuje się on “grubymi ogonami” oraz dużym skupieniem wokół średniej. Należy więc przypuszczać, że do jego prawidłowego zobrazowania powinno się posługiwać raczej rozkładem t-Studenta aniżeli rozkładem normalnym.

4 Analiza portfela

4.1 Testy efektów heteroskedastycznych

Poniżej zaprezentowany jest wykres funkcji ACF kwadratów zwrotów z portfela (ang. Autocorrelation Function), z którego można odczytać czy stopa zwrotów podlega autokorelacji.

lag <- 1:length(acf(dR^2, plot=FALSE)$lag[-1])
acf <- acf(dR^2, plot=FALSE)$acf[-1]
C <- 1.96/sqrt(length(dR^2))

ggplot(NULL, aes(x=lag, y=acf)) +
  geom_bar(stat="identity", position= "identity") +
  geom_hline(yintercept=C, color="red3", size=0.5)+
  geom_hline(yintercept=-C, color="red3", size=0.5)

Z wykresu można wywnioskować, że logarytmiczna stopa zwrotu podlega stosunkowo silnej autokorealcji. Aby sprawdzić to przypuszczenie, wykonałem testy Ljungi-Boxa na występowanie autokorelacji kwadratów reszt dla pierwszych pięciu opóźnień.

Ljung_Box.statistic <- c(Box.test(dR^2, lag=1, type="Ljung-Box")$statistic, Box.test(dR^2, lag=2, type="Ljung-Box")$statistic, Box.test(dR^2, lag=3, type="Ljung-Box")$statistic, Box.test(dR^2, lag=4, type="Ljung-Box")$statistic, Box.test(dR^2, lag=5, type="Ljung-Box")$statistic) 

Ljung_Box.pvalue <- c(Box.test(dR^2, lag=1, type="Ljung-Box")$p.value, Box.test(dR^2, lag=2, type="Ljung-Box")$p.value, Box.test(dR^2, lag=3, type="Ljung-Box")$p.value, Box.test(dR^2, lag=4, type="Ljung-Box")$p.value, Box.test(dR^2, lag=5, type="Ljung-Box")$p.value)

Ljung_Box <- cbind(Ljung_Box.statistic, Ljung_Box.pvalue)
rownames(Ljung_Box) <- c("lag=1", "lag=2", "lag=3", "lag=4", "lag=5")
colnames(Ljung_Box) <- c("X-squared", "p-value")

kable(Ljung_Box)

	X-squared	p-value
lag=1	3.266610	0.0707036
lag=2	7.087257	0.0289082
lag=3	21.168989	0.0000971
lag=4	23.839545	0.0000860
lag=5	29.794723	0.0000162

Dla każdego z badanych opóźnień, pomijając pierwsze, występuje zjawisko autokorelacji. Następnie został przeprowadzony test mnożników Lagrange’a na występowanie efektów ARCH.

Arch_Test.statistic <- c(ArchTest(dR, lags=1)$statistic, ArchTest(dR, lags=2)$statistic, ArchTest(dR, lags=3)$statistic, ArchTest(dR, lags=4)$statistic, ArchTest(dR, lags=5)$statistic)

Arch_Test.pvalue <- c(ArchTest(dR, lags=1)$p.value, ArchTest(dR, lags=2)$p.value, ArchTest(dR, lags=3)$p.value, ArchTest(dR, lags=4)$p.value, ArchTest(dR, lags=5)$p.value)

Arch_Test <- cbind(Arch_Test.statistic, Arch_Test.pvalue)
rownames(Arch_Test) <- c("lag=1", "lag=2", "lag=3", "lag=4", "lag=5")
colnames(Arch_Test) <- c("Chi-squared", "p-value")

kable(Arch_Test)

	Chi-squared	p-value
lag=1	3.264077	0.0708129
lag=2	6.640547	0.0361430
lag=3	19.057598	0.0002660
lag=4	20.257913	0.0004441
lag=5	23.693631	0.0002486

Z tabeli znajdującej się powyżej należy wyciągnąć wniosek, że dla pierwszego opóźnienia nie można odrzucić \(H_{0}\) o niewystępowaniu efektów ARCH. Natomiast dla każdego z kolejnych opóźnień p-value osiąga wartości mniejsze od 0.05, a więc należy odrzucić \(H_{0}\) o niewystępowaniu efektów ARCH. W takim wypadku prawidłowym działaniem w celu modelowania zmienności logarytmicznych stóp zwrotu z portfela będzie użycie modeli klasy GARCH.

4.2 Podział próby na in-sample i out-of-sample

Przystępując do dalszej części analizy, dokonałem podziału próby na okres in-sample i out-of-sample. Mając na uwadze zachowanie odpowiednich proporcji przy dostatecznie długim okresie in-sample względem wyjściowej próby, ustaliłem, że okres ten będzie przypadał na przedział 10.08.2015 - 21.05.2017 (1 rok 9 miesięcy 11 dni). Natomiast okres out-of-sample będzie się mieścił w granicach 22.05.2017 - 22.07.2017 (2 miesiące). Podział ten daje stosunek okresu in-sample do okresu out-of-sample na poziomie bliskim 2:21 \(\approx\) 1:11.

portfel <- cbind(rev(bitcoin$Date), data.frame(V_portfel))
colnames(portfel) <- c("Date", "Value")
insample.df <- portfel[portfel$Date < as.Date("2017-05-22"),] #10.08.2015 - 21.05.2017 (1 rok 9 miesięcy 11 dni)
outsample.df <- portfel[portfel$Date >= as.Date("2017-05-22"),] #22.05.2017 - 22.07-2017 (2 miesiące)

insample <- xts(insample.df$Value, insample.df$Date)
colnames(insample) <- "Value"
dR.in <- dailyReturn(insample, subset=NULL, type="log")

4.3 Wybór parametrów modelu

W celu właściwego wykorzystania modli klasy ARIMA i GARCH, kluczowym elementem jest określenie rzędu p części autoregresyjnej modelu AR oraz rzędu q części związanej ze średnią ruchomą modelu MA.

4.3.1 ARIMA(p,q)

Poniżej przedstawię działanie oraz wynik algorytmu pozwalającego na wybór parametrów p i q w modelu ARIMA. Do oceny, który model lepiej dopasowuje się do użytych danych wykorzystam kryterium informacyjne Akaike.

cl <- makePSOCKcluster(10)
AC <- autoarfima(as.numeric(dR.in), ar.max=2, ma.max=2, criterion='AIC', method='partial', arfima=FALSE, include.mean=NULL, distribution.model='std', solver='solnp', cluster=cl)

kable(head(AC$rank.matrix))

AR	MA	Mean	ARFIMA	AIC	converged
2	2	0	0	-3.129875	1
0	0	1	0	-3.126914	1
1	1	1	0	-3.125816	1
1	1	0	0	-3.124229	1
0	1	1	0	-3.124067	1
1	0	1	0	-3.124058	1

Z uwagi na fakt, iż rozkład logarytmicznych stóp zwrotu z portfela jest silnie leptokurtyczny, rozkład t-Studenta wydaje się lepiej dopasowanym do wyników rozkładem niż rozkład normalny. Tak więc powyższy wybór parametrów p i q przeprowadziłem zakładając, że logarytmiczne stopy zwrotu mają rozkład t-Studenta. Na podstawie kryterium AIC jako najlepszy model wyłoniony został ARIMA(2,2).

4.3.2 GARCH(p,q)

Następnie przeprowadziłem estymację modelu GARCH z różnymi wariantami parametrów p i q oraz z zastosowaniem modelu średniej ARIMA(2,2). Kolejne estymacje wykonywałem na wzór modelu GARCH(0,1):

garchspec01 <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(0,1)), mean.model=list(armaOrder=c(2,2)), distribution="std") 
garchfit01 <- ugarchfit(spec=garchspec01, data=dR.in)

Następnie wyłoniłem najlepszy model na podstawie kryteriów informacyjnych:

kryteria_info <- data.frame(infocriteria(garchfit01), infocriteria(garchfit10), infocriteria(garchfit11), infocriteria(garchfit12), infocriteria(garchfit21), infocriteria(garchfit22))
colnames(kryteria_info) <- c("GARCH(0,1)", "GARCH(1,0)", "GARCH(1,1)", "GARCH(1,2)", "GARCH(2,1)", "GARCH(2,2)")

kable(kryteria_info)

	GARCH(0,1)	GARCH(1,0)	GARCH(1,1)	GARCH(1,2)	GARCH(2,1)	GARCH(2,2)
Akaike	-3.122058	-3.191647	-3.260829	-3.259613	-3.257757	-3.256563
Bayes	-3.067023	-3.136612	-3.198914	-3.190819	-3.188962	-3.180889
Shibata	-3.122355	-3.191945	-3.261204	-3.260076	-3.258219	-3.257122
Hannan-Quinn	-3.100713	-3.170302	-3.236815	-3.232932	-3.231075	-3.227213

Wszystkie kryteria informacyjne wskazują, że najlepiej dopasowanym do danych modelem jest GARCH(1,1).

4.4 Model GARCH

W celu utwierdzenia się w przekonaniu, że w moim przypadku lepszym wyborem w modelach klasy GARCH jest rozkład t-Studenta niż rozkład normalny, przeprowadziłem analizę porównawczą modelu GARCH(1,1) z odpowiednimi rozkładami, bazującą na kryteriach informacyjnych.

garchspec11_norm <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(2,2)), distribution="norm")
garchfit11_norm <- ugarchfit(spec=garchspec11_norm, data=dR.in)

kryteria_info.garch11 <- data.frame(infocriteria(garchfit11_norm), infocriteria(garchfit11))
colnames(kryteria_info.garch11) <- c("GARCH(1,1)", "GARCH-t(1,1)")

kable(kryteria_info.garch11)

	GARCH(1,1)	GARCH-t(1,1)
Akaike	-3.172356	-3.260829
Bayes	-3.117321	-3.198914
Shibata	-3.172653	-3.261204
Hannan-Quinn	-3.151011	-3.236815

Wszystkie kryteria informacyjne przybierają wartości niższe dla modelu GARCH-t(1,1), co jednoznacznie wskazuje na przewagę rozkładu t-Studenta nad rozkładem normalnym. W związku z tym w pozostałej części analizy, wszystkie kolejne modele będą bazowały na rozkładzie t-Studenta, a dla przejrzystości wyników przedrostek “t-” będzie pomijany.

Poniżej znajduje się wydruk dopasowania modelu GARCH(1,1):

show(garchfit11)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.001422    0.001488    0.95514 0.339505
## ar1    -0.814171    0.001998 -407.51832 0.000000
## ar2    -0.998210    0.004038 -247.21800 0.000000
## ma1     0.812397    0.000181 4494.05374 0.000000
## ma2     1.010069    0.000599 1686.75922 0.000000
## omega   0.000207    0.000101    2.04354 0.040999
## alpha1  0.252967    0.091735    2.75760 0.005823
## beta1   0.725789    0.078893    9.19965 0.000000
## shape   3.959975    0.716072    5.53014 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.001422    0.001729    0.82202 0.411064
## ar1    -0.814171    0.005145 -158.24585 0.000000
## ar2    -0.998210    0.010135  -98.49532 0.000000
## ma1     0.812397    0.000396 2052.89020 0.000000
## ma2     1.010069    0.001338  755.14577 0.000000
## omega   0.000207    0.000126    1.63680 0.101673
## alpha1  0.252967    0.131359    1.92577 0.054133
## beta1   0.725789    0.104559    6.94141 0.000000
## shape   3.959975    0.748573    5.29003 0.000000
## 
## LogLikelihood : 1070.4 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.2608
## Bayes        -3.1989
## Shibata      -3.2612
## Hannan-Quinn -3.2368
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       3.073 0.07958
## Lag[2*(p+q)+(p+q)-1][11]     7.425 0.01283
## Lag[4*(p+q)+(p+q)-1][19]    14.234 0.04943
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.004241  0.9481
## Lag[2*(p+q)+(p+q)-1][5]  0.727967  0.9175
## Lag[4*(p+q)+(p+q)-1][9]  1.500974  0.9555
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2885 0.500 2.000  0.5912
## ARCH Lag[5]    1.1490 1.440 1.667  0.6893
## ARCH Lag[7]    1.3785 2.315 1.543  0.8454
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.6321
## Individual Statistics:              
## mu     0.29237
## ar1    0.04342
## ar2    0.12928
## ma1    0.02903
## ma2    0.06613
## omega  0.12692
## alpha1 0.14293
## beta1  0.20376
## shape  0.19409
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.1 2.32 2.82
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          1.58141 0.1143    
## Negative Sign Bias 1.12017 0.2631    
## Positive Sign Bias 0.03079 0.9754    
## Joint Effect       3.13729 0.3709    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     31.06      0.03978
## 2    30     55.77      0.00201
## 3    40     62.06      0.01085
## 4    50     66.43      0.04919
## 
## 
## Elapsed time : 0.445683

persistence(garchfit11)

## [1] 0.978756

Na wydruku można zaobserwować, że suma parametrów modelu jest mniejsza od 1 (funkcja persistence), parametry modelujące zmienność wariancji są dodatnie i istotne statystycznie. P-value testu Ljungi-Boxa dla kwadratów wystandaryzowanych reszt wskazuje na fakt, iż nie stwierdzono występowania autokorelacji w tych szeregach. Natomiast wartości p-value testu mnożników Lagrange’a na występowanie efektów ARCH wśród wystandaryzowanych reszt świadczą o tym, że nie można odrzucić hipotezy o niewystępowaniu efektów ARCH. Z kolei wynik Sign Bias Test mówi nam o tym, że w tym przypadku nie stwierdzono występowania efektu dźwigni.

Z uwagi na to, że kwadraty wystandaryzowanych reszt nie podlegają autokorelacji, wśród reszt nie występują efekty ARCH oraz parametry modelujące zmienność wariancji są dodatnie i istotne statystycznie (oraz mniejsze od 1), można stwierdzić, iż model ten wydaje się być dobry. Poniżej zaprezentowane są wykresy obrazujące efekty wykonanej estymacji.

Na szczególną uwagę zaslugują wykresy Empirical Density of Standirez Residuals i News Impact Curve. Pierwszy z nich przedstawia rozkład wystandaryzowanych reszt i obrazuje przewagę użytego w estymacji rozkłady t-Studenta nad rozkładem normalnym. Z kolei drugi z nich przedstawia jaki wpływ na osiąganą stopę zwrotu mają wpływ negatywne i pozytywne informacje z rynku. W standardowych modelach GARCH, tzw. efekt dźwigni nie jest uwzględniony, dzięki czeu możemy zaobserwować, iż wpływ przeciwstawnych informacji na wariancję zwrotów jest tutaj jednakowy.

4.5 Model TGARCH

Bazując na osiągniętych wyżej wynikach dla modelu GARCH(1,1), przeprowadziłem odpowiednią analizę dla modelu TGARCH(1,1).

tgarchspec11 <- ugarchspec(variance.model=list(model="fGARCH", submodel="TGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(2,2)), distribution="std") 
tgarchfit11 <- ugarchfit(spec=tgarchspec11, data=dR.in) 

show(tgarchfit11)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : fGARCH(1,1)
## fGARCH Sub-Model : TGARCH
## Mean Model   : ARFIMA(2,0,2)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.002571    0.001458     1.7638 0.077767
## ar1    -0.815156    0.002486  -327.9279 0.000000
## ar2    -0.997584    0.000551 -1809.7403 0.000000
## ma1     0.812342    0.000275  2952.7317 0.000000
## ma2     1.009945    0.000508  1989.5465 0.000000
## omega   0.002376    0.000884     2.6868 0.007213
## alpha1  0.103663    0.025277     4.1010 0.000041
## beta1   0.875704    0.031538    27.7670 0.000000
## eta11  -0.672320    0.217624    -3.0894 0.002006
## shape   4.228062    0.752495     5.6187 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.002571    0.001618    1.58952 0.111943
## ar1    -0.815156    0.006503 -125.34722 0.000000
## ar2    -0.997584    0.005090 -196.00214 0.000000
## ma1     0.812342    0.000473 1716.71439 0.000000
## ma2     1.009945    0.001432  705.07988 0.000000
## omega   0.002376    0.001227    1.93733 0.052705
## alpha1  0.103663    0.121749    0.85145 0.394522
## beta1   0.875704    0.081261   10.77648 0.000000
## eta11  -0.672320    0.992860   -0.67715 0.498307
## shape   4.228062    0.855810    4.94042 0.000001
## 
## LogLikelihood : 1075.217 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.2726
## Bayes        -3.2038
## Shibata      -3.2730
## Hannan-Quinn -3.2459
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.6186  0.4316
## Lag[2*(p+q)+(p+q)-1][11]    5.0428  0.9502
## Lag[4*(p+q)+(p+q)-1][19]   10.8044  0.3417
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.6108  0.4345
## Lag[2*(p+q)+(p+q)-1][5]    0.9530  0.8702
## Lag[4*(p+q)+(p+q)-1][9]    1.8446  0.9227
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2126 0.500 2.000  0.6448
## ARCH Lag[5]    0.5330 1.440 1.667  0.8738
## ARCH Lag[7]    1.1245 2.315 1.543  0.8925
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.6522
## Individual Statistics:              
## mu     0.24459
## ar1    0.03425
## ar2    0.14413
## ma1    0.04009
## ma2    0.05906
## omega  0.02749
## alpha1 0.03214
## beta1  0.03873
## eta11  0.06727
## shape  0.18773
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.29 2.54 3.05
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias            1.026 0.30525    
## Negative Sign Bias   2.505 0.01251  **
## Positive Sign Bias   0.188 0.85095    
## Joint Effect         6.365 0.09514   *
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     32.53      0.02720
## 2    30     41.86      0.05784
## 3    40     51.86      0.08157
## 4    50     68.74      0.03281
## 
## 
## Elapsed time : 1.259575

persistence(tgarchfit11)

## [1] 0.9498691

Suma parametrów modelu jest mniejsza od 1. P-value Weighted Ljung-Box Test on Standardized Squared Residuals mówi nam o tym, że nie mamy podstaw aby odrzucić \(H_{0}\) jakoby kwadraty wystandaryzowanych reszt były białym szumem. Z kolei Weighted ARCH LM Tests wskazuje na brak występowania efektów ARCH. Wyniki Sign Bias Test świadczą o występowaniu efektu dźwigni: odrzucamy \(H_0\) o występowaniu nieistotnych negatywnych reakcji na szoki. Model wydaje się być poprawny. Opisane zależności można zweryfikować na poniższych wykresach.

W tym wypadku należy zwrócić uwagę na wykres News Impact Curve. Widać wyraźnie, że efekt dźwigni został uwzględniony, co charakteryzuje się odmiennym wpływem informacji pozytywnych i negatywnych na wariancję zwrotów z portfela.

4.6 Model EGARCH

Analogiczną analizę przeprowadziłem dla modelu EGARCH(1,1).

egarchspec11 <- ugarchspec(variance.model=list(model="eGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(2,2)), distribution="std") 
egarchfit11 <- ugarchfit(spec=egarchspec11, data=dR.in) 

show(egarchfit11)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : eGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu       0.00254    0.001663     1.5277 0.126590
## ar1     -0.81507    0.000771 -1056.9561 0.000000
## ar2     -0.99724    0.002521  -395.5130 0.000000
## ma1      0.81234    0.000223  3642.1680 0.000000
## ma2      1.00998    0.000154  6575.4755 0.000000
## omega   -0.25687    0.028172    -9.1178 0.000000
## alpha1   0.12199    0.027085     4.5041 0.000007
## beta1    0.95860    0.004442   215.8161 0.000000
## gamma1   0.13822    0.019991     6.9144 0.000000
## shape    4.16973    0.736574     5.6610 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu       0.00254    0.002387    1.0639 0.287358
## ar1     -0.81507    0.001809 -450.5480 0.000000
## ar2     -0.99724    0.004266 -233.7766 0.000000
## ma1      0.81234    0.000267 3039.1683 0.000000
## ma2      1.00998    0.000210 4810.2621 0.000000
## omega   -0.25687    0.039023   -6.5823 0.000000
## alpha1   0.12199    0.029474    4.1391 0.000035
## beta1    0.95860    0.005634  170.1486 0.000000
## gamma1   0.13822    0.036509    3.7860 0.000153
## shape    4.16973    0.798807    5.2199 0.000000
## 
## LogLikelihood : 1075.306 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.2728
## Bayes        -3.2040
## Shibata      -3.2733
## Hannan-Quinn -3.2461
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.2362  0.6270
## Lag[2*(p+q)+(p+q)-1][11]    5.6193  0.7296
## Lag[4*(p+q)+(p+q)-1][19]   11.7077  0.2243
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.034  0.3091
## Lag[2*(p+q)+(p+q)-1][5]     1.186  0.8163
## Lag[4*(p+q)+(p+q)-1][9]     1.913  0.9151
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.02774 0.500 2.000  0.8677
## ARCH Lag[5]   0.15232 1.440 1.667  0.9763
## ARCH Lag[7]   0.46622 2.315 1.543  0.9813
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.8528
## Individual Statistics:              
## mu     0.23050
## ar1    0.03608
## ar2    0.15352
## ma1    0.03629
## ma2    0.05944
## omega  0.05354
## alpha1 0.07150
## beta1  0.04486
## gamma1 0.25033
## shape  0.31647
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.29 2.54 3.05
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           1.1450 0.252622    
## Negative Sign Bias  2.7803 0.005588 ***
## Positive Sign Bias  0.4192 0.675201    
## Joint Effect        7.9237 0.047616  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     32.41      0.02809
## 2    30     42.13      0.05462
## 3    40     56.77      0.03275
## 4    50     62.59      0.09184
## 
## 
## Elapsed time : 1.053798

persistence(egarchfit11)

## [1] 0.9585967

W tym przypadku podobnie jak w poprzendich modelach, suma parametrów jest mniejsza od 1. P-value Weighted Ljung-Box Test on Standardized Squared Residuals mówi nam o tym, że kwadraty wystandaryzowanych reszt nie podlegają autokorelacji. Z kolei Weighted ARCH LM Tests wskazuje na brak występowania efektów ARCH. Natomiast wyniki Sign Bias Test świadczą o występowaniu efektu dźwigni. Model ten również wydaje się być poprawny. Poniżej przedstawiono wykresy, obrazujące opisane wyniki.

Tak jak w poprzednim przypadku (TGARCH) News Impact Curve wskazuje na wsytępowanie efektu dźwigni, lecz charakteryzuje się nieco odmiennym przebiegiem.

4.7 Wybór najlepszego modelu

Dla oszacowanych powyżej modeli, wyodrębniłem poszczególne kryteria informacyjne, w celu oceny, który z nich najlepiej dopasowuje się do danych.

kryteria_info2 <- data.frame(infocriteria(garchfit11), infocriteria(tgarchfit11), infocriteria(egarchfit11))
colnames(kryteria_info2) <- c("GARCH(1,1)", "TGARCH(1,1)", "EGARCH(1,1)")

kable(kryteria_info2)

	GARCH(1,1)	TGARCH(1,1)	EGARCH(1,1)
Akaike	-3.260829	-3.272556	-3.272830
Bayes	-3.198914	-3.203762	-3.204036
Shibata	-3.261204	-3.273019	-3.273293
Hannan-Quinn	-3.236815	-3.245874	-3.246148

Według wszystkich kryteriów informacyjnych najlepszym modelem jest EGARCH(1,1).

5 Analiza warunkowej wariancji i VaR w okresie in-sample

5.1 Warunkowa i bezwarunkowa wariancja

W poniższej tabeli przedstawione są wartości bezwarunkowej wariancji dla poszczególnych, oszacowanych wcześniej modeli.

unc_variance <- cbind(uncvariance(garchfit11), uncvariance(tgarchfit11), uncvariance(egarchfit11))
colnames(unc_variance) <- c("GARCH(1,1)", "TGARCH(1,1)", "EGARCH(1,1)")

kable(unc_variance)

GARCH(1,1)	TGARCH(1,1)	EGARCH(1,1)
0.0097219	0.0022468	0.0020214

Najmniejszą bezwarunkową wariancję przyjmuje model EGARCH(1,1), który został wyłoniony wcześniej jako model najlepiej opisujący posiadane dane. Zbliżoną, również niską wartością bezwarunkowej wariancji cechuje się model TGARCH(1,1). Natomiast model wyjściowy GARCH(1,1) charakteryzuje się największą bezwarunkową wariancją.

Następnie przystąpiłem do wyznaczenia warunkowych wariancji w rozpatrywanych modelach. Wykresy warunkowej wariancji, wraz z naniesioną na nie wartością bezwarunkowej wariancji oraz użyty w tym celu schemat kodu, znajdują się poniżej.

var_garch11 <- xts(cbind(data.frame(garchfit11@fit$var), uncvariance(garchfit11)), rev(insample.df$Date))
colnames(var_garch11) <- c("conditional.variance", "unconditional.variance")

dygraph(var_garch11, main="Wariancje dla modelu GARCH(1,1): in-sample") %>% 
  dyRangeSelector(fillColor = "white") %>%
  dyAxis("y", label = "Variance") %>%
  dySeries("conditional.variance", color = "green") %>%
  dySeries("unconditional.variance", color="red", strokePattern="dashed")

Analizując powyższe wykresy można zaobserwować, że zdecydowanie największe szoki warunkowej wariancji pojawiają się na wykresie obrazującym model GARCH(1,1). Z kolei warunkowe wariancje dla dwóch pozostałych modeli charakteryzują się zbliżonym przebiegiem.

5.2 VaR

Jednym z kluczowych zastosowań modeli klasy GARCH jest szacowanie prawdopodobieństwa wystąpienia wysokich strat. Taką wysoką stratę zwykło się określać jako VaR (ang. Value at Risk). Jest to strata jaka zostaje przekroczona z danym prawdopodobieństwem równym \(\alpha\) (zazwyczaj \(\alpha = 0.01\)). Rozkład VaR zmienia się w czasie w związku ze zmianami warunkowej wariancji.

Wyznaczyłem wartości VaR w okresie in-sample dla rozpatrywanych modeli według schematu postępowania:

zt_garch11 <- residuals(garchfit11, standardize=TRUE)
sigma_garch11 <- sigma(garchfit11)
quantile_garch11 <- quantile(zt_garch11, 0.01, na.rm=TRUE)  
VaR_garch11 <- quantile_garch11 * sigma_garch11
colnames(VaR_garch11) <- c("VaR")

dRVaR_garch11 <- cbind(dR.in, VaR_garch11)
dygraph(dRVaR_garch11, main="VaR dla modelu GARCH(1,1): in-sample") %>% 
  dyRangeSelector(fillColor = "white") %>%
  dyAxis("y", label = "Log return") %>%
  dySeries("daily.returns", color = "#663300") %>%
  dySeries("VaR", color = "red")

Aby lepiej zrozumieć otrzymane wyniki, wyznaczyłem średnie wartości VaR dla rozpatrywanych modeli.

mean_VaR <- cbind(mean(VaR_garch11), mean(VaR_tgarch11), mean(VaR_egarch11))
colnames(mean_VaR) <- c("GARCH(1,1)", "TGARCH(1,1)", "EGARCH(1,1)")
kable(mean_VaR)

GARCH(1,1)	TGARCH(1,1)	EGARCH(1,1)
-0.1215381	-0.1180136	-0.1146738

Największą średnią wartość narażoną na ryzyko odnotowuje się dla modelu GARCH(1,1), natomiast najmniejszą dla modelu EGARCH(1,1) (bierzemy pod uwagę wartości bezwzględne zamieszczonych w powyższej tabeli wielkości). W celu porównania przebiegów VaR dla analizowanych modeli stworzyłem ich wspólny wykres.

VaR <- cbind(VaR_garch11, VaR_tgarch11, VaR_egarch11)
colnames(VaR) <- c("VaR.garch11", "VaR.tgarch11", "VaR.egarch11")
dygraph(cbind(dR.in, VaR), main="VaR dla poszczególnych modeli: in-sample") %>% 
  dyRangeSelector(fillColor = "white") %>%
  dyAxis("y", label = "Log return") %>%
  dySeries("daily.returns", color="brown") %>%
  dySeries("VaR.garch11", color = "red") %>%
  dySeries("VaR.tgarch11", color = "magenta") %>%
  dySeries("VaR.egarch11", color = "purple")

Ważną kwestią w ocenie efektywności wyznaczonych szeregów VaR jest sprawdzenie odsetka przypadków, w których straty przekroczyły zakładany poziom wartości narażonej na ryzyko.

VaR_garch11.przekrocz <- sum(dR.in < VaR_garch11) / length(VaR_garch11)
VaR_tgarch11.przekrocz <- sum(dR.in < VaR_tgarch11) / length(VaR_tgarch11)
VaR_egarch11.przekrocz <- sum(dR.in < VaR_egarch11) / length(VaR_egarch11)
VaR.przekrocz <- cbind(VaR_garch11.przekrocz, VaR_tgarch11.przekrocz, VaR_egarch11.przekrocz)
colnames(VaR.przekrocz) <- c("GARCH(1,1)", "TGARCH(1,1)", "EGARCH(1,1)")
rownames(VaR.przekrocz) <- "VaR przekroczone"
kable(VaR.przekrocz)

	GARCH(1,1)	TGARCH(1,1)	EGARCH(1,1)
VaR przekroczone	0.0107527	0.0107527	0.0092166

VaR zostało przekroczone najmniejszą ilość razy dla modelu EGARCH(1,1).

6 Analiza warunkowej wariancji i VaR w okresie out-of-sample

Kolejnym krokiem było wykonanie podobnej analizy co poprzednio, lecz tym razem w okresie out-of-sample. Próba out-of-sample obejmowała okres 2 miesięcy (22.05.2017 - 22.07.2017). W celu estymacji modeli zastosowałem funkcję ugarchroll z pakietu rugarch według wzorca:

roll_garch11 <- ugarchroll(spec=garchspec11, data=dR, n.ahead=1, forecast.length=62, refit.every=2, refit.window='recursive', window.size=NULL, solver="hybrid", calculate.VaR=TRUE, VaR.alpha=c(0.01, 0.05), keep.coef=TRUE)

6.1 Warunkowa wariancja

Wykresy przedstawiające warunkową wariancję w okresie out-of-sample wykonałem na wzór poniższego kodu:

var_garch11.out <- data.frame(roll_garch11@forecast$density[,2] ^ 2)
var_garch11.out <- xts(var_garch11.out, outsample.df$Date)
colnames(var_garch11.out) <- "var_garch11.out"

dygraph(var_garch11.out, main="Warunkowa wariancja dla modelu GARCH(1,1): out-of-sample") %>% 
  dyRangeSelector(fillColor = "white") %>%
  dyAxis("y", label = "Variance") %>%
  dySeries(color = "blue")

Warunkowe wariancje dla analizowanych modeli w okresie out-of-sample różnią się od siebie. W przypadku modelu TGARCH(1,1) i EGARCH(1,1) różnice te nie są wielkie, natomiast jeżeli porównamy je z modelem GARCH(1,1), uda nam się zaobserwować znacznie odbiegające od siebie wartości. Z podobną sytuacją mieliśmy do czynienia dla warunkowej wariancji w okresie in-sample: wyniki modeli EGARCH(1,1) i TGARCH(1,1) są zbliżone, ze względu na podobieństwa w założeniach. Aby lepiej przedstawić różnice w warunkowej wariancji, sporządziłem również wykres wspólny dla wszystkich trzech modeli.

var.out <- cbind(var_garch11.out, var_tgarch11.out, var_egarch11.out)
dygraph(var.out, main="Warunkowa wariancja dla wszystkich modeli: out-of-sample") %>% 
  dyRangeSelector(fillColor = "white") %>%
  dyAxis("y", label = "Variance") %>%
  dySeries("var_garch11.out", color = "blue", label="GARCH(1,1)") %>%
  dySeries("var_tgarch11.out", color = "#6495ED", label="TGARCH(1,1)") %>%
  dySeries("var_egarch11.out", color = "#48D1CC", label="EGARCH(1,1)")

6.2 VaR

Sposób wyznaczenia wartości narażonej na ryzyko w okresie out-of-sample:

VaR_garch11.out <- data.frame(roll_garch11@forecast$VaR[,1])
VaR_garch11.out <- xts(VaR_garch11.out, outsample.df$Date)
colnames(VaR_garch11.out) <- "VaR_garch11.out"

Aby poprawnie przedstawić VaR należało również obliczyć zrealizowane stopy zwrotu dla próby out-of-sample:

dR.out <- data.frame(roll_garch11@forecast$density[,6])
dR.out <- xts(dR.out, outsample.df$Date)
colnames(dR.out) <- "daily.return.out"

Poniżej wykresy przedstawiające relację zrealizowany stóp zwrotu oraz wartości narażonej na ryzyko.

Przebiegi wartości narażonej na ryzyko dla rozpatrywanych modeli w okresie out-of-sample różnią się między sobą. Ponownie największe odchylenia zaobserwować można pomiędzy modelem GARCH(1,1), a modelami uwzględniającymi efekt dźwigni.

mean_VaR.out <- cbind(mean(VaR_garch11.out), mean(VaR_tgarch11.out), mean(VaR_egarch11.out))
colnames(mean_VaR.out) <- c("GARCH(1,1)", "TGARCH(1,1)", "EGARCH(1,1)")
kable(mean_VaR.out)

GARCH(1,1)	TGARCH(1,1)	EGARCH(1,1)
-0.2108087	-0.191641	-0.1953751

Analogicznie jak podczas analizy próby in-sample, także tym razem wyznaczyłem odsetek strat, które przekroczyły wartość narażoną na ryzyko.

VaR_garch11.out.przekrocz <- sum(dR.out < VaR_garch11.out) / length(VaR_garch11.out)
VaR_tgarch11.out.przekrocz <- sum(dR.out < VaR_tgarch11.out) / length(VaR_tgarch11.out)
VaR_egarch11.out.przekrocz <- sum(dR.out < VaR_egarch11.out) / length(VaR_egarch11.out)
VaR.out.przekrocz <- cbind(VaR_garch11.out.przekrocz, VaR_tgarch11.out.przekrocz, VaR_egarch11.out.przekrocz)
colnames(VaR.out.przekrocz) <- c("GARCH(1,1)", "TGARCH(1,1)", "EGARCH(1,1)")
rownames(VaR.out.przekrocz) <- "VaR przekroczone"

kable(VaR.out.przekrocz) 

	GARCH(1,1)	TGARCH(1,1)	EGARCH(1,1)
VaR przekroczone	0.016129	0.016129	0

Dla modelu EGARCH(1,1) VaR nie zostało przekroczone ani razu.

7 Analiza wrażliwości

Przeprowadzona poniżej analiza wrażliwości ma na celu ukazanie jak przyjęcie początkowych założeń analizy, takich jak okno próby, wpływa na estymację parametrów i końcowe wyniki. W związku z tym zdecydowałem się na oszacowanie wyłonionego wczęsniej jako najlepszy, modelu EGARCH(1,1), bazując na kilku różnych wariantach długości próby out-of-sample: 5 miesięcy, 3 miesiące, 1 miesiąc, 2 tygodnie. Otrzymane wyniki porównałem z analizowanym uprzednio modelem EGARCH(1,1) na próbie out-of-sample o długości 2 miesięcy. Całość analizy przeprowadziłem analogicznie jak w poprzednim rozdziale, dlatego też pozwolę sobie przedstawić jedynie osiągnięte wyniki.

Pierwszym rozpatrywanym oknem długości próby out-of-sample było 5 miesięcy (22.02.2017 - 22.07.2017).

Następnie rozpatrzyłem okres out-of-sample o długości 3 miesięcy (22.04.2017 - 22.07-2017).

Kolejną analizowaną długością próby out-of-sample był 1 miesiąc (22.06.2017 - 22.07-2017).

Ostatnią rozpatrywaną długością próby out-of-sample były 2 tygodnie (09.07.2017 - 22.07-2017).

W celu zestawienia ze sobą osiągniętych powyżej wyników, porównałem wartości średnie VaR w modelu EGARCH(1,1) dla poszczególnych długości próby out-of-sample.

5 mies.	3 mies.	2 mies.	1 mies.	2 tyg.
-0.2033518	-0.2165026	-0.1953751	-0.1752499	-0.2288379

Z uwagi na fakt, że wartość narażona na ryzyko silnie zależy od warunkowej wariancji zwrotów w danym okresie, poszczególne średnie VaR różnią się od siebie, lecz nie są to znaczące różnice. Największą można zaobserwować pomiędzy średnimi VaR dla 1 miesiąca i 2 tygodni. Pozostałe VaR cechują się zbliżonymi wartościami, osculującymi wokół 0.2, czyli wokół VaR dla pierwotnie rozpatrywanej dwumiesięcznej długości próby out-of-sample (0.1954).

Wyznaczyłem również odsetek przypadków, w których straty przewyższyły wartość narażoną na ryzyko.

	5 mies.	3 mies.	2 mies.	1 mies.	2 tyg.
VaR przekroczone	0	0	0	0	0

Jak widać, dla każdej rozpatrywanej długości próby out-of-sample straty ani razu nie przekroczyły VaR. Z wykresów zaobserwować można jednak, że w jednym przypadku osiągnęły one wielkość równą VaR (10.07.2017 r. straty z portfela i VaR zanotowały wartość 0.14).

8 Podsumowanie

W przeprowadzonej analizie porównałem oszacowania funkcji warunkowej wariancji oraz wartości narażonej na ryzyko w dwóch rozszerzeniach modeli klasy GARCH: TGARCH i EGARCH. Z uwagi na fakt, że wśród modeli GARCH o różnej kombinacji parametrów najlepszym okazał się model GARCH(1,1) z rozkładem t-Studenta, to do niego odnosiłem wyniki osiągnięte za pomocą dwóch rozszerzeń o analogicznych parametrach i rozkładzie: TGARCH(1,1) oraz EGARCH(1,1). Spośród wszystkich trzech modeli najlepiej dopasowanym do danych okazał się model EGARCH(1,1), lecz osiągane przez niego przebiegi, zarówno warunkowej wariancji jak i wartości narażonej na ryzyko, były zbliżone do modelu TGARCH(1,1). Wiąże się to z faktem, iż modele te uwzględniają możliwość asymetrycznej reakcji wariancji zwrotów na dodatnie i ujemne odchylenia losowe, co jest założeniem zgodnym z funkcjonowaniem rynków finansowych. Nieco różniące się wyniki, choć znacząco nie odbiegające, można było otrzymać za pomocą estymacji modelu wyjściowego GARCH(1,1).
W przeprowadzonej na końcu badania analizie wrażliwości dla najlepiej dopasowanego modelu, EGARCH(1,1), zbadałem jak zmiana okna próby może wpłynąć na końcowe wyniki. Otrzymane wyniki nie wykazały znaczących różnic, odnosząc się do wyjściowego, dwumiesięcznego okna estymacji.