Analiza porównawcza ryzyka w modelach klasy GARCH dla portfela wybranych aktywów
Wojciech Snarski
25 czerwca 2017
Streszczenie
Celem poniższej pracy jest analiza ryzyka równoważonego portfela aktywów. Modelowanie ryzyka dokonano za pomocą oszacowania funkcji warunkowej wariancji z wykorzystaniem modeli klasy GARCH. Analizę porównawczą przeprowadzono zestawiając ze sobą standardowy model GARCH oraz dwa rozszerzenia tego modelu: EGARCH i GARCH-t. Oszacowano warunkową wariancję oraz wartość narażoną na ryzyko w okresach in-sample i out-of-sample. Dokonano również analizy wrażliwości uzyskanych oszacowań w zależności od przyjętych początkowych parametrów.
Budowa portfela
Na początku zaimportowano wszystkie niezbędne pakiety:
library(xts)
library(fBasics)
library(tseries)
library(car)
library(FinTS)
library(fGarch)
library(rugarch) Do przeprowadzenia analizy wybrano 6 aktywów, których notowania zaczerpnięto z serwisu https://stooq.com/. Wybrano 2 kursy walut:
-Kurs kryptowaluty Bitcoin wyrażony w polskich złotych
-Kurs korony norweskiej wyrażony w PLN
Notowania 3 spółek giełdowych:
-Netia S.A.
-CD Projekt S.A.
-Mostostal Zabrze S.A.
Oraz:
-Kontrakt terminowy 3-miesięcznej ceny cynku.
W zaimportowanych danych wybrano jedynie ceny zamknięcia i wyznaczono logarytmiczne ciągłe stopy zwrotu dla każdego z aktywów. Poniżej procedura dla jednego aktywa:
url_btc <- "https://stooq.com/q/d/l/?s=btcpln&i=d"
bitcoin <- read.csv(url_btc,header = TRUE,
sep = ",", dec = ".", stringsAsFactors = F)
bitcoin$Date <- as.Date(bitcoin$Date)
bitcoin <- bitcoin[, c("Date", "Close")]
colnames(bitcoin) <- c("Date", "bitcoin")
bitcoin$rbitcoin <- diff.xts(log(bitcoin$bitcoin)) Następnie wyznaczono początkową datę dla analizowanego portfela jako wartość najwyżsżą z początkowych dat każdego z szeregów.
s_date<-max(min(mzab$Date), min(nok$Date), min(netia$Date),
min(bitcoin$Date), min(cdpr$Date),min(zinc$Date)) Aby wyniki nie były zaburzone poprzez pojawianie się nowych obserwacji w szeregach, przyjęto ograniczenie górne badanego zakresu.
f_date <- "2017-06-16"Ograniczono każdy z szeregów za pomocą wyznaczonych dat w taki sposób, aby każdy z nich miał taką samą długość. Dla każdego z szeregów zastosowano formułę:
bitcoin <- bitcoin[bitcoin$Date <= as.Date(f_date), ]
bitcoin <- bitcoin[bitcoin$Date > as.Date(s_date), ] Należy się tutaj mała uwaga. Notowania dla szrergów netia, cdpr i mzab zawierają mniej rekordów, ponieważ nie występują notowania w takie dni jak 11 listopada czy 6 stycznia, czyli polskie dni świąteczne. Z kolei dla pozostałych szeregów dla tych dni występują obserwacje. W przypadku tego badania zdecydowano się pominąć te dni w analizie, to znaczy przy budowie portfela uwzględniono jedynie te dni, dla których każdy z szeregów miał wartość niezerową.
Połączono szeregi do wspólnej tablicy danych.
merg1<-merge(bitcoin,netia, by="Date")
merg2<-merge(nok,cdpr, by="Date")
merg3<-merge(merg1,merg2, by="Date")
merg4<-merge(merg3,mzab, by="Date")
x1<- merge(merg4,zinc, by="Date")Dodano numer obserwacji.
x1$obs<-1:length(x1$bitcoin)Portfel zbudowano jako sumę zwrotów z aktywów z odpowiednimi wagami. W przypadku portfela równoważonego każde aktywo ma taką samą wagę, zatem zastosowano następującą formułę:
x1$r<-((x1$rbitcoin+x1$rnetia+x1$rnok+x1$rcdpr+x1$rmzab+x1$rzinc)/6) Zatem badany okres obejmował obserwacje od 2010-07-20 do 2017-06-16
Analiza graficzna portfela
Poniżej przedstawiono wykresy notowań dla każdego z aktywów osobno. Szybkie spostrzeżenia są takie, że aktywa raczej nie są powiązane podobnymi mechanizami kształtowania cen. Ze wstępnej analizy wykresów można stwierdzić, że zmienne te nie są stacjonarne - stałość średniej i wariancji nie jest zachowana.
plot(x1$Date, x1$bitcoin, type = "l",
main = "Dzienne ceny zamknięcia kursu bitcoin/PLN",
xlab = "", ylab = "")
Poniżej wykres dziennych zwrotów z portfela. Wariancja wydaje się nie być stała w czasie. Występują okresy wyższych wahań (np. 2011-2012) i niższych (np. połowa 2016 do połowy 2017)
plot(x1$Date, x1$r, type = "l", col = "red1",
main = "Dzienne zwroty z portfolio",
xlab = "", ylab = "zwroty")
abline(h = 0, col = "black")
Badanie normalności zwrotów
Wygrenerowano statystyki opisowe zwrotów z portfela. Średnia rozkładu wyniosła 0.001, co jest wartością zbliżoną do zera. Kurtoza miała wartość 4.987, co różni się od wartości 3 dla rozkładu normalnego.
Poniżej przedstawiono histogram z naniesioną gęstością rozkładu normalnego z parametrami oszacowanymi na podstawie pełnej próbki zwrotów. Z analizy wykresu można stwierdzić, że rozkład wygląda na leptokurtyczny: zbyt duże skupienie wokół średniej i zbyt grube ogony.
hist(x1$r, prob = T, breaks = 120, main="Histogram")
curve(dnorm(x, mean = mean(x1$r, na.rm = T),
sd = sd(x1$r, na.rm = T)),
col = "darkblue", lwd = 2, add = TRUE)
Za pomocą statystyki Jarque-Bery można stwierdzić, że hipoteza o normalności zwrotów jest silnie odrzucana.
jarque.bera.test(na.omit(x1$r))##
## Jarque Bera Test
##
## data: na.omit(x1$r)
## X-squared = 1770.4, df = 2, p-value < 2.2e-16Wykres funkcji autokrelacji dla zwrotów z portfela wskazuje na występowanie istotnych opóźnień. Są to opóźnienia: 1,3,9,11,12,16,18 i 23. Oznacza to, że wartości zmiennej zależą istotnie od jej niektórych przeszłych wartości.
acf(x1$r, lag.max = 48, na.action = na.pass,
ylim = c(-0.05, 0.15),
col = "blue3", lwd = 6,
main = "ACF zwrotów z portfolio")
Wykres ACF kwadratów zwrotów wskazuje na występowanie większości istotnych opóźnień. 
Za pomocą statystyki Durbina-Watsona stwierdzono występowanie istotnych (na poziomie istotności 5%) opóźnień 1 i 3.
durbinWatsonTest(lm(formula = x1$r ~ 1),max.lag = 5)## lag Autocorrelation D-W Statistic p-value
## 1 0.06621972 1.867369 0.006
## 2 0.04729435 1.901896 0.038
## 3 0.06059950 1.858351 0.002
## 4 0.02623190 1.921314 0.120
## 5 0.02433288 1.923832 0.134
## Alternative hypothesis: rho[lag] != 0Wynik testu na występowanie efektów ARCH wśród zwrotów sugeruje odrzucenie hipotezy o braku wystepowania efektów ARCH.
ArchTest(x1$r, lags = 5)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: x1$r
## Chi-squared = 128.68, df = 5, p-value < 2.2e-16Wybór parametrów modeli
Standardowy model GARCH
Przy wyborze parametrów modelu należy kierować się następującymi wytycznymi:
-warunek stacjonarności w wariancji procesów GARCH to: \[ \sum_{i=1}^{q}\alpha_i+ \sum_{i=1}^{p}\beta_i<1 \]
-aby wariancja była dodatnia, oszacowania parametrów powinny być nieujemne
-kwadraty wystandaryzowanych reszt nie mogą podlegać autokorelacji, ponieważ w p.p. występują wśród nich efekty ARCH
-wśród reszt nie mogą występować efekty ARCH
-parametry modelu powinny być istotne
-model powinień być oszczędny w parametry
-dopasowanie modeli porównuje się za pomocą kryteriów informacyjnych, czyli wybiera się model z najniższymi wartościami statystyk AIC i BIC
Na początku sprawdzono model GARCH(1,1):
garch11 <- garchFit(formula = ~ garch(1, 1),
data = na.omit(x1$r),
include.mean = F,
cond.dist = "norm",
trace = FALSE)
summary(garch11)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = na.omit(x1$r), cond.dist = "norm",
## include.mean = F, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x000000001d59a430>
## [data = na.omit(x1$r)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## omega alpha1 beta1
## 1.1571e-05 1.6376e-01 7.8247e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 1.157e-05 2.780e-06 4.161 3.16e-05 ***
## alpha1 1.638e-01 2.742e-02 5.972 2.34e-09 ***
## beta1 7.825e-01 3.329e-02 23.507 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 5043.514 normalized: 2.970267
##
## Description:
## Mon Jun 26 12:12:08 2017 by user: Administrator
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 1037.289 0
## Shapiro-Wilk Test R W 0.9657199 0
## Ljung-Box Test R Q(10) 34.08877 0.000178411
## Ljung-Box Test R Q(15) 48.14791 2.407061e-05
## Ljung-Box Test R Q(20) 56.04206 2.864801e-05
## Ljung-Box Test R^2 Q(10) 4.03268 0.9458607
## Ljung-Box Test R^2 Q(15) 5.77326 0.9833098
## Ljung-Box Test R^2 Q(20) 10.96093 0.9472295
## LM Arch Test R TR^2 4.127944 0.9810154
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -5.937001 -5.927394 -5.937007 -5.933444Suma parametrów modelu jest mniejsza od 1, a wszystkie parametry są dodatnie i istotne. P-value testu Ljung-Boxa dla kwadratów wystandaryzowanych reszt wskazuje na brak podstaw do odrzucenia hipotezy, że kwadraty wystandaryzowanych reszt są białym szumem. P-value testu LM ARCH wskazuje na brak występowania efektów ARCH wśród reszt modelu. Model jest bardzo oszczędny w parametry, a wartość AIC wynosi -5.937. Model ten spełnia wszystkie warunki i wydaje się być dobry.
Następnie sprawdzono model GARCH(1,2).
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 2), data = na.omit(x1$r), cond.dist = "norm",
## include.mean = F, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 2)
## <environment: 0x0000000010a3a0b8>
## [data = na.omit(x1$r)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## omega alpha1 beta1 beta2
## 0.00001316 0.20100309 0.32462467 0.41285283
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 1.316e-05 3.426e-06 3.842 0.000122 ***
## alpha1 2.010e-01 3.396e-02 5.918 3.26e-09 ***
## beta1 3.246e-01 1.085e-01 2.993 0.002767 **
## beta2 4.129e-01 1.042e-01 3.963 7.39e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 5046.682 normalized: 2.972133
##
## Description:
## Mon Jun 26 12:12:09 2017 by user: Administrator
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 961.663 0
## Shapiro-Wilk Test R W 0.9671943 0
## Ljung-Box Test R Q(10) 34.37014 0.0001598376
## Ljung-Box Test R Q(15) 49.31991 1.554341e-05
## Ljung-Box Test R Q(20) 56.64094 2.326448e-05
## Ljung-Box Test R^2 Q(10) 3.315462 0.972996
## Ljung-Box Test R^2 Q(15) 5.234923 0.9899436
## Ljung-Box Test R^2 Q(20) 9.503367 0.9763099
## LM Arch Test R TR^2 3.71358 0.9880679
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -5.939555 -5.926747 -5.939566 -5.934813Model ten również spełnia wszystkie wymagania, a wartość AIC jest nieznacznie wyższa niż w modelu GARCH(1,1). Zawiera on więcej parametrów, zatem nie ma podstaw do twierdzenia, że jest on lepszy niż model (1,1).
Sprawdzono model GARCH(2,2):
W modelu tym parametry \(\alpha_2\) i \(\beta_1\) okazały się nieistone statystycznie. Również w modelu GARCH(2,1) parametr \(\alpha_2\) był nieistotny.
W związku z tym najlepiej dopasowanym do analizowanych danych modelem z klasy GARCH jest model GARCH(1,1), który w literaturze empirycznej jest dosyć popularny, m.in. dzięki swojej prostocie. Uzyskano następujące parametry modelu:
omega=0
alpha1=0.164
beta1=0.782
Model EGARCH
W przypadku modeli EGARCH modelowany jest logarytm wariancji. Model umożliwia zidentyfikowanie asymetrycznej reakcji wariancji na szoki. W przypadku modeli EGARCH oszacowania parametrów nie muszą być dodatnie.
Oszacowano model EGARCH(1,1):
spec = ugarchspec(variance.model = list(model ="eGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = F),
distribution.model = "norm")
egarch11 <- ugarchfit( spec = spec, data = na.omit(x1$r))
egarch11##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega -0.528124 0.124968 -4.22607 0.000024
## alpha1 -0.003222 0.014374 -0.22416 0.822630
## beta1 0.937914 0.014326 65.46714 0.000000
## gamma1 0.293648 0.035332 8.31104 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega -0.528124 0.231447 -2.28184 0.022499
## alpha1 -0.003222 0.030949 -0.10411 0.917079
## beta1 0.937914 0.026783 35.01894 0.000000
## gamma1 0.293648 0.070142 4.18645 0.000028
##
## LogLikelihood : 5044.358
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.9368
## Bayes -5.9240
## Shibata -5.9368
## Hannan-Quinn -5.9321
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 9.214 0.0024013
## Lag[2*(p+q)+(p+q)-1][2] 11.114 0.0009708
## Lag[4*(p+q)+(p+q)-1][5] 14.810 0.0004998
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5143 0.4733
## Lag[2*(p+q)+(p+q)-1][5] 1.4555 0.7509
## Lag[4*(p+q)+(p+q)-1][9] 2.3497 0.8594
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.094 0.500 2.000 0.2956
## ARCH Lag[5] 1.151 1.440 1.667 0.6888
## ARCH Lag[7] 1.827 2.315 1.543 0.7539
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.2475
## Individual Statistics:
## omega 0.91313
## alpha1 0.18194
## beta1 0.93748
## gamma1 0.03103
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.3062 0.1917
## Negative Sign Bias 0.2634 0.7923
## Positive Sign Bias 0.2155 0.8294
## Joint Effect 3.5807 0.3104
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 82.19 7.779e-10
## 2 30 103.55 2.620e-10
## 3 40 122.78 1.342e-10
## 4 50 113.48 4.984e-07
##
##
## Elapsed time : 0.328017Uzyskano nieistotne oszacowanie parametru \(\alpha_1\). Brak podstaw do odrzucenia hipotezy, że kwadraty wystandaryzowanych reszt są białym szumem. Wśród reszt nie występują efekty ARCH. Wartość statystyki AIC wyniosła -5.937. Model ten spełnia większość warunków, ale przez nieistotność jednego z parametrów nie może zostać uznany za poprawny.
Następnie oszacowano model EGARCH(1,2)
spec = ugarchspec(variance.model = list(model ="eGARCH",
garchOrder = c(1, 2)),
mean.model = list(armaOrder = c(0, 0), include.mean = F),
distribution.model = "norm")
egarch12 <- ugarchfit( spec = spec, data = na.omit(x1$r))
egarch12##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,2)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega -0.473371 0.154965 -3.05470 0.002253
## alpha1 0.002846 0.016941 0.16799 0.866593
## beta1 0.434133 0.041796 10.38692 0.000000
## beta2 0.509772 0.043438 11.73560 0.000000
## gamma1 0.347203 0.047007 7.38614 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega -0.473371 0.341520 -1.386071 0.165725
## alpha1 0.002846 0.036574 0.077811 0.937979
## beta1 0.434133 0.029330 14.801801 0.000000
## beta2 0.509772 0.039740 12.827822 0.000000
## gamma1 0.347203 0.107986 3.215261 0.001303
##
## LogLikelihood : 5048.045
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.9400
## Bayes -5.9240
## Shibata -5.9400
## Hannan-Quinn -5.9341
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 8.84 0.0029471
## Lag[2*(p+q)+(p+q)-1][2] 10.92 0.0010924
## Lag[4*(p+q)+(p+q)-1][5] 14.71 0.0005318
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.09459 0.7584
## Lag[2*(p+q)+(p+q)-1][8] 1.86321 0.8794
## Lag[4*(p+q)+(p+q)-1][14] 3.19578 0.9440
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.04912 0.500 2.000 0.8246
## ARCH Lag[6] 0.19067 1.461 1.711 0.9707
## ARCH Lag[8] 1.54538 2.368 1.583 0.8311
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.0623
## Individual Statistics:
## omega 0.65124
## alpha1 0.19074
## beta1 0.66619
## beta2 0.66852
## gamma1 0.03607
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.2779 0.2015
## Negative Sign Bias 0.2262 0.8211
## Positive Sign Bias 0.2442 0.8071
## Joint Effect 4.6121 0.2025
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 73.12 2.767e-08
## 2 30 93.27 1.132e-08
## 3 40 99.03 3.991e-07
## 4 50 107.24 3.106e-06
##
##
## Elapsed time : 0.4020231Uzyskane właściwości i wyniki testów w modelu są zbliżone do poprzedniego. Parametr \(\alpha_1\) również okazał się być nieistotny statystycznie. Wielkość statystyki AIC wyniosła -5.934, co oznacza, że model ten jest gorzej dopasowany niż EGARCH(1,1).
W następnej kolejności sprawdzono model EGARCH(2,1):
spec = ugarchspec(variance.model = list(model ="eGARCH",
garchOrder = c(2, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = F),
distribution.model = "norm")
egarch21 <- ugarchfit( spec = spec, data = na.omit(x1$r))
egarch21##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(2,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega -0.064126 0.001986 -32.2846 0.000000
## alpha1 -0.105961 0.032947 -3.2161 0.001299
## alpha2 0.129669 0.033452 3.8763 0.000106
## beta1 0.992359 0.000293 3385.8709 0.000000
## gamma1 0.327942 0.022781 14.3954 0.000000
## gamma2 -0.227642 0.012378 -18.3910 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega -0.064126 0.005935 -10.8041 0.000000
## alpha1 -0.105961 0.056140 -1.8874 0.059102
## alpha2 0.129669 0.057310 2.2626 0.023661
## beta1 0.992359 0.000541 1834.4038 0.000000
## gamma1 0.327942 0.044689 7.3384 0.000000
## gamma2 -0.227642 0.056867 -4.0031 0.000063
##
## LogLikelihood : 5057.781
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.9503
## Bayes -5.9311
## Shibata -5.9503
## Hannan-Quinn -5.9432
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 7.269 0.007016
## Lag[2*(p+q)+(p+q)-1][2] 8.937 0.003652
## Lag[4*(p+q)+(p+q)-1][5] 12.037 0.002724
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.05223 0.8192
## Lag[2*(p+q)+(p+q)-1][8] 6.20702 0.2151
## Lag[4*(p+q)+(p+q)-1][14] 8.11499 0.3712
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.08753 0.500 2.000 0.7673
## ARCH Lag[6] 1.86569 1.461 1.711 0.5211
## ARCH Lag[8] 2.55223 2.368 1.583 0.6301
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.8148
## Individual Statistics:
## omega 0.20006
## alpha1 0.11990
## alpha2 0.11026
## beta1 0.19528
## gamma1 0.09221
## gamma2 0.14331
##
## 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.76490 0.4444
## Negative Sign Bias 0.05949 0.9526
## Positive Sign Bias 0.36383 0.7160
## Joint Effect 0.73716 0.8644
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 74.63 1.540e-08
## 2 30 95.85 4.461e-09
## 3 40 106.76 3.236e-08
## 4 50 130.21 2.716e-09
##
##
## Elapsed time : 0.6680391Tym razem uzyskano wyłącznie istotne oszacowania parametrów modelu. Pozostałe wymagania również są spełnione przez ten model. Wartość AIC wyniosła -5.947, co jest wielkością najniższą wśród analizowanych modeli EGARCH. Zatem najlepiej dopasowanym do danych modelem będzie EGARCH(2,1). Parametry modelu:
omega=-0.064
alpha1=-0.106
alpha2=0.13
beta1=0.992
gamma1=0.328
gamma2=-0.228
Dodatnia wartość parametru \(\gamma_1\) świadczy o braku występowania asymetrii w reakcji na szoki.
Model GARCH-t
W modelach GARCH-t przyjmuje się, że wariancja ma rozkład t-Studenta zamiast rozkładu normalnego.
Oszacowano model GARCH-t(1,1):
spec <- ugarchspec(variance.model = list(model = "sGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "std")
garcht11 <- ugarchfit(spec=spec, data=na.omit(x1$r))
garcht11##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000011 0.000001 12.3532 0
## alpha1 0.128655 0.010801 11.9119 0
## beta1 0.811559 0.018723 43.3453 0
## shape 5.108278 0.549471 9.2967 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000011 0.000001 9.953 0
## alpha1 0.128655 0.012857 10.006 0
## beta1 0.811559 0.018039 44.990 0
## shape 5.108278 0.587093 8.701 0
##
## LogLikelihood : 5121.622
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.0278
## Bayes -6.0150
## Shibata -6.0278
## Hannan-Quinn -6.0231
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 9.202 0.0024170
## Lag[2*(p+q)+(p+q)-1][2] 11.164 0.0009414
## Lag[4*(p+q)+(p+q)-1][5] 14.723 0.0005274
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4849 0.4862
## Lag[2*(p+q)+(p+q)-1][5] 1.2470 0.8017
## Lag[4*(p+q)+(p+q)-1][9] 2.0740 0.8960
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.9101 0.500 2.000 0.3401
## ARCH Lag[5] 0.9868 1.440 1.667 0.7369
## ARCH Lag[7] 1.6308 2.315 1.543 0.7948
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 29.9301
## Individual Statistics:
## omega 7.6953
## alpha1 1.5363
## beta1 1.3177
## shape 0.3949
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.2903 0.1971
## Negative Sign Bias 0.4227 0.6726
## Positive Sign Bias 0.4535 0.6503
## Joint Effect 3.5519 0.3141
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 34.25 0.01717
## 2 30 49.00 0.01156
## 3 40 55.03 0.04590
## 4 50 63.37 0.08143
##
##
## Elapsed time : 0.467025Suma parametrów jest mniejsza od 1, parametry są dodatnie i istotne. Wynik testu Ljung-Boxa wskazuje na brak podstaw do odrzucenia hipotezy o tym, że kwadraty wystandaryzowanych reszt są białym szumem. P-value testu na występowanie efektów ARCH nie daje podstaw do odrzucenia hipotezy o braku występowania efektów ARCH wśród reszt. Wielkość wskaźnika AIC wyniosła -6.027, co jest najniższą uzyskaną do tej pory wartością wśród wszystkich badanych modeli. Zatem zastąpienie rozkładu normalnego rozkładem t-Studenta polepszyło dopasowanie modelu. Oszacowanie parametru shape wskazuje na liczbę stopni swobody w rozkładzie t-Studenta. Uzyskana wartość 5.108 świadczy o tym, że lepiej założyć, że reszty pochodzą z rozkładu t-Studenta, niż z rozkładu normalnego. Model ten jest oszczędny w parametry, co dodatkowo podnosi jego atrakcyjność.
Oszacowano również modele GARCH-t(1,2) i GARCH-t(2,1). Jednakże w pierwszym z nich parametr \(\beta_2\) okazał się nieistotny statystycznie, a w drugim \(\alpha_2\) było nieistotne. W związku z tym najlepszym modelem okazał się prosty model GARCH-t(1,1).
Oszacowania parametrów:
omega=0
alpha1=0.129
beta1=0.812
shape=5.108
Prognoza warunkowej wariancji
Za pomocą uzyskanych modeli można oszacować prognozę warunkowej wariancji na przyszłe okresy. Dla modelu GARCH(1,1) wyznaczono bezwarunkową wariancję korzystając ze wzoru: \[E\epsilon_t^2=\frac{\alpha_0}{1-\alpha_1-\beta_1} \].
var_uncond_1 <- garch11@fit$par[1] / (1 - garch11@fit$par[2]
- garch11@fit$par[3])
names(var_uncond_1) <- "unconditional variance"Następnie oszacowano prognozę warunkowej wariancji na kolejne 300 okresów:
fore100_1 <- predict(garch11, n.ahead = 300)Na poniższym wykresie przedstawiono oszacowania warunkowej i bezwarunkowej wariancji zwrotów dla modelu GARCH(1,1): 
Z wykresu wynika, że prognozy warunkowej wariancji zbiegają w długim okresie do poziomu wariancji bezwarunkowej.
Kolejny wykres przedstawia oszacowania warunkowej wariancji dla modelu EGARCH(2,1): 
Następny wykres odnosi się do prognoz warunkowej wariancji w modelu GARCH-t(1,1). 
Zestawienie prognoz z wszystkich trzech modeli przedstawiono poniżej: 
Prognoza warunkowej wariancji dla każdego z trzech modeli stabilizuje się na konkretnym poziomie w różnych okresach wprzód. Najwcześniej w modelu GARCH-t, bo około 50 okresu, prognoza się stabilizuje. Dla modelu GARCH jest to około 80 okresu. Z kolei w przedstawionym przedziale 300 okresów naprzód prognoza z modelu EGARCH nie stabilizuje się i nadal rośnie. W modelach GARCH i GARCH-t prognoza jest bardziej “wygładzona”, w przeciwieństwie do modelu EGARCH, gdzie jest bardziej “kanciasta”. Pronozowana wariancja warunkowa jest najniższa dla modelu GARCH-t, a najwyższa dla EGARCH, z wyjątkiem okresów od ok. 20 do ok. 190 gdzie najwyższe wartości przyjmuje prognoza z modelu GARCH.
Szacowanie wartości narażonej na ryzyko w okresie in-sample
Wyznaczenie wartości narażonej na ryzko (VaR) polega na oszacowaniu wielkości straty portfela, która wystąpi z określonym prawdopodobieństwem. Analizę należy rozpocząć od określenia przedziałów in-sample i out-of-sample. W tym przydaku zdecydowano się podzielić obserwacje w stosunku 1:13, to znaczy okres out-of-sample będzie stanowiło niecałe 6 miesięcy, od początku 2017 roku, do połowy czerwca tego roku. Aby posiadać wystarczająco dużo rekordów do szacowania modeli, za okres in-sample uznano obserwacje od połowy 2010 roku do końca 2016 roku, czyli obejmujące 6,5 roku.
Ograniczono podzbiór in-sample:
x1_ins <- x1[x1$Date <= as.Date("2017-01-02"), ]Następnie dokonano standaryzacji zwrotów z portfela:
x1_ins$rstd <- (x1_ins$r - mean(x1_ins$r, na.rm=T)) / sd(x1_ins$r ,na.rm = T)Wyznaczono 1% kwantyl empiryczny rozkładu zwrotów:
q01 <- quantile(x1_ins$rstd, 0.01, na.rm = T)Uzyskana wartość w wysokości -2.846 różniła się od 1% kwantyla rozkładu normalnego o wartości -2.326.
Następnie oszacowano modele GARCH(1,1), EGARCH(2,1) oraz GARCH-t(1,1) na danych w okresie in-sample.
ins_garch11 <- garchFit(formula = ~ garch(1, 1),
data = na.omit(x1_ins$r),
include.mean = F,
cond.dist = "norm",
trace = FALSE)
spec = ugarchspec(variance.model = list(model ="eGARCH",
garchOrder = c(2, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = F),
distribution.model = "norm")
ins_egarch21 <- ugarchfit( spec = spec, data = na.omit(x1_ins$r))
spec <- ugarchspec(variance.model = list(model = "sGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "std")
ins_garcht11 <- ugarchfit(spec=spec, data=na.omit(x1_ins$r))Wyznaczono wartość narażoną na ryzyko jako iloczyn wektora oszacowań odchylenia standardowego i 1% kwantyla empirycznego rozkładu.
x1_ins$var1 <- q01 * ins_garch11@sigma.t
x1_ins$var2 <- q01 * ins_egarch21@fit$sigma
x1_ins$var3 <- q01 * ins_garcht11@fit$sigmaNa poniższych wykresach przedstawiono zwroty wraz z oszacowaniami wartości narażonej na ryzyko dla poszczególnych modeli.
plot(x1_ins$Date, x1_ins$r, col = "red", lwd = 1, type = 'l', ylim = c(-0.1, 0.1),ylab="",
xlab="Date", main="Zwroty i VaR dla modelu GARCH(1,1)")
abline(h = 0, lty = 2)
lines(x1_ins$Date, x1_ins$var1, type = 'l', col = "green")
Dla modelu GARCH(1,1) straty przekroczyły zakładany poziom VaR w 0.567% przypadków. 
Dla modelu EGARCH(2,1) straty przekroczyły zakładany poziom VaR w 0.757% przypadków. 
Dla modelu GARCH-t(1,1) straty przekroczyły zakładany poziom VaR w 0.631% przypadków.
Na podstawie uzyskanych oszacowań można stwierdzić, że straty przekroczyły zakładany poziom wartości narażonej na ryzyko najczęściej w przypadku modelu EGARCH(2,1), a najrzadziej w modelu standardowym GARCH(1,1).
Szacowanie wartości narażonej na ryzyko w okresie out-of-sample
Prognoza wartości narażonej na ryzyko w okresie out-of-sample bazuje na jednookresowych prognozach funkcji warunkowej wariancji. Wyznaczono podzbiór:
start <- x1$obs[x1$Date == as.Date("2017-01-03")]
finish <- x1$obs[x1$Date == as.Date("2017-06-16")]
x1_a <-x1[start:finish, ]Dla modelu GARCH(1,1) należy zaimplementować funkcję:
VaR <- rep(NA, times = finish - start + 1)
for (k in start:finish) {
tmp.data <- x1[x1$obs <= (k - 1), ]
tmp.data <- tmp.data[as.Date("2010-08-01") <= tmp.data$Date, ]
tmp.data$rstd <- (tmp.data$r - mean(tmp.data$r, na.rm = T)) /
sd(tmp.data$r, na.rm = T)
q01 <- quantile(tmp.data$rstd, 0.01, na.rm = T)
t.garch11 <- garchFit(formula = ~ garch(1, 1),
data = na.omit(tmp.data$r),
include.mean = F,
cond.dist = "norm",
trace = FALSE)
sigma.forecast <- predict(t.garch11, n.ahead = 1)
sigma.forecast2 <- sigma.forecast$standardDeviation
VaR[k - start + 1] <- q01 * sigma.forecast2
}
x1_a$var1 <- VaRPoniżej przedstawiono wykres prognozy:
plot(x1_a$Date, x1_a$r, col = "red", lwd = 1, type = 'l',
ylim = c(-0.05, 0.05),xlab="Time",ylab="",
main="Zwroty z portfela i VaR w okresie out-of-sample
w modelu GARCH(1,1)")
abline(h = 0, lty = 2)
lines(x1_a$Date, x1_a$var1, type = 'l', col = "green")
Dla modelu GARCH(1,1) straty przekroczyły zakładany poziom VaR w 0% przypadków.
Następnie oszacowano prognozę z zastosowaniem modelu EGARCH(2,1):
Poniżej przedstawiono wykres prognozy: 
Dla modelu EGARCH(2,1) straty przekroczyły zakładany poziom VaR w 0.893% przypadków.
Przeprowadzono prognozę dla modelu GARCH-t(1,1):
Poniżej przedstawiono wykres prognozy: 
Dla modelu GARCH-t(1,1) straty przekroczyły zakładany poziom VaR w 0% przypadków.
Zestawienie prognoz wartości narażonej na ryzyko dla wszystkich trzech modeli kształtuje się następująco:
plot(x1_a$Date, x1_a$var1, col = "green", lwd = 1, type = 'l',
ylim = c(-0.05, -0.02),xlab="Date",ylab="",
main="Value at Risk comparison between models
in out-of-sample period ")
lines(x1_a$Date, x1_a$var2, type = 'l', col = "brown")
lines(x1_a$Date, x1_a$var3, type = 'l', col = "blue")
legend("bottomleft", legend=c("EGARCH(2,1)","GARCH-t(1,1)",
"GARCH(1,1)"), lwd=c(2,2), col=c("brown","blue","green"))
Jedynie dla modelu EGARCH(2,1) straty przekroczyły zakładany poziom VaR, ale zaledwie 1 raz w badanym okresie. Dla pozostałych modeli straty nigdy nie przeproczyły tego poziomu. Oznacza to, że wyznaczone wartości narażone na ryzyko zawierają w sobie zbyt duży margines straty, gdyż takie wielkości strat z portfela dla modeli GARCH i GARCH-t w badanym okresie wystąpiły z zerowym prawdopodobieństwem.
Analiza wrażliwości
W tej częsci badania postanowiono przeanalizwać prognozy wartości narażonej na ryzyko w zależności od długości przedziału out-of-sample. Do analizy wybrano jedynie model GARCH-t(1,1), ponieważ pod względem kryteriów informacyjnych był to model najlepiej dopasowany do danych. Również funkcja wyznaczająca prognozy VaR najszybciej sobie radzi w przypadku tego modelu. Wybrano okresy 3,12 i 24 miesięczne out-of-sample oraz 6-miesięczny, uwzględniony w poprzedniej części badania.
Wyznaczono prognozę 12-miesięczną:
start <- x1$obs[x1$Date == as.Date("2016-07-01")]
finish <- x1$obs[x1$Date == as.Date("2017-06-16")]
x1_w1 <-x1[start:finish, ]
VaR <- rep(NA, times = finish - start + 1)
for (k in start:finish) {
tmp.data <- x1[x1$obs <= (k - 1), ]
tmp.data <- tmp.data[as.Date("2010-08-01") <= tmp.data$Date, ]
tmp.data$rstd <- (tmp.data$r - mean(tmp.data$r, na.rm = T)) /
sd(tmp.data$r, na.rm = T)
q01 <- quantile(tmp.data$rstd, 0.01, na.rm = T)
spec <- ugarchspec(variance.model = list(model = "sGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "std")
t.garcht11 <- ugarchfit(spec=spec, data=na.omit(tmp.data$r))
sigma.forecast <- ugarchforecast(t.garcht11, n.ahead = 1)
sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
VaR[k - start + 1] <- q01 * sigma.forecast2
}
x1_w1$var <- VaR
Dla modelu GARCH-t(1,1) w przypadku 12-miesięcznego okresu out-of-sample straty przekroczyły zakładany poziom VaR w 0% przypadków.
Dla prognozy 3-miesięcznej:
Dla modelu GARCH-t(1,1) w przypadku 3-miesięcznego okresu out-of-sample straty przekroczyły zakładany poziom VaR w 0% przypadków.
Wyznaczono jeszcze prognozę 24-miesięczną: 
Dla modelu GARCH-t(1,1) w przypadku 24-miesięcznego okresu out-of-sample straty przekroczyły zakładany poziom VaR w 0.402% przypadków.
Dla analizowanego aktywa, jedynie podczas analizy 24-miesięcznego okresu out-of-sample wielkość straty w ogóle przekroczyła zakładany poziom wartości narażonej na ryzyko. W tym przypadku stosunek okresu prognozowanego (out-of-sample) do okresu wykorzystanego do estymacji modelu (in-sample) wyniósł 2:5, co jest wartością najwyższą w porównaniu do pozostałych wariantów. Można przypuszczać, że w miarę wydłużania się analizowanego okresu prognozy rośnie prawdopodobieństwo pojawienia się nieprzewidzianych zjawisk i szoków.
Podsumowanie
Przeprowadzone badanie umożliwiło sformułowanie kilku wniosków. Po pierwsze, wśród modeli klasy GARCH i ich rozszerzeń, najprostsze modele, czyli takie z najmniejszą liczbą parametrów, okazały się być najlepiej dopasowane do analizowanych danych. Wariancja zwrotów ze zbudowanego portfela była najlepiej opisywana poprzez model GARCH-t(1,1), uwzględniający t-Studenta rozkład reszt, zamiast rozkładu normalnego. Różnice pomiędzy poszczególnymi modelami zidentyfikowano na przykładzie oszacowań warunkowej wariancji. Oszacowanie wartości narażonej na ryzyko pozwoliło na określenie maksymalnego poziomu straty inwestora z analizowanego portfela aktywów. Brak założeń wynikających z teorii ekonomii w przypadku modelowania warunkowej wariancji modelami klasy GARCH umożliwia analizę portfela zbudowanego z zupełnie różnych aktywów, zazwyczaj niepowiązanych wspólnymi mechanizmami rynkowymi. Dzięki temu można modelować zmienność dowolnych portfeli inwestycyjnych.