Analiza porównawcza ryzyka dla portfela akcji z NYSE i NASDAQ w modelach klasy GARCH
Vladyslav Sukhovetskyi
Czerwiec, 2017
Streszczenie
Niniejsza praca jest raportem z przeprowadzonej analizy porównawczej ryzyka rozumianego jako oszacowanie funkcji warunkowej wariancji w modelach klasy GARCH. Analiza porównawcza została przeprowadzona za pomocą zestawienia ze sobą standardowego modelu GARCH oraz jego rozszerzeń - EGARCH i TARCH. Do danych wejściowych należą akcji spółek NASDAQ i NYSE. Dla portfela zrównoważonego przeprowadzono analizę porównawczą oszacowań warunkowej wariancji, oraz oszacowań wartości narażonej na ryzyko w okresie in-sample i out-sample. Także została przeprowadzona krótka analiza wrażliwości estymacji parametrów modelu w zależności od długości próby.
Przygotowanie danych
Najpierw zostały wczytane niezbędne pakiety oraz dodatkowe funkcje, które zostały wykorzystane w analizie:
library(fBasics)
library(tseries)
library(car)
library(FinTS)
library(fGarch)
library(rugarch)
library(xts)
library(plyr)
library(dygraphs)
source("functions.R")Do analizy zostało wykorzystano dane dotyczące akcji amerykańskich spółek NASDAQ i NYSE. Notowania zostały pobrane z serwisu “Yahoo Finance”. Portfel zrównoważony składa się ze 6 spółek giełdowych, które pochodzą z rożnych branż. Są to:
- Spółki z branży technologicznej:
- Apple (AAPL)
- Microsoft (MSFT)
- Spółki z branży opieki zdrowotnej:
- Johnson & Johnson (JNJ)
- HUMANA INC (HUM)
- Spółki z branży finansowej:
- Bank of America Corporation (BAC)
- JPMorgan Chase & Co. (JPM)
Następnie zaimportowano surowe dane dotyczące notowań poszczególnych spółek za pomocą funkcji import. Wybrano jedynie ceny zamknięcia i datę:
import <- function(plik){
url <- plik
plik <- read.csv(url,
header = TRUE,
sep = ",",
dec = ".",
stringsAsFactors = F)
plik$Date <- as.Date(as.character(plik$Date),
format = "%m/%d/%Y")
plik <- plik[, c("Date", "Close")]
plik <- plik[as.Date("1986-03-13") <= plik$Date, ]
plik <- plik[plik$Date <= as.Date("2017-06-20"), ]
plik$r <- diff.xts(log(plik$Close))
plik <- plik[, c("Date","Close", "r")]
return(plik)
}
AAPL <- import("data/AAPL.csv")
BAC <- import("data/BAC.csv")
HUM <- import("data/HUM.csv")
JNJ <- import("data/JNJ.csv")
JPM <- import("data/JPM.csv")
MSFT <- import("data/MSFT.csv")Powyższa funkcja także uwzględnia obcięcie okresu szeregu z dołu i z góry. Wyznaczono początkową datę 1986-03-13 jako wartość najwyższą z początkowych dat każdego z szeregów czasowych. Także, zostały policzone logarytmiczne stopy zwrotu dla poszczególnych spółek.
Analiza graficzna wejściowych danych
W tej części przedstawiono wykresy dziennych notowań dla poszczególnych akcji. Wstępna analiza wykresów świadczy o nie stacjonarności szeregów. Można zobaczyć zmienność średniej i wariancji we wszystkich szeregach.
par(mfrow = c(2, 1))
plot(AAPL$Date, AAPL$Close,
type = "l", col="red", lwd = 1,
main = "Dzienne notowania Apple"
, xlab = "", ylab = "")
plot(MSFT$Date, MSFT$Close,
type = "l", col="red", lwd = 1,
main = "Dzienne notowania Microsoft"
, xlab = "", ylab = "")
plot(HUM$Date, HUM$Close,
type = "l", col="red", lwd = 1,
main = "Dzienne notowania HUMANA INC"
, xlab = "", ylab = "")
plot(JNJ$Date, JNJ$Close,
type = "l", col="red", lwd = 1,
main = "Dzienne notowania Johnson & Johnson"
, xlab = "", ylab = "")
plot(BAC$Date, BAC$Close,
type = "l", col="red", lwd = 1,
main = "Dzienne notowania Bank of America Corporation"
, xlab = "", ylab = "")
plot(JPM$Date, JPM$Close,
type = "l", col="red", lwd = 1,
main = "Dzienne notowania JPMorgan Chase & Co."
, xlab = "", ylab = "")
Patrząc na wykresy analizowanych zmiennych, widać, że nie są one stacjonarne (zarówno średnia jak i wariancja nie wydają się być stałe w czasie). Z powyższych grafów jest widoczne podobne poruszanieszie się dla spółek Apple, HUMANA i Johnson & Johnson. Nieco inny kształt przyjmują szeregi branży finansowej oraz Microsoft. Dla wszystkich spółek widoczny charakterystyczny spadek wartości akcji przy 2010 i 2000 roku, które były spowodowane kryzysem.
Budowa portfela inwestycyjnego
W poniższej części raportu przedstawiony proces formułowania portfolio inwestycyjnego obejmującego wyższe zaznaczone spółki. Najpierw zmieniono nazwy cen zamknięcia dla poszczególnych giełd a potem wrzucono ich do funkcji join_all, która włączy spotrzebowane tablice.
portfel <- join_all(list(AAPL,BAC,HUM,JNJ,JPM,MSFT), by = 'Date', type = 'full')
portfel$Mean <- rowMeans(portfel[,3:8])Zostały zmieniony nazwy kolumn i dodano numer obserwacji:
portfolio <- data.frame(portfel$Date,portfel$Mean, 1:length(portfel$Mean))
colnames(portfolio) <- c("Date","r","obs")Przedstawiono fragment struktury danych:
## 'data.frame': 7884 obs. of 3 variables:
## $ Date: Date, format: "1986-03-13" "1986-03-14" ...
## $ r : num NA 0.024 -0.00651 -0.00628 -0.00269 ...
## $ obs : int 1 2 3 4 5 6 7 8 9 10 ...Analiza graficzna portfela
Przedstawiono wykres dziennych zwrotów z portfela. Jak widać, największa wahanie wśród zwrotów występuje w 1996-2003 i 2008-2010. Takie wahanie spowodowane są gwałtownymi zmianami wartości akcji poszczególnych spółek w odpowiadającym okresie czasowym.
library(xts)
portfolio.xts <- xts(portfolio[, c("r", "Date")], order.by = portfolio$Date)
dygraph(portfolio.xts$r) %>% dyAxis("y", label = "Stopa zwrotu", valueRange = c(-0.13, 0.15)) %>% dyRangeSelector(height = 40 )Diagnostyka portfela
Przeprowadzono testy statystyczne i diagnostyczne w celu wykrycia w analizowanym szeregu efektów (cech charakterystycznych) modelu GARCH.
Najpierw statystyki opisowe:
## X..portfolio.r
## nobs 7884.000000
## NAs 1.000000
## Minimum -0.237457
## Maximum 0.116876
## 1. Quartile -0.006785
## 3. Quartile 0.008349
## Mean 0.000557
## Median 0.000974
## Sum 4.393439
## SE Mean 0.000171
## LCL Mean 0.000222
## UCL Mean 0.000893
## Variance 0.000231
## Stdev 0.015191
## Skewness -0.770271
## Kurtosis 14.877237Jak widać, średnia rozkładu jest zbliżoną do zera. Kurtoza wynosi 14.87, co różni ją od wartości 3 dla rozkładu normalnego.
Poniżej przedstawiono histogram z nałożoną gęstością rozkładu normalnego dla badanego portfolio.
hist(portfolio$r, prob = T, breaks = 100, xlab = "Zwoty", main = "Rozkład zwrotów portfela")
curve(dnorm(x,
mean = mean(portfolio$r, na.rm = T),
sd = sd(portfolio$r, na.rm = T)),
col = "darkblue", lwd = 2, add = TRUE)
Z analizy graficznej wykresu można stwierdzić, że rozkład portfelu ma zbyt duże skupienie wokół średniej i grube ogony. Także są widoczne autlajery, niektóre wartości zwrotów przyjmują wartość do -0,25. Więc, można stwierdzić o wystąpieniu leptokurtyczności dla analizowanego portfelu.
Następnie przedstawiono wykresy funkcji autokorelacji dla zwrotów z portfela.
acf(portfolio$r, lag.max = 36, na.action = na.pass,
xlab = "Opóżnienia",
ylim = c(-0.1, 0.1),
col = "darkblue", lwd = 7,
main = "Wykres ACF zwrotów portfela")
Jak widać z wykresu, występują istotne opóźnienia. Są to opóźnienia: 2,7,15,16,18,21,31,34. Więc, wartości zwrotów zależą istotnie od jej przeszłych wartości. Występuje autokorelacja wśród reszt - reszty nie są białym szumem.
Wykres funkcji autokorelacji dla kwadratów zwrotów.
acf(portfolio$r^2, lag.max = 100, na.action = na.pass,
xlab = "Opóżnienia",
ylim = c(0, 0.5),
col = "darkblue", lwd = 7,
main = "Wykres ACF kwadratów zwrotów portfela")
Analogicznie do poprzedniego wykresu, wśrów kwadratów logarytmicznych stop zwrotów także występuje silna autokorelacja, co wskazuje na występowanie zależności wyższego rzędu. Więc kwadraty reszt nie są białym szumem, co świadczy o wystąpieniu efektów GARCH dla analizowanego portfela.
Wykres quantile:
qqnorm(portfolio$r)
qqline(portfolio$r, col = 2)
Rozkład cechuję się leptokurtucznością, ponieważ wykres kwantuliu przyjmuje formę litery “s”.
Inne statystyki testowe
- Statystyka Jarque-Bery
jarque.bera.test(na.omit(portfolio$r))##
## Jarque Bera Test
##
## data: na.omit(portfolio$r)
## X-squared = 73523, df = 2, p-value < 2.2e-16Przy poziomie istotności 5%, odrzucamy \(H_0\) o normalności zwrotów [2.2e-16 < 0,05].
- Statystyka Durbina-Watsona
durbinWatsonTest(lm(formula = portfolio$r ~ 1),
max.lag = 5)## lag Autocorrelation D-W Statistic p-value
## 1 -0.017312447 2.034268 0.146
## 2 -0.022036317 2.043587 0.058
## 3 -0.015745225 2.030979 0.170
## 4 0.010159627 1.979154 0.366
## 5 -0.008236644 2.015804 0.468
## Alternative hypothesis: rho[lag] != 0Według testu sprawdźmy autokorelację dla 5 pierwszych opóźnień. Przy poziomie istotności 5%, nie da się odrzucić \(H_0\) że pierwsze 5 współczynników autokorelacji wśród reszt równe zeru [p-value > 0,05]. Podtwierdzono fakt autokorelacji wśród reszt.
- Test na występowanie efektów ARCH wśród zwrotów
ArchTest(portfolio$r, lags = 5)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: portfolio$r
## Chi-squared = 729.35, df = 5, p-value < 2.2e-16Wynik testu na występowanie efektów ARCH wśród zwrotów sugeruje odrzucenie hipotezy (na poziomie 5%) o braku występowania efektów ARCH.
- Statystykę DW dla kwadratów zwrotów
durbinWatsonTest(lm(formula = portfolio$r ^ 2 ~ 1),
max.lag = 5)## lag Autocorrelation D-W Statistic p-value
## 1 0.1630374 1.673905 0
## 2 0.1933192 1.613337 0
## 3 0.1510067 1.697949 0
## 4 0.1140503 1.771848 0
## 5 0.2148838 1.570173 0
## Alternative hypothesis: rho[lag] != 0Sprawdźmy autokorelację dla 5 pierwszych opóźnień kwadratów zwrotów. Według statystyki wszystkie opóźnienia silnie istotne.
Wybór właściwej specyfikacji modelu
W tym podrozdziale, dla analizowanego portfela została przeprowadzona analiza porównawcza oszacowań warunkowej wariancji. Zostały porównywane dwa rozszerzenia GARZH w odniesieniu do modelu bazowego jakim jest model GARCH(p, q).
Zdecydowano na wybór długości badanego okresu z 2011-12-20 po 2016-12-20, ponieważ ten okres cechuje się bardziej aktualnymi danymi, które są ciekawsze dla właściwej analizy.
portfolio2 <- portfolio[portfolio$Date <= as.Date("2016-12-20"), ]
portfolio2 <- portfolio2[as.Date("2011-12-20") <= portfolio2$Date, ]Wykres nowej próby:
Dalej przeprowadzona standaryzacja zwrotów
portfolio2$rstd <- (portfolio2$r - mean(portfolio2$r, na.rm=T)) /
sd(portfolio2$r ,na.rm = T)
tail(portfolio2)## Date r obs rstd
## 7755 2016-12-13 7.079220e-03 7755 0.62141289
## 7756 2016-12-14 -1.167523e-03 7756 -0.19913119
## 7757 2016-12-15 7.506113e-03 7757 0.66388834
## 7758 2016-12-16 -5.446321e-03 7758 -0.62486794
## 7759 2016-12-19 -8.347357e-06 7759 -0.08379414
## 7760 2016-12-20 3.410425e-03 7760 0.25637080Poniżej przedstawiono histogram z nałożoną gęstością rozkładu normalnego dla standaryzowanych zwrotów.
hist(portfolio2$rstd, prob = T, main = 'Rozkład standartyzowanych zwrotów', breaks = 40, xlab = "Standartyzowane zwroty", ylab = "Density")
curve(dnorm(x, mean = mean(portfolio2$rstd, na.rm = T),
sd = sd(portfolio2$rstd, na.rm = T)),
col = "darkblue", lwd = 2, add = TRUE)
Z analizy graficznej wykresu można stwierdzić, że rozkład portfelu ma zbyt duże skupienie wokół średniej i grube ogony. Także są widoczne autlajery, niektóre wartości zwrotów przyjmują wartość do -4 i +4. Więc, można stwierdzić o wystąpieniu leptokurtyczności dla analizowanego portfelu.
Wartość 1% kwantyliu empirycznego oraz standardowego rozkładu normalnego
q01 <- quantile(portfolio2$rstd, 0.01, na.rm = T)
q01## 1%
## -2.602548qnorm(0.01, 0, 1)## [1] -2.326348Pozostałe testy statystyczne i diagnostyczne przeprowadzone w celu wykrycia w analizowanej skróconej próbie efektów modelu GARCH. Wszystkie testy wskazują na wystąpienia autokorelacji wśród stop zwrotów i o występowaniu efektów GARCH.
Model zwykły GARCH
Jest to uogólnienie modelu ARCH. Model GARCH dopuszcza, aby wartości warunkowej wariancji zależały także od swoich własnych opóźnień.
Podstawowy model GARCH(q,p): \[ \sigma^2_t = \alpha_0 +
\sum_{i=1}^{q} \alpha_i \epsilon^2_{t-i} + \sum_{j=1}^{p} \beta_j \sigma^2_{t-j}
\] Wymagania prawidłowej specyfikacji modelu GARCH(p,q):
- warunek stacjonarności procesów GARCH:
\[\sum_{i=1}^{q} \alpha_i + \sum_{j=1}^{p} \beta_j < 1\]
- aby wariancja była dodatnia, oszacowania parametrów powinny być nieujemne
- kwadraty wystandaryzowanych reszt nie mogą podlegać autokorelacji
- parametry modelu powinny być istotne
- dopasowanie modeli porównuje się za pomocą kryteriów informacyjnych AIC i BIC
Wady modelu GARCH(q,p):
- Nie uwzględnia asymetrycznego wpływ dobrych i złych wiadomości.
- Nie daje możliwości modelowania efektu dźwigni, gdyż warunkowa wariancja zależy jedynie od wartości bezwzględnych wcześniejszych realizacji (nie uwzględnia znaków).
- Nie pozwala na uwzględnienie efektu długiej pamięci.
Rozpoczęto szacowanie modelu GARCH(1,1) z jednym opóźnieniem szoku (q = 1) i jednym opóźnieniem warunkowej wariancji (p = 1).
spec <- ugarchspec(variance.model = list(model = "sGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "norm")
garch11 <- ugarchfit(spec = spec,
data = na.omit(portfolio2$r))
garch11##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega 0.00001 0.000000 26.742 0
## alpha1 0.16410 0.016166 10.150 0
## beta1 0.73707 0.021660 34.030 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.00001 0.000000 23.2448 0
## alpha1 0.16410 0.016527 9.9291 0
## beta1 0.73707 0.023936 30.7940 0
##
## LogLikelihood : 4064.775
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.4524
## Bayes -6.4401
## Shibata -6.4524
## Hannan-Quinn -6.4478
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4871 0.4852
## Lag[2*(p+q)+(p+q)-1][2] 0.4963 0.6953
## Lag[4*(p+q)+(p+q)-1][5] 1.2855 0.7924
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1251 0.7235
## Lag[2*(p+q)+(p+q)-1][5] 0.3055 0.9831
## Lag[4*(p+q)+(p+q)-1][9] 3.0163 0.7563
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1706 0.500 2.000 0.6796
## ARCH Lag[5] 0.3079 1.440 1.667 0.9381
## ARCH Lag[7] 3.3503 2.315 1.543 0.4500
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 23.6179
## Individual Statistics:
## omega 7.8632
## alpha1 0.1303
## beta1 0.2274
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 0.846 1.01 1.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.3468 0.17829
## Negative Sign Bias 0.1668 0.86753
## Positive Sign Bias 0.9352 0.34986
## Joint Effect 7.6291 0.05433 *
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 39.47 0.003833
## 2 30 53.07 0.004137
## 3 40 57.82 0.026553
## 4 50 80.83 0.002827
##
##
## Elapsed time : 0.6000812Uzyskano oceny parametrów: \(\omega=\) 1.039293e-05 \(\alpha_1=\) 0.1640958 \(\beta_1=\) 0.7370698
Suma parametrów modelu jest mniejsza od 1, oraz wszystkie parametry są dodatnie i istotne. P-value testu Ljung-Boxa dla kwadratów wystandaryzowanych reszt wskazuje na brak podstaw do odrzucenia hipotezy(przy poziomie istotności 5%), ż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 ten spełnia wszystkie warunki i wydaje się być dobry.
Oszacowanie modeli GARCH(1,2) i GARCH(2,1) doprowadziło do uzyskania nieistotnych parametrów przy drugich opóźnieniach(\(\alpha_2 i \beta_2\)). Także uzyskano większe wartości kryterium informacyjnych. W związku z tym najlepiej dopasowanym do analizowanych danych modelem z klasy GARCH jest model GARCH(1,1).
Model GARCH-t
Właściwy model zakłada, że \(\mu\) ma rozkład t-Studenta zamiast rozkładu normalnego.
- Rozkład t-Studenta z małą liczbą stopni swobody (tutaj 2) ma grubsze ogony.
- Im więcej stopni swobody, tym bardziej rozkład t-Studenta jest zbliżony do rozkładu normalnego.
Wymagania prawidłowej specyfikacji modelu GARCH(p,q) są takie same jak i w zwykłego modelu GARCH(p,q)
Rozpoczęto szacowanie modelu GARCH-t(1,1) z jednym opóźnieniem szoku (q = 1) i jednym opóźnieniem warunkowej wariancji (p = 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")
k.garcht11 <- ugarchfit(spec = spec,
data = na.omit(portfolio2$r))
k.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.000009 0.000001 13.2059 0
## alpha1 0.145650 0.017406 8.3680 0
## beta1 0.768138 0.023935 32.0926 0
## shape 7.990292 1.577892 5.0639 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000009 0.000001 12.4302 0e+00
## alpha1 0.145650 0.014972 9.7284 0e+00
## beta1 0.768138 0.022959 33.4572 0e+00
## shape 7.990292 1.628041 4.9079 1e-06
##
## LogLikelihood : 4080.613
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.4760
## Bayes -6.4596
## Shibata -6.4760
## Hannan-Quinn -6.4698
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4623 0.4965
## Lag[2*(p+q)+(p+q)-1][2] 0.4673 0.7091
## Lag[4*(p+q)+(p+q)-1][5] 1.2568 0.7993
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.459e-05 0.9947
## Lag[2*(p+q)+(p+q)-1][5] 2.425e-01 0.9892
## Lag[4*(p+q)+(p+q)-1][9] 2.986e+00 0.7613
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.2668 0.500 2.000 0.6055
## ARCH Lag[5] 0.4435 1.440 1.667 0.9002
## ARCH Lag[7] 3.5671 2.315 1.543 0.4133
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 29.2177
## Individual Statistics:
## omega 6.6634
## alpha1 0.1390
## beta1 0.2464
## shape 0.1463
##
## 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.2511 0.21112
## Negative Sign Bias 0.1415 0.88752
## Positive Sign Bias 0.7301 0.46548
## Joint Effect 6.9680 0.07293 *
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 26.54 0.11573
## 2 30 40.01 0.08382
## 3 40 40.60 0.39960
## 4 50 65.19 0.06069
##
##
## Elapsed time : 0.71684Uzyskano oceny parametrów: \(shape=\) 7.990292 \(\omega=\) 9.062614e-06 \(\alpha_1=\) 0.1456503 \(\beta_1=\) 0.7681382
Suma parametrów modelu jest mniejsza od 1, wszystkie parametry są dodatnie i istotne. P-value testu Ljung-Boxa dla kwadratów wystandaryzowanych reszt wskazuje na brak podstaw do odrzucenia hipotezy(przy poziomie istotności 5%), ż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 ten spełnia wszystkie warunki i wydaje się być dobry.
Parametr \(shape\) określa liczbę stopni swobody w rozkładzie t-Studenta. W naszym przypadku jest to wartość 7.99.Więc, lepiej założyć, że reszty pochodzą z rozkładu t-Studenta aniżeli z rozkładu normalnego.
Oszacowanie modeli GARCH-t(1,2) i GARCH-t(2,1) doprowadziło do uzyskania nieistotnych parametrów przy drugich opóźnieniach(\(\alpha_2 i \beta_2\)). Także uzyskano większe wartości kryterium informacyjnych. W związku z tym najlepiej dopasowanym do analizowanych danych modelem z klasy GARCH-t jest model GARCH-t(1,1).
Model TGARCH
Funkcja warunkowej wariancji w tym modelu ma postać: \[ \sigma^2_t = \alpha_0 + \sum_{i=1}^{q} \alpha_i \epsilon^2_{t-i} + \sum_{j=1}^{p} \beta_j \sigma^2_{t-j} + \sum_{i=1}^{r} \lambda_k D_{t-k} \epsilon^2_{t-i} \] gdzie: \(D_{t-k}\) to zmienna zerojedynkowa przyjmująca wartość 1 gdy reszta w okresie (t − k) była ujemna i 0 w pozostałych przypadkach. Istotność parametru \(\lambda_k\) świadczy o asymetrycznej reakcji warunkowej wariancji na wiadomości napływające na rynek.
Najpierw szacujemy model TGARCH(2,1) z jednym opóźnieniem szoku (q = 2) i jednym opóźnieniem warunkowej wariancji (p = 1).
spec <- ugarchspec(
variance.model = list(model ="fGARCH",
garchOrder = c(2,1),
submodel = "TGARCH"),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "norm")
tgarch11 <- ugarchfit( spec = spec,
data = na.omit(portfolio2$r))
tgarch11##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : fGARCH(2,1)
## fGARCH Sub-Model : TGARCH
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega 0.001210 0.000174 6.9572 0.000000
## alpha1 0.080151 0.014280 5.6126 0.000000
## alpha2 0.106404 0.037710 2.8217 0.004777
## beta1 0.736225 0.022822 32.2594 0.000000
## eta11 1.000000 0.319717 3.1278 0.001761
## eta12 -0.335714 0.141620 -2.3705 0.017762
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.001210 0.000321 3.76398 0.000167
## alpha1 0.080151 0.037751 2.12313 0.033743
## alpha2 0.106404 0.073045 1.45670 0.145200
## beta1 0.736225 0.015219 48.37565 0.000000
## eta11 1.000000 0.560578 1.78387 0.074444
## eta12 -0.335714 0.359800 -0.93306 0.350791
##
## LogLikelihood : 4091.649
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.4903
## Bayes -6.4658
## Shibata -6.4904
## Hannan-Quinn -6.4811
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.3244 0.5690
## Lag[2*(p+q)+(p+q)-1][2] 0.6243 0.6380
## Lag[4*(p+q)+(p+q)-1][5] 1.7550 0.6773
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4496 0.5025
## Lag[2*(p+q)+(p+q)-1][8] 2.3515 0.8016
## Lag[4*(p+q)+(p+q)-1][14] 6.0962 0.6259
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.2256 0.500 2.000 0.6348
## ARCH Lag[6] 0.5466 1.461 1.711 0.8789
## ARCH Lag[8] 4.8707 2.368 1.583 0.2629
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.6319
## Individual Statistics:
## omega 0.2135
## alpha1 0.1898
## alpha2 0.1535
## beta1 0.1629
## eta11 0.5639
## eta12 0.1065
##
## 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.56270 0.5737
## Negative Sign Bias 0.32779 0.7431
## Positive Sign Bias 0.09149 0.9271
## Joint Effect 0.37546 0.9453
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 43.64 0.001060
## 2 30 56.12 0.001832
## 3 40 59.67 0.018159
## 4 50 80.59 0.002981
##
##
## Elapsed time : 1.06561Uzyskano oceny parametrów: \(\eta_1=\) -0.3357142 \(\eta_2=\) -0.3357142 \(\omega=\) 0.00120992 \(\alpha_1=\) 0.0801506 \(\beta_1=\) 0.7362247
Suma parametrów modelu jest mniejsza od 1, wszystkie parametry są dodatnie i istotne. P-value testu Ljung-Boxa dla kwadratów wystandaryzowanych reszt wskazuje na brak podstaw do odrzucenia hipotezy(przy poziomie istotności 5%), ż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 ten spełnia wszystkie warunki i wydaje się być dobry.
Parametr \(\eta_2\) oświadczy o asymetrycznej reakcji warunkowej wariancji na wiadomości napływające na rynek. Ten spółczynnik jest większy od zera i istotny, więc model uwzględnia asymetryczną reakcję warunkowej wariancji.
Oszacowanie modeli TGARCH(2,1) doprowadziło do uzyskania nieistotnych parametrów przy drugich opóźnieniach(\(\alpha_2\)). Także uzyskano większe wartości kryterium informacyjnych. W modelu TGARCH(1,1) uzyskano istotność wszystkich parametrów w modelu, ale wartości kryterium informacyjnych włąściwego modelu były niższe. W związku z tym najlepiej dopasowanym do analizowanych danych modelem z klasy TGARCH jest model TGARCH(2,1).
Prognoza warunkowej wariancji
Za pomocą uzyskanych modeli można oszacować prognozę warunkowej wariancji na przyszłe okresy. Oszacowano prognozę warunkowej wariancji na kolejne 300 okresów:
var_uncond <- garch11@fit$coef[1] / (1 - garch11@fit$coef[2]
- garch11@fit$coef[3])
names(var_uncond) <- "unconditional variance"
k.fore100 <- ugarchforecast(garch11, n.ahead = 100)
# wykres oszacowań i prognoz warunkowej wariancji w długim okresie
plot(k.fore100@forecast$sigmaFor ^ 2, type = "l", xlab = "Okres prognozy", ylab = "Wartość warunkowej wariancji")
abline(h = var_uncond, col = "red", lty = 2)
title(main = "Warunkowa i bezwarunkowa wariancja zwrotów portfela")
Oszacowania wartości narażonej na ryzyko (VaR) dla wybranych modeli w okresie in-sample
Wartość narażone na ryzyko (VaR) polega na oszacowaniu wielkości straty portfela, która wystąpi z określonym prawdopodobieństwem.
Przed rozpoczęciem oszacowania wartości narażonej na ryzyko (VaR) dla wybranych modeli najpierw określono przedział in-sample i out-of-sample. Zdecydowano się uwzględnić końcową część próbki. Więc, okres in-sample będzie stanowił okres czasu - 1.5 roku, od maja 2015 roku, do grudnia 2016 roku. Za okres out-of-sample zdecydowano się uwzględnić obserwacje od grudnia 2016 roku do 20 czerwca 2017 roku, czyli 6 miesiące. Taki wybór okresów spowodowany chęcią zbadać zachowanie się VaR na aktualnych danych, i jakich zmian wartości narażonej na ryzyko należy się spodziewać w najbliższej przyszłości.
Najpierw skrócono próbę do okresu in-sample:
portfolio3 <- portfolio[portfolio$Date <= as.Date("2016-12-20"), ]
portfolio3 <- portfolio3[as.Date("2015-05-20") <= portfolio3$Date, ]Dodano standaryzowane reszty i policzono kwantyl empiryczny dla portfelu w okresie in-sample:
portfolio3$rstd <- (portfolio3$r - mean(portfolio3$r, na.rm = T)) /
sd(portfolio3$r, na.rm = T)
q012 <- quantile(portfolio3$rstd, 0.01, na.rm = T)
q012## 1%
## -2.661102Jak widać, empiryczny kwantyl (-2.66) różni się od kwantuliu standardowego rozkładu normalnego o wartości -2.33.
Następnie oszacowano modele GARCH(1,1), GARCH-t(1,1) oraz TGARCH(2,1) na danych w okresie in-sample.
- GARCH(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 = "norm")
k2.garch11 <- ugarchfit(spec = spec, data = na.omit(portfolio3$r))Dla modelu GARCH(1,1)wyznaczono VaR jako iloczyn wektora oszacowań odchylenia standardowego i 1% kwantyla empirycznego rozkładu. Także zostały przedstawione zwroty wraz z oszacowaniami wartości narażonej na ryzyko dla modelu.
portfolio3$VaR <- q012 * k2.garch11@fit$sigma
plot(portfolio3$Date, portfolio3$r, col = "red",
lwd = 1, type = 'l', ylim = c(-0.07, 0.045), main = 'Zwroty i oszacowania VaR modelu GARCH(1,1)', xlab = "Rok", ylab = "Wartość zwrotów")
abline(h = 0, lty = 2)
lines(portfolio3$Date, portfolio3$VaR, type = 'l', col = "green")
Dla modelu GARCH(1,1) straty przekroczyły zakładany poziom VaR w 0.007462% przypadków - 3 razy.
- 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")
k.garcht11 <- ugarchfit(spec = spec,
data = na.omit(portfolio3$r))Dla modelu GARCH-t(1,1) wyznaczono VaR. Także zostały przedstawione zwroty wraz z oszacowaniami wartości narażonej na ryzyko dla modelu.
Dla modelu GARCH-t(1,1) straty przekroczyły zakładany poziom VaR w 0.007461% przypadków - 3 razy.
- TGARCH(2,1)
spec <- ugarchspec(
variance.model = list(model ="fGARCH",
garchOrder = c(2,1),
submodel = "TGARCH"),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "norm")
tgarch11 <- ugarchfit( spec = spec,
data = na.omit(portfolio3$r))Dla modelu TGARCH(2,1) wyznaczono VaR. Także zostały przedstawione zwroty wraz z oszacowaniami wartości narażonej na ryzyko dla modelu.
Dla modelu GARCH-t(1,1) straty przekroczyły zakładany poziom VaR w 0.0049% przypadków - 2 razy.
Na podstawie uzyskanych oszacowań VaR dla poszczególnych modeli można stwierdzić, że dla modelu TGARCH(2,1) straty przekroczyły zakładany poziom wartości narażonej na ryzyko najmniejszą ilość razy - 2. Więc można stwierdzić że model TGARCH(2,1) lepszy dopasowany.
O jakości dopasowania modelu można stwierdzać, kiedy porównujemy kryterium informacyjne poszczególnych modeli.
infocriteria(k2.garch11) #GARCH(1,1)##
## Akaike -6.311527
## Bayes -6.281702
## Shibata -6.311637
## Hannan-Quinn -6.299718infocriteria(k.garcht11) #GARCH-t(1,1)##
## Akaike -6.345643
## Bayes -6.305878
## Shibata -6.345839
## Hannan-Quinn -6.329899infocriteria(tgarch11) #TGARCH(2,1)##
## Akaike -6.374208
## Bayes -6.314560
## Shibata -6.374645
## Hannan-Quinn -6.350591Według kryterium informacyjnych model TGARCH(2,1) jest lepsze dopasowany.
Oszacowania wartości narażonej na ryzyko (VaR) dla wybranych modeli w okresie out-of-sample
W tym podrozdziale została przeprowadzona analiza porównawcza oszacowań wartości narażonej na ryzyko (VaR)w okresie out-of-sample, bazując na jednookresowych prognozach funkcji warunkowej wariancji.
Najpierw zmieniono próbę na okres out-of-sample:
start <- portfolio$obs[portfolio$Date == as.Date("2016-12-21")]
finish <- portfolio$obs[portfolio$Date == as.Date("2017-06-20")]
portfolio4 <- portfolio[start:finish, ]Następnym krokiem jest wykonanie pętli, która prognozuje wartości narażonej na ryzyko w okresie out-of-sample.
- GARCH(1,1)
VaR <- rep(NA, times = finish - start + 1)
for (k in start:finish) {
tmp.data <- portfolio[portfolio$obs <= (k - 1), ]
tmp.data <- tmp.data[as.Date("2011-12-20") <= 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 = "norm")
tmp.garch11 <- ugarchfit(spec = spec, data = na.omit(tmp.data$r))
sigma.forecast <- ugarchforecast(tmp.garch11, n.ahead = 1)
sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
VaR[k - start + 1] <- q01 * sigma.forecast2
}
# przypisujemy oszacowania VaRów do zbioru
portfolio4$VaR <- VaRNastępnie przedstawiono prognozę wartości narażonej na ryzyko(VaR) dla modeli GARCH(1,1) w okresie out-of-sample. Przedstawiono wykres prognozy:
plot(portfolio4$Date, portfolio4$r, col = "red", lwd = 1, type = 'l',
main = 'Zwroty i oszacowania VaR modelu GARCH(1,1) dla okresu in-sample', xlab = "Rok", ylab = "Wartość zwrotów",
ylim = c(-0.035, 0.025))
abline(h = 0, lty = 2)
lines(portfolio4$Date, portfolio4$VaR, type = 'l', col = "green")
W 0.016129% przypadkach straty przekroczyły zakładany poziom VaR.
- GARCH-t(1,1)
Wydruk wykonania pętli która prognozuje wartości narażonej na ryzyko w okresie out-of-sample dla modeli GARCH-t(1,1) został ukryty, ponieważ zmiany stosują się wyłącznie specyfikacji modeli.
Następnie przedstawiono prognozę wartości narażonej na ryzyko(VaR) dla modeli GARCH-t(1,1) w okresie out-of-sample. Przedstawiono wykres prognozy:
plot(portfolio4$Date, portfolio4$r, col = "red", lwd = 1, type = 'l',
main = 'Zwroty i oszacowania VaR modelu GARCH-t(1,1) dla okresu in-sample', xlab = "Rok", ylab = "Wartość zwrotów",
ylim = c(-0.035, 0.025))
abline(h = 0, lty = 2)
lines(portfolio4$Date, portfolio4$VaR_t, type = 'l', col = "green")
W 0.016129% przypadkach straty przekroczyły zakładany poziom VaR.
- TGARCH(2,1)
Następnie przedstawiono prognozę wartości narażonej na ryzyko(VaR) dla modeli TGARCH(2,1) w okresie out-of-sample. Przedstawiono wykres prognozy:
plot(portfolio4$Date, portfolio4$r, col = "red", lwd = 1, type = 'l',
main = 'Zwroty i oszacowania VaR modelu TGARCH(2,1) w okresie in-sample', xlab = "Rok", ylab = "Wartość zwrotów",
ylim = c(-0.035, 0.025))
abline(h = 0, lty = 2)
lines(portfolio4$Date, portfolio4$VaR_T, type = 'l', col = "green")
W 0.016129% przypadkach straty przekroczyły zakładany poziom VaR.
Analiza porównawcza prognoz wartości narażonej na ryzyko dla poszczególnych modeli w okresie out-sample wygląda następująco:
plot(portfolio4$Date, portfolio4$r, col = "red", lwd = 1, type = 'l',
main = 'Zwroty i oszacowania VaR modelej GARCH w okresie in-sample', xlab = "Rok", ylab = "Wartość zwrotów",
ylim = c(-0.035, 0.016))
abline(h = 0, lty = 2)
lines(portfolio4$Date, portfolio4$VaR, type = 'l', col = "yellow")
lines(portfolio4$Date, portfolio4$VaR_t, type = 'l', col = "blue")
lines(portfolio4$Date, portfolio4$VaR_T, type = 'l', col = "green")
legend("topright", legend=c("TGARCH(2,1)","GARCH-t(1,1)",
"GARCH(1,1)"), lwd=c(2,2), col=c("green","red","yellow"))
Z analizy graficznej widać że straty przekroczyły zakładany poziom VaR w dwóch miejscach. Prognoza wartości narażonej na ryzyko dla modeli GARCH(1,1) i GARCH-t przybliżono są podobne. Dla prognozy TGARCH kształt prognozy lepiej reaguje na wahanie się stop zwrotów, tak jak jest lepiej dopasowany.
Analiza wrażliwości estymacji parametrów modelu w zależności od długości próby
W tej częsci raportu zostały przedstawione wyniki prognozy wartości narażonej na ryzyko w zależności od długości przedziału out-of-sample. Do analizy wybrano wszystkie poprzednio zaznaczone modele GARCH. Wybrano okresy 3 i 12 miesięczne out-of-sample.
Prognoza VaR poszczególnych modeli dla 12 miesięcy:
plot(portfolio6$Date, portfolio6$r, col = "red", lwd = 1, type = 'l',
main = 'Zwroty i oszacowania VaR modelej GARCH w okresie in-sample(12 miesięcy)', xlab = "Rok", ylab = "Wartość zwrotów",
ylim = c(-0.035, 0.016))
abline(h = 0, lty = 2)
lines(portfolio6$Date, portfolio6$VaR, type = 'l', col = "yellow")
lines(portfolio6$Date, portfolio6$VaR_t, type = 'l', col = "blue")
lines(portfolio6$Date, portfolio6$VaR_T, type = 'l', col = "green")
legend("topright", legend=c("TGARCH(2,1)","GARCH-t(1,1)",
"GARCH(1,1)"), lwd=c(2,2), col=c("green","red","yellow"))
Prognoza VaR poszczególnych modeli dla 3 miesięcy:
plot(portfolio5$Date, portfolio5$r, col = "red", lwd = 1, type = 'l',
main = 'Zwroty i oszacowania VaR modelej GARCH w okresie in-sample(3 miesięcy)', xlab = "Rok", ylab = "Wartość zwrotów",
ylim = c(-0.035, 0.016))
abline(h = 0, lty = 2)
lines(portfolio5$Date, portfolio5$VaR, type = 'l', col = "yellow")
lines(portfolio5$Date, portfolio5$VaR_t, type = 'l', col = "blue")
lines(portfolio5$Date, portfolio5$VaR_T, type = 'l', col = "green")
legend("topright", legend=c("TGARCH(2,1)","GARCH-t(1,1)",
"GARCH(1,1)"), lwd=c(2,2), col=c("green","red","yellow"))
Podsumowanie
W tym badaniu przeprowadzono analizę porównawczą oszacowań warunkowej wariancji, oraz oszacowań wartości narażonej na ryzyko (VaR). Zostały porównywane modeji GARCH-t i TGARCH w odniesieniu do modelu bazowego GARCH(p, q) w okresie in-sample i out-of-sample. Na podstawie estymacji danych okazało się, że najlepiej dopasowanym modelem jest TGARCH. Jest to model który posiada najniższe wartości kryterium informacyjnych. Także z prognozowano wartości narażone na ryzyko dla okresu 3,6 i 12 miesięcznego out-of-sample. Podsumowując, takie prognozowania dają możliwość na określenie maksymalnego poziomu straty inwestora z analizowanego portfela aktywów.