Analiza porównawcza ryzyka portfela zrównoważonego oszacowanego w modelach klasy GARCH
Szymon Lesiak, Wojciech Mazur
30 czerwca 2017
Wstep
Dywersyfiacja portfela inwestycyjnego jest istotnym czynnikiem wpływającym na ograniczenie strat wynikłych z wahań rynku oraz szoków. Ponadto, przy odpowiednim doborze aktywów możliwy jest hedging ograniczający potencjalne straty. Odpowiednimi metodami do analizy stopnia ryzyka portfela jest oszacowanie jego warunkowej wariancji oraz wartości narażnej na ryzyko. Ponadto, analizy te mogą być przeprowadzane iteracyjnie wraz z napływem danych. Istotnym elementem analizy jest wybór odpwiedniej długości okresu in-sample i out-of sample. W przypadku pierwszej próby, dłuższy okres pozwala uzyskać więcej informacji o badanym zjawisku, jednak nie ma gwarancji, że informacja ta jest aktualna. Z kolei w przypadku okresu prognozy, apriori wiadomo, że im dłuższy tym mniej dokładna predykcja. W celu zbadania danego zagadnienia wykorzystane zostaną modele klasy Garch oraz Var estymowane na zrównoważonym portfelu inwestyctyjnym. Pierwsza klasa modeli pozwoli uzyskać oszaocwanie warunkowej wariancji potfela, z kolei druga umożliwi estymować wartość narażoną na ryzyko. Weryfikacją poprawności oraz stopania restrykcyjności modelu będzie empiryczne sprawdzenie częstości niedoszacowania warunkowej wariancji.
Metodologia
Pierwszym krokiem wykonanym w celu analizy była budowa portfela zrównoważonego składającego się z 5 instrumentów finansowych. Dane zostały w całości zaczerpnięte ze strony http://stooq.pl/ i obejmują dzienne ceny zamknięcia. W skład portfela zrównoważonego użytego w analizie wchodzą: notowania cen akcji APPLE w dolarach, ceny ropy BRENT w dolarach, notowania cen akcji PKN ORLEN w zł,ceny akcji PKO BP i notowania kursu walutowego UDS/PLN. Dane pochodzą z okresu od 2007-06-01 i zostały obcięte do 2017-06-01. Ze względu na braki danych usunięto 137 pustych obserwacji. Następnie dla utworzonego portfela został oszacowany model GARCH z dwoma rozszerzeniami: EGARCH (Exponential GARCH) i TGARCH (Treshold GARCH) Parametry modeli klasy GARCH zostały oszacowane na próbce in-sample obejmującej okres od 2007-06-01 do 2016-06-01 czyli 9 lat i na ich podstawie wykonana została prognoza jednookresowa na okres 1 roku tzw. out-of-sample. Dla oszacowanych modeli została przeprowadzona analiza porównawcza oszacowań warunkowej wariancji i oszacowań wartości narażonej na ryzyko (ang. Value-at-Risk, VaR) w okresie in-sample. Ostatnim część pracy zawiera identyczną analizę przeprowadzoną na okresie out-of-sample, bazującą na jednookresowych prognozach funkcji warunkowej wariancji.
Przed dokonaniem analizy wczytajmy przydatne biblioteki i funkcje.
Wymagane biblioteki:
library(xts)
library(fBasics) # e.g. basicStats()
library(tseries)# e.g. jarque.bera.test()
library(car) # e.g. durbinWatsonTest()
library(FinTS) # e.g. ArchTest()
library(fGarch) # e.g. garchFit()
library(rugarch) # e.g. ugarchfit()
library(plotrix)Wybór wskaznych aktywów do portfela spowodowany był chęcią holistycznego ujęcia rynku płynnych paliw kopalnych-stąd notowania ropy BRENT, Orlenu oraz kursu dolara (obrót ropą rowadzony jest właśnie w tej walucie). Oczekuje się, że dobór takch notowań pozwoli na ograniczenie wariancji poprzez hedging. Ponadto, do budowy portfela wykorzystano akcje Apple, jako przedstawiciela sektora technologicznego oraz Banku PKO z sektora usług finansowych. Ostatni wybór był spowodowany oczekiwanym sygnalingiem sektora bankowego w warunkach zwiększonej niepewności.
## 'data.frame': 2433 obs. of 6 variables:
## $ Date : Date, format: "2007-06-01" "2007-06-04" ...
## $ AAPL : num 15.3 15.7 15.8 16 16.1 ...
## $ OIL : num 69.1 70.3 70.4 71 68.5 ...
## $ PKN : num 42.7 43.2 44.2 44.3 45.4 ...
## $ PKO : num 173 175 173 166 165 ...
## $ USDPLN: num 2.83 2.82 2.83 2.84 2.87 ...## Date AAPL OIL PKN PKO USDPLN
## 1 2007-06-01 15.277 69.09 42.656 172.69 2.8281
## 2 2007-06-04 15.655 70.30 43.162 174.95 2.8190
## 3 2007-06-05 15.828 70.39 44.176 173.40 2.8263
## 4 2007-06-06 15.952 70.96 44.260 166.41 2.8373
## 5 2007-06-08 16.064 68.50 45.359 165.44 2.8679
## 6 2007-06-11 15.508 69.24 45.400 169.46 2.8696
## 7 2007-06-12 15.533 68.73 43.415 164.91 2.8755
## 8 2007-06-13 15.161 69.84 44.346 165.57 2.8754
## 9 2007-06-14 15.323 71.28 44.766 166.22 2.8693
## 10 2007-06-15 15.547 71.47 45.570 171.06 2.8402## Date AAPL OIL PKN PKO USDPLN
## 2424 2017-05-18 152.54 52.49 109.70 131.12 3.80006
## 2425 2017-05-19 153.06 53.77 110.05 130.19 3.74324
## 2426 2017-05-22 153.99 53.74 110.90 130.84 3.73495
## 2427 2017-05-23 153.80 54.19 106.05 130.19 3.75487
## 2428 2017-05-24 153.34 53.85 105.90 130.65 3.72254
## 2429 2017-05-25 153.87 51.33 108.50 131.87 3.72401
## 2430 2017-05-26 153.61 52.26 105.10 130.28 3.74550
## 2431 2017-05-30 153.67 52.19 107.90 128.18 3.73300
## 2432 2017-05-31 152.76 50.97 106.15 128.79 3.71934
## 2433 2017-06-01 153.18 50.38 104.60 126.92 3.73905Wstepna analiza danych
Po wczytaniu odpowiednich danych i utworzeniu z nich portfela należy prześledzić przebieg zmienności notowań w czasie. W tym celu zaprezentowane zostaną wykresy z cenami kilku noto
twoord.plot(lx=c(seq(1:2433)),ly=portfel$USDPLN,rx=c(seq(1:2433)),
ry=portfel$OIL,
xlab="Date",
ylab="USDPLN",rylab="OIL",lcol=4,
main="Wykres cen ropy i kursu Dolara",
type=c("l","l"),
xticklab = portfel$Date)
Zgodnie z oczekiwaniami cena ropy jest skorelowana kursem dolara. Ujęcie w portfelu wskazanych instumentów pozwala na ograniczenie strat wynikających ze spadku cen jednego z aktywów. Jednak tym samym ogranicza potencjalne zyski w przypadku wzrostu jednego z nich. Jednak w strategii długookresowej wskazana relacja może okazać się elementem stabiliuzjącym wariancję portfela
twoord.plot(lx=c(seq(1:2433)),ly=portfel$AAPL,rx=c(seq(1:2433)),ry=portfel$OIL,
xlab="Date",
ylab="AAPL",rylab="OIL",lcol=4,
main="Wykres cen ropy i notowań Apple",
type=c("l","l"),
xticklab = portfel$Date)
W przypadku notowań Apple utrzymuje się rosnący trend, lecz zauważalne jest podobne zachowanie indeksów w pewnych okresach, np. po kryzysie finansowym do 2012 roku.
twoord.plot(lx=c(seq(1:2433)),ly=portfel$PKN,rx=c(seq(1:2433)),ry=portfel$OIL,
xlab="Date",
ylab="PKN",rylab="OIL",lcol=4,
main="Wykres cen ropy i notowań Orlen",
type=c("l","l"),
xticklab = portfel$Date)
Natomiast w przypadku akcji PPKN Orlen zauważalne jest, że tendencja wzrostowa jest silniejsza w warunkach taniego surowca. Mechanizm ten zgadza się z intuicją, ponieważ niższe ceny ropy naftowej to szok popytowy dla całej gospodarki,w tym sektora rafineryjnego.
twoord.plot(lx=c(seq(1:2433)),ly=portfel$USDPLN,rx=c(seq(1:2433)),ry=portfel$PKN,
xlab="Date",
ylab="USDPLN",rylab="OIL",lcol=4,
main="Wykres cen ropy i notowań Orlen",
type=c("l","l"),
xticklab = portfel$Date)
Przypuszczać można, że niższy kurs dolara również będzie stanowił pozytywny szok popytowy dla sektora rafinaryjnego, jednak hipotezy tej nie można odrzucać ani przyjmować na podstawie analizowanego okresu. Kurs dolara oraz cena ropy Brent utrzymują zbliżony trend. Znaczące fluktuacje kursu do roku 2010 nie oddziałują na notowania Orlenu, z kolei szybki wzrost cen akcji polskiej społki jest spowodowany spadkiem cen ropy. Jednak, jak już wcześniej wspomniano, istnieje zależność między kursem dolara a ceną ropty Brent, która to z kolei wpływa na wyniki Orlenu.
twoord.plot(lx=c(seq(1:2433)),ly=portfel$PKN,rx=c(seq(1:2433)),
ry=portfel$PKO,
xlab="Date",
ylab="PKN",rylab="PKO",lcol=4,
main="Wykres notowań PKO i Orlen",
type=c("l","l"),
xticklab = portfel$Date)
Kolejnym analizowanym aktywem są akcje banku PKO. Powodem uwzględnienia notowań banku w badaniu jest przypuszczenie, że zwiększona niepewność na rynku będzie silniej oddziaływać na akcje instytucji finansowych. To przypuszczenie potwierdza się, ponieważ zauważalne są dwa znaczące topnięcia przed kryzysem finansowym oraz w 2015 roku, czyli w czasie ostaniej “fali” kryzysu. Skala spadków cen akcji jest większa niż w przypadku Orelnu i obserwowana wcześniej. Obserwacja ta jest ważna z perspektywy monitoringu portfela.
twoord.stackplot(lx=portfel$Date,rx=portfel$Date,ldata=portfel$USDPLN,
rdata=cbind(portfel$OIL, portfel$PKN,portfel$PKO,portfel$AAPL),
lcol=c("red"),
rcol=c("black", "green","blue","yellow"), ltype="l", rtype=c("l","l"),
border="grey80",xlab="", lylab="USDPLN", rylab="PKN,OIL,PKO,AAPL",
main="Wykres notowań dolara,akcji Orlen oraz cen ropy Brent")
par(xpd=TRUE) #extend the area of plotting
par(new=TRUE) #to add new graph "layers"
plot(0:1, 0:1, type="n", xlab="",ylab="", axes=FALSE) #redo the x/y limits
#first legend
legend(-0.18, 1.2, leg=c("USDPLN","OIL","PKN","PKO","AAPL"), fill=c("red", "black","green","blue","yellow"))
Podsumowując, zaprezentowane instrumenty tworzą portfel, który w pewnym stopniu można hedgować, co ogranicza zmienność jego wariancji,a zarazem odpowiednio zdywersyfikowany, co ogranicza ryzyko sektorowe . Z perspektywy niniejszego badania sa to istotne elementy, zatem wybrane aktywa są poprawne pod względem metodologicznym.
Logarytmiczne stopy zwrotu
Teraz zostaną utworzone logarytmiczne (ciągłe) stopy zwrotu
portfel$r1 <- diff.xts(log(portfel$AAPL))
portfel$r2 <- diff.xts(log(portfel$OIL))
portfel$r3 <- diff.xts(log(portfel$PKN))
portfel$r4 <- diff.xts(log(portfel$PKO))
portfel$r5 <- diff.xts(log(portfel$USDPLN))Ponieważ portfel równoważony jest tworzony w ten sposób, że każdemu instrumentowi przypisuje taka samą wagę zadeklarujemy nową zmienną, która stanowi średnią arytmetyczną stóp zwrotu wszystkich 5 szeregów
portfel$r <- (portfel$r1+portfel$r2+portfel$r3+portfel$r4+portfel$r5)/5Wykres logarytmicznych stóp zwrotu portfela prezentuje sie następująco:
Zauważalny jest silny efekt kryzysu finansowego na zmienność stopy zwrotu. Ponadto, występuje zjawisko grupowania wariancji. Jednak powstały portfel spełnia zadaną funkcję ograniczając zmienność o czym świadzą wartości wariancji ciągłych stóp z próby dla poszczególnych aktywów.
## varr varr1 varr2 varr3 varr4
## [1,] 0.0001250657 0.0004417526 0.0004975646 0.0005117123 0.0005433937
## varr5
## [1,] 0.0001045786Wariancja ciągłej stopy zwrotu portfela jest znacznie mniejsza od wszystkich ciągłych stóp zwrotu aktywów, oprócz kursu dolara (najmniejsza jednostka miary). Kolejnym etapem analizy jest zbadanie, czy zmienna podlega autokorelacji oraz charakterysty jej rozkładu.
Wykres wartości ACF dla zwrotów:
Korelogram ciągłych zwrotów wskazuje, że badana zmienna podlega autokorelacji. Bardziej widoczny efekt autokorelacji widoczny będzie dla kwadratów zwrotów.
Korelogram kwadratów zwrotów z portfela potwierdza wcześniej wysnuty wniosek. Ciągła stopa zwrotu podlega autokorelacji.
Teraz zostanie przeprowadzony test Durbin’a-Watson’a dla reszt i kwadratów reszt stóp zwrotu portfela. Wyniki testu prezentują się następująco:
Reszty:
## lag Autocorrelation D-W Statistic p-value
## 1 0.036902269 1.925533 0.072
## 2 -0.093316477 2.185589 0.000
## 3 0.002639894 1.993617 0.904
## 4 0.005004943 1.988817 0.826
## 5 -0.014408053 2.027625 0.468
## Alternative hypothesis: rho[lag] != 0Z wyników testu możemy zaobserwować, że p-value dla wszystkich opóźnień z wyjątkiem drugiego przekracza założony poziom istotności 5%, co oznacza, że nie mamy podstaw do odrzucenia hipotezy zerowej o braku autokorelacji reszt dla 1,3,4 i 5 opóźnienia. Wartości statystyk Durbin’a-Watson’a także są bliskie wartości 2.
Kwadraty reszt:
## lag Autocorrelation D-W Statistic p-value
## 1 0.1228113 1.754349 0.000
## 2 0.1860194 1.627896 0.000
## 3 0.1098358 1.780166 0.002
## 4 0.1648246 1.670093 0.002
## 5 0.1479647 1.703705 0.004
## Alternative hypothesis: rho[lag] != 0Dla kwadratów reszt możemy zaobserwować bardzo niskie wartości p-value. Oznacza to, że dla każdego opóźnienia odrzucamy hipotezę zerową o braku autokorelacji. Zatem występuje autokorelacja kwadratów reszt stóp dla 5 pierwszych opóźnień. Wnioski wyciągniete z analizy korelogramów kwadratów reszt znajdują więc potwierdzenie w testach.
W dalszej części pracy przydatną informacją okażą się własności rozkładu badanej zmiennej, dlatego poniżej zostaną predstawione oraz omówione odpowiednie statystyki oraz testy.
Statystyki opisowe:
## X..portfel.r
## nobs 2433.000000
## NAs 1.000000
## Minimum -0.092131
## Maximum 0.067130
## 1. Quartile -0.005766
## 3. Quartile 0.005942
## Mean 0.000235
## Median 0.000378
## Sum 0.571542
## SE Mean 0.000227
## LCL Mean -0.000210
## UCL Mean 0.000680
## Variance 0.000125
## Stdev 0.011183
## Skewness -0.264806
## Kurtosis 5.393258Ze względu na to, że badana zmienna jest charakterystyką finansową, najistotniejszymi obserwacjami na temat rozkładu są skośność oraz nadwyżkowa kurtoza. Dla ciągłej stopy zwrotu z analizowanego portfela wynoszą one odpowiednio -0.265806 oraz 5.393258 i różnią się od wartości przyjmowanych przez zmianną z rozkładu normalnego( w obu przypadkach 0). Ponadto, nadwyżkowa kurtoza świadczy o leptokurtozie rozkładu badanej zmiennej. Analizę rozkładu może wspomóc również jej histogram.
Histogram wzbogacony o funkcję gestości rozkładu normalnego z parametrami oszacowanymi na podstawie pełnej próbki zwrotów:
Histogram badanej zmiennej potwierdza, że pochodzi z leptokurtycznego rozkładu, o czym świadczy nadmiarowa koncentracja obserwacji wokół średniej oraz grubsze ogony w porównaniu do rozkładu normalnego. Ostatecznym potwierdzeniem, czy zmienna pochodzi z rozkładu normalnego jest test Jarque-Bery.
Test Jarque-Bery dla zwrotów:
##
## Jarque Bera Test
##
## data: na.omit(portfel$r)
## X-squared = 2983.5, df = 2, p-value < 2.2e-16Na podstawie testu należy odrzucić hipotezę zerową o rozkladzie normalnym zmiennej (statystyka testowa=2983.5, p-value=0).
Wyniki testu na występowanie efektów ARCH wśród 5 pierwszych opóźnień zostaną zaprezentowane poniżej:
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: portfel$r
## Chi-squared = 176.9, df = 5, p-value < 2.2e-16Bardzo niska wartość p-value zmusza nas do odrzucenia hipotezy zerowej o braku efektów ARCH. Przyjmujemy zatem hipotezę alternatywną, że efekty te występują.
Poszukiwania parametrów modelu klasy GARCH(p,q)
Przed przystąpieniem do dopasowania odpowiedniego modelu do danych próbka została podzielona na okresy: in-sample od 2007-06-01 do 2016-06-01 i out-od-sample od 2016-06-02 do 2017-06-01
W pierszej kolejności został oszacowany model GARCH(1,0) czyli zwykły ARCH(1). Wyniki estymacji tego modelu prezentują się następująco:
## [1] " Estimate Std. Error t value Pr(>|t|) "
## [2] "mu 3.281e-04 2.339e-04 1.403 0.161 "
## [3] "omega 1.031e-04 4.211e-06 24.497 < 2e-16 ***"
## [4] "alpha1 2.440e-01 3.803e-02 6.415 1.41e-10 ***"P-value testu t dla zmiennej “mu” przekracza założony poziom istotności 5%. Oznacza to, że stała w równaniu średniej jest nieistotna. Oszacowany zostanie zatem model bez stałej:
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 0), data = na.omit(insample$r),
## cond.dist = "norm", include.mean = F, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 0)
## <environment: 0x0000000019d40260>
## [data = na.omit(insample$r)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## omega alpha1
## 0.00010369 0.23817384
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 1.037e-04 4.196e-06 24.708 < 2e-16 ***
## alpha1 2.382e-01 3.720e-02 6.402 1.53e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 6713.186 normalized: 3.069587
##
## Description:
## Thu Jul 06 21:38:37 2017 by user: wojtu
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2242.83 0
## Shapiro-Wilk Test R W 0.9554431 0
## Ljung-Box Test R Q(10) 33.98332 0.0001859032
## Ljung-Box Test R Q(15) 43.5209 0.0001306483
## Ljung-Box Test R Q(20) 45.52206 0.0009370152
## Ljung-Box Test R^2 Q(10) 240.2708 0
## Ljung-Box Test R^2 Q(15) 325.1399 0
## Ljung-Box Test R^2 Q(20) 420.4945 0
## LM Arch Test R TR^2 173.8059 0
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.137344 -6.132140 -6.137346 -6.135442P-value testu Ljunga-Boxa dla reszt i kwadratów reszt dla pierwszych 10, 15 i 20 opóźnień jest bardzo niskie więc należy odrzucić hipotezę zerową o braku autokorelacji. Należy także odrzucić hipotezę zerową o braku efektów ARCH, gdyż p-value testu LM Arch Test wynosi 0. Ponieważ model ARCH(1) nie wyjaśnia efektów ARCH i występuję autokorelacja reszt i kwadratów reszt oszacowany zostanie model ARCH(3) czyli tzw. GARCH(3,0).
Test Ljunga-Boxa: H0: pierwsze q wsp. autokorelacji jest rowne 0
H1: wystepuje autokorelacjaWyniki estymacji modelu ARCH(3):
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(3, 0), data = na.omit(insample$r),
## cond.dist = "norm", include.mean = F, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(3, 0)
## <environment: 0x0000000021b7bc98>
## [data = na.omit(insample$r)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## omega alpha1 alpha2 alpha3
## 6.8268e-05 1.9219e-01 1.5085e-01 1.8151e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 6.827e-05 4.222e-06 16.171 < 2e-16 ***
## alpha1 1.922e-01 3.578e-02 5.372 7.80e-08 ***
## alpha2 1.509e-01 2.877e-02 5.243 1.58e-07 ***
## alpha3 1.815e-01 3.272e-02 5.548 2.89e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 6791.575 normalized: 3.10543
##
## Description:
## Thu Jul 06 21:38:37 2017 by user: wojtu
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2233.458 0
## Shapiro-Wilk Test R W 0.9654204 0
## Ljung-Box Test R Q(10) 21.36705 0.01867458
## Ljung-Box Test R Q(15) 26.23592 0.03561489
## Ljung-Box Test R Q(20) 27.70394 0.1165951
## Ljung-Box Test R^2 Q(10) 57.43166 1.105835e-08
## Ljung-Box Test R^2 Q(15) 81.57842 3.589429e-11
## Ljung-Box Test R^2 Q(20) 107.4596 5.695444e-14
## LM Arch Test R TR^2 75.10259 3.5137e-11
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.207202 -6.196794 -6.207208 -6.203398Na podstawie odpowiednich testów diagnostycznych można smiało stwierdzić, że w modelu tym także występuje autokorelacja reszt dla pierwszych 15 opóźnień i autokorelacja kwadratów reszt dla pierwszych 20 opóźnień. Model ARCH(3) nie radzi sobie także z wytłumaczeniem efektów ARCH.
Teraz zostanie oszacowany model GARCH(1,1). W modelu tym stała w równaniu średniej jest istotna jednak zostanie zaprezentowany model bez stałej, który osiąga niższa wartości kryterium BIC. Wyniki estymacji poniżej.
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = na.omit(insample$r),
## cond.dist = "norm", include.mean = F, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x000000001e15cc60>
## [data = na.omit(insample$r)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## omega alpha1 beta1
## 8.8477e-07 5.5385e-02 9.3788e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 8.848e-07 3.252e-07 2.721 0.00651 **
## alpha1 5.538e-02 8.191e-03 6.762 1.36e-11 ***
## beta1 9.379e-01 8.800e-03 106.582 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 6947.623 normalized: 3.176782
##
## Description:
## Thu Jul 06 21:38:38 2017 by user: wojtu
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 119.8954 0
## Shapiro-Wilk Test R W 0.991538 5.547527e-10
## Ljung-Box Test R Q(10) 11.29215 0.3352143
## Ljung-Box Test R Q(15) 13.56268 0.5589208
## Ljung-Box Test R Q(20) 14.32762 0.8135038
## Ljung-Box Test R^2 Q(10) 16.76817 0.07965367
## Ljung-Box Test R^2 Q(15) 19.06528 0.210791
## Ljung-Box Test R^2 Q(20) 21.44264 0.3715034
## LM Arch Test R TR^2 17.28158 0.1393086
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.350822 -6.343016 -6.350825 -6.347968Wartość p-value testu Ljunga-Boxa przekracza poziom istotności 5% co oznacza, że nie mamy podstaw do odrzucenia hipotezy zerowej o braku autokorelacji reszt i kwadratóW reszt 10, 15 i 20 opóźnień. Z testu LM Arch Test wynika także, że nie mamy podstaw do odrzucenia hipotezy zerowej o występowaniu efektów ARCH. Model GARCH(1,1) wytłumaczył zatem efekty arch. Parametry modelu są dodatnie, co oznacza, że spełnione jest założenie o dodatniej wariancji. Proces jest także stacjonarny ponieważ suma parametrów jest mniejsza od jedności.
Postanowiono także sprawdzić modele GARCH(2,1) i GARCH(1,2) jednak jak można zaobserwować dodatkowe parametry tj. odpowiednio beta2 i alpha2 są nieistotne
GARCH(1,2)
## Estimate Std. Error t value Pr(>|t|)
## omega 1.069090e-06 4.251084e-07 2.514866 1.190777e-02
## alpha1 7.060045e-02 1.501109e-02 4.703220 2.560899e-06
## beta1 6.349800e-01 2.286994e-01 2.776482 5.495062e-03
## beta2 2.864458e-01 2.169742e-01 1.320183 1.867738e-01Zmienna beta2 nieistotna - model niepoprawny.
GARCH(2,1)
## Estimate Std. Error t value Pr(>|t|)
## omega 8.846178e-07 3.456821e-07 2.559050e+00 0.01049587
## alpha1 5.542114e-02 2.243002e-02 2.470846e+00 0.01347937
## alpha2 1.000000e-08 2.422502e-02 4.127963e-07 0.99999967
## beta1 9.378579e-01 9.779341e-03 9.590195e+01 0.00000000Alpha2 nieistotna.
Szacowane modele zakładały, że stopy zwrotu pochodzą z rozkładu normalnego. Postanowiono sprawdzić model GARCH(1,1) bez stałej który zadkłada, że zwroty pochodzą z rozkładu t-studenta. Wyniki estymacji prezentują się następująco:
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = na.omit(insample$r),
## cond.dist = "std", include.mean = F, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x000000001f981038>
## [data = na.omit(insample$r)]
##
## Conditional Distribution:
## std
##
## Coefficient(s):
## omega alpha1 beta1 shape
## 1.2042e-06 5.4940e-02 9.3495e-01 8.3041e+00
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 1.204e-06 4.620e-07 2.606 0.00915 **
## alpha1 5.494e-02 1.066e-02 5.155 2.54e-07 ***
## beta1 9.349e-01 1.232e-02 75.914 < 2e-16 ***
## shape 8.304e+00 1.376e+00 6.034 1.60e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 6974.763 normalized: 3.189192
##
## Description:
## Thu Jul 06 21:38:39 2017 by user: wojtu
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 128.4862 0
## Shapiro-Wilk Test R W 0.9912992 3.550646e-10
## Ljung-Box Test R Q(10) 11.93822 0.2892128
## Ljung-Box Test R Q(15) 14.23329 0.5079142
## Ljung-Box Test R Q(20) 14.92459 0.7807062
## Ljung-Box Test R^2 Q(10) 17.56091 0.06283776
## Ljung-Box Test R^2 Q(15) 19.35233 0.1982096
## Ljung-Box Test R^2 Q(20) 21.51878 0.3671831
## LM Arch Test R TR^2 17.82834 0.1210028
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.374726 -6.364318 -6.374733 -6.370922W modelu zakładającym rozkład t-studenta także brak autokorelacji standaryzowanych reszt i ich kwadratóW dla 10, 15 i 20 pierwszych opóźnień. Brak efektów ARCH, co oznacza, że zostały one wytłumaczone przez model.
Warto jeszcze zobaczyć wartości kryteriów informacyjnych dla wyestymowanych modeli:
## AIC BIC SIC HQIC model
## 1 -6.137344 -6.132140 -6.137346 -6.135442 k.arch1
## 2 -6.207202 -6.196794 -6.207208 -6.203398 k.arch3
## 3 -6.350822 -6.343016 -6.350825 -6.347968 k.garch11
## 4 -6.352097 -6.341689 -6.352103 -6.348293 k.garch11stala
## 5 -6.374726 -6.364318 -6.374733 -6.370922 k.garch11student
## 6 -6.350365 -6.339957 -6.350371 -6.346560 k.garch12
## 7 -6.349927 -6.339520 -6.349934 -6.346123 k.garch21Wszystkie kryteria informacyjne są najmniejsze dla modelu GARCH(1,1) bez stałej zakładającego, że zwroty z portfela pochodzą z rozkładu t-studenta. Postanowiono wybrać ten model za ostateczny.
Zobaczmy ACF wystandaryzowanych reszt i ich kwadratow dla wybranego modelu
plot(k.garch11student, which=10)
plot(k.garch11student, which=11)
Wypustki wystandaryzowanych reszt i kwadratóW reszt są nieistotne gdyż nie wychodzą poza obszar obszar krytyczny. Oznacza to, że nie występuje autoorelacja reszt i kwadratów reszt. Potwierdzają to także testy Ljung-Box. Wypustki kwadratow wystandaryzowanych reszt są nieistotne więc nie występują efekty arch co także potwierdzaja testy.
Dla ostatecznie wybranego modelu GARCH (1,1):
Na wykresie można zaobserwować, że największe wzrosty warunkowej wariancji pojawiają się w miejscach tzw. grupowania wariancji zwrotów portfela szczególnie w okresie kryzysu.
Następnie zostanie sprawdzone czy wystandaryzowane reszty w modelu GARCH(1,1) pochodzą z rozkładu normalnego. Histogram z funkcją gętości rozkładu normalnego prezentuje się następująco:
Warto także przeprowadzić test Jarque-Bery:
# test Jarque-Bery
jarque.bera.test(stdres)##
## Jarque Bera Test
##
## data: stdres
## X-squared = 119.9, df = 2, p-value < 2.2e-16Wystandaryzowane reszty na histogramie przypominają w pewnym stopniu rozkład normalny, jednak wartość p-value testu Jarque-Bery jest bliska 0, co oznacza, że należy odrzucić hipotezę zerową o rozkładzie normalnym.
Test Durbina-Watsona na autokorelację kwadratów wystandaryzowanych reszt:
# test Durbina-Watsona
durbinWatsonTest(lm(formula = stdres^2 ~ 1),
max.lag = 5) # dla 5 pierwszych opóźnień## lag Autocorrelation D-W Statistic p-value
## 1 0.013128060 1.973588 0.520
## 2 -0.014434809 2.028614 0.484
## 3 -0.015395622 2.030273 0.458
## 4 0.045084049 1.909145 0.066
## 5 -0.003959422 2.007039 0.872
## Alternative hypothesis: rho[lag] != 0# DW oscyluje wokol 2 wiec chyba okDW oscyluje wokół wartości 2. P-value także wysokie, co oznacza, że nie mamy podstaw do odrzucenia hipotezy zerowej o braku autokorelacji kwadratów wystandaryzowanych reszt wśród wszystkich opóźnień.
Podsumowując model GARCH(1,1) bez stałej zakładający, że zwroty pochodzą z rozkładu t-studenta okazał się być najlepszy. Jest stacjonarny gdyż suma parametrów alpha i beta jest mniejsza od jedności. Parametry modelu są istotne. Model wyjaśnia efekty ARCH i nie występuje autokorelacja wystandaryzowanych reszt i ich kwadratów. Wystandaryzowane reszty nie pochodzą jednak z rozkładu normalnego. Jest to jednak częsty przypadek dla danych finansowych.
Rozszerzenia modeli klasy GARCH
EGARCH
Jak już wspomniano, klasyczny GARCH jest metodą pozwalającą wyeliminować problem leptokurtyczności oraz grupowania wariancji. Z kolei model Exponential Garch umożliwia uwzględnienie w modelu efektu dźwigni, powstałego na skutek niesymetrycznej reakcji na szoki. Analiza rozpoczęta będzie od modelu EGARCH(1,1), czyli podstawowego modelu tej klasy.
## [1] "Optimal Parameters"
## [2] "------------------------------------"
## [3] " Estimate Std. Error t value Pr(>|t|)"
## [4] "omega -0.051602 0.000889 -58.045 0"
## [5] "alpha1 -0.075096 0.006455 -11.635 0"
## [6] "beta1 0.994124 0.000031 31929.470 0"
## [7] "gamma1 0.074609 0.005510 13.540 0"
## [8] ""
## [9] "Robust Standard Errors:"
## [10] " Estimate Std. Error t value Pr(>|t|)"
## [11] "omega -0.051602 0.001313 -39.292 0"
## [12] "alpha1 -0.075096 0.006777 -11.081 0"
## [13] "beta1 0.994124 0.000047 21007.653 0"
## [14] "gamma1 0.074609 0.006182 12.068 0"
## [15] ""
## [16] "LogLikelihood : 6981.559 "
## [17] ""
## [18] "Information Criteria"
## [19] "------------------------------------"
## [20] " "
## [21] "Akaike -6.3809"
## [22] "Bayes -6.3705"
## [23] "Shibata -6.3809"
## [24] "Hannan-Quinn -6.3771"
## [25] ""
## [26] "Weighted Ljung-Box Test on Standardized Residuals"
## [27] "------------------------------------"
## [28] " statistic p-value"
## [29] "Lag[1] 2.463 0.11658"
## [30] "Lag[2*(p+q)+(p+q)-1][2] 5.469 0.03047"
## [31] "Lag[4*(p+q)+(p+q)-1][5] 7.461 0.04004"
## [32] "d.o.f=0"
## [33] "H0 : No serial correlation"
## [34] ""
## [35] "Weighted Ljung-Box Test on Standardized Squared Residuals"
## [36] "------------------------------------"
## [37] " statistic p-value"
## [38] "Lag[1] 0.4034 0.5254"
## [39] "Lag[2*(p+q)+(p+q)-1][5] 2.5941 0.4862"
## [40] "Lag[4*(p+q)+(p+q)-1][9] 4.6910 0.4762"
## [41] "d.o.f=2"
## [42] ""
## [43] "Weighted ARCH LM Tests"
## [44] "------------------------------------"
## [45] " Statistic Shape Scale P-Value"
## [46] "ARCH Lag[3] 0.1798 0.500 2.000 0.6715"
## [47] "ARCH Lag[5] 3.4364 1.440 1.667 0.2324"
## [48] "ARCH Lag[7] 3.8246 2.315 1.543 0.3726"Wszystkie zmienne w modelu są istotne na poziomie istotności 5%, zatem nie ma potrzeby ograniczać go. Ponadto, parametr alpha1 jest ujemny, co świadczy o występowaniu efektu dźwigni. Ponadto, na podstawie testu Ljung-Box nie ma podstaw do odrzucenia hipotezy o braku autokorelacji standaryzowanych reszt tylko dla jednego opóźnienia (0.1166>0.05). z kolei w przypadku tego samego testu dla kwadratów standaryzowanych reszt nie ma podstaw do odrzucenia hipotezy zerowej o braku autokorelacji pierwszego, drugiego, trzecigo rzędu. Ponadto, nie ma podstaw do odrzucenia hipotezy zerowej o braku efektów ARCH w resztach modelu dla opóźnień do siódmego rzędu (p-value dla siedmiu opóźnień=0.3724>0.05). Oszacowany model jest poprawny, ale należy sprawdzić, czy dodanie nowych parametrów nie poporawi dopasowania. Przed wyborem propozycji modelu, pomocna może okazać się analiza korelogramów reszt.
plot(k.egarch11, which=11)
Jedyna istotna utokorelacja reszt drugiego rzędu sugeruje, że wynika ona z procesu średniej ruchomej, dlatego w dalszym etapie analizy oszacowanie zostanie model EGARCH(1,2):
##
## *---------------------------------*
## * 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.069181 0.007046 -9.8187 0
## alpha1 -0.109063 0.011140 -9.7901 0
## beta1 0.518573 0.000043 12162.9342 0
## beta2 0.473524 0.000103 4608.3093 0
## gamma1 0.105999 0.018538 5.7179 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega -0.069181 0.004456 -15.5260 0e+00
## alpha1 -0.109063 0.015193 -7.1783 0e+00
## beta1 0.518573 0.000198 2624.2532 0e+00
## beta2 0.473524 0.000181 2616.0418 0e+00
## gamma1 0.105999 0.022818 4.6455 3e-06
##
## LogLikelihood : 6984.292
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.3825
## Bayes -6.3695
## Shibata -6.3825
## Hannan-Quinn -6.3778
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.336 0.12642
## Lag[2*(p+q)+(p+q)-1][2] 5.398 0.03183
## Lag[4*(p+q)+(p+q)-1][5] 7.431 0.04072
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.006322 0.9366
## Lag[2*(p+q)+(p+q)-1][8] 3.763791 0.5471
## Lag[4*(p+q)+(p+q)-1][14] 8.476549 0.3324
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 4.981 0.500 2.000 0.02563
## ARCH Lag[6] 5.216 1.461 1.711 0.10140
## ARCH Lag[8] 5.329 2.368 1.583 0.21589
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.2095
## Individual Statistics:
## omega 0.83642
## alpha1 0.05258
## beta1 0.82622
## beta2 0.82615
## gamma1 0.42523
##
## 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 0.4656 0.6416
## Negative Sign Bias 0.2113 0.8326
## Positive Sign Bias 1.0633 0.2877
## Joint Effect 2.8369 0.4175
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 22.68 0.2515
## 2 30 32.44 0.3010
## 3 40 36.19 0.5987
## 4 50 49.01 0.4728
##
##
## Elapsed time : 0.5067449Model EGARCH(2,1), pomimo istotnych wszystkich zmiennych nie wydaje się być istotnie lepszym od EGARCH(1,1). Świadczy o tym wyższa wartośc kryterium BiC. Wygenerowano również korelogram stnadaryzowanych reszt:
plot(k.egarch12, which=11)
W GARCH(1,2) nie występuje korelacja drugiego rzędu. Jednak istotne są korelacje czwarta i dziewiąta. Jak widać, wystąpiło przesunięcie niewyjaśnionej zmienności na dalsze korelacje. Na podstawie wskazanych argumentów lepszym modelem jest EGARCH(1,1). Ostatnią propozycją polepszenia modelu, jest zmiana zakładanego rozkładu ciągłej stopy zwrotu z normalnego na t-studenta:
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega -0.062157 0.001188 -52.3346 0
## alpha1 -0.070647 0.009148 -7.7229 0
## beta1 0.993164 0.000048 20494.2798 0
## gamma1 0.075748 0.003864 19.6051 0
## shape 10.526203 2.025624 5.1965 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega -0.062157 0.001386 -44.8451 0
## alpha1 -0.070647 0.010248 -6.8941 0
## beta1 0.993164 0.000052 18990.1159 0
## gamma1 0.075748 0.005438 13.9296 0
## shape 10.526203 2.003656 5.2535 0
##
## LogLikelihood : 6997.837
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.3949
## Bayes -6.3819
## Shibata -6.3949
## Hannan-Quinn -6.3902
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.553 0.11009
## Lag[2*(p+q)+(p+q)-1][2] 5.623 0.02773
## Lag[4*(p+q)+(p+q)-1][5] 7.652 0.03594
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.3003 0.5837
## Lag[2*(p+q)+(p+q)-1][5] 2.4877 0.5085
## Lag[4*(p+q)+(p+q)-1][9] 4.5997 0.4905
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1833 0.500 2.000 0.6685
## ARCH Lag[5] 3.4466 1.440 1.667 0.2312
## ARCH Lag[7] 3.8488 2.315 1.543 0.3689
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.2196
## Individual Statistics:
## omega 0.77789
## alpha1 0.08563
## beta1 0.75517
## gamma1 0.01606
## shape 0.19629
##
## 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 0.6507 0.5153
## Negative Sign Bias 0.7076 0.4792
## Positive Sign Bias 0.9193 0.3580
## Joint Effect 5.7395 0.1250
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 19.32 0.4365
## 2 30 19.98 0.8935
## 3 40 26.79 0.9309
## 4 50 36.43 0.9080
##
##
## Elapsed time : 1.344542Wartość kryteriów informacyjnych jest lepsza od modelu EGARCH(1,1) zakładającego rozład normalny zwrotów. Ostatecznym modelem tej klasy jest EGARCH(1,1) zakładający stopy zwrotu pochodzące z rozkładu t-studenta. Jak już wspomniano, w modelu występuje efekt dźwigni. Warto zobaczyć wizualizację tego zjawiska:
plot(k.egarchs11, which=12)
Dla badanej zmiennej asymetryczna reakcja na szoki wygląda bardzo ciekawie. Na podstawie analizowanego modelu można stwierdzić, że warunkowa wiariancja silniej reaguje na spadki cen instumentów, natomiast w minimalnym stopniu na wzrosty. Niestety taka cecha z perspektywy inwestora jest niekorzystna, ponieważ portfel nie jest odporny na niepwność na rynku.
Wzrost warunkowej wariancji obserwowany jest w czasie wzrostu niepewności na rynku. Największe wartości przyjmuje w czasie kryzysu.
Kolejnym wybranym na podstawie analizy rozszerzeniem GARCH jest model Threshold GARCH. W przypadku tego modelu również można uwzględnic efekt asymetrycznej reakcji na szoki. Wyniki specyfikacji przedstawiają sie następująco:
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : fGARCH(1,1)
## fGARCH Sub-Model : TGARCH
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000074 0.000019 3.9963 6.4e-05
## alpha1 0.042959 0.003421 12.5581 0.0e+00
## beta1 0.960517 0.002193 437.9118 0.0e+00
## eta11 0.928319 0.121744 7.6252 0.0e+00
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000074 0.000020 3.7358 0.000187
## alpha1 0.042959 0.002988 14.3748 0.000000
## beta1 0.960517 0.001133 847.4080 0.000000
## eta11 0.928319 0.123159 7.5376 0.000000
##
## LogLikelihood : 6981.301
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.3807
## Bayes -6.3703
## Shibata -6.3807
## Hannan-Quinn -6.3769
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.679 0.10165
## Lag[2*(p+q)+(p+q)-1][2] 5.806 0.02478
## Lag[4*(p+q)+(p+q)-1][5] 7.856 0.03200
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.3312 0.5650
## Lag[2*(p+q)+(p+q)-1][5] 2.6986 0.4649
## Lag[4*(p+q)+(p+q)-1][9] 4.7781 0.4627
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1826 0.500 2.000 0.6692
## ARCH Lag[5] 3.4418 1.440 1.667 0.2318
## ARCH Lag[7] 3.8045 2.315 1.543 0.3756
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.2693
## Individual Statistics:
## omega 1.0558
## alpha1 0.9574
## beta1 1.0549
## eta11 0.2206
##
## 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 0.5415 0.5882
## Negative Sign Bias 0.4317 0.6660
## Positive Sign Bias 1.0823 0.2792
## Joint Effect 4.8344 0.1843
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 26.34 0.1209
## 2 30 29.86 0.4211
## 3 40 38.72 0.4827
## 4 50 52.25 0.3487
##
##
## Elapsed time : 0.9781981Wszystkie parametry są istotne, również Eta11, zatem T-Garch również wskazuje,że występuje efekt dźwigni w badanej zmiennej. Podobnie jakw EGARCH nie wsytępuje autokorelacja standaryzowanych kwadratów reszt oraz efekt ARCH, jednak juz dla drugiego opóźnienia istotna jest korelacja drugiego rzędu standaryzowanych reszt. Przprowadzoną modyfikacją jest zmiana zakładanego rozkładu zmiennej objaśnianej na t-studenta.:
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : fGARCH(1,1)
## fGARCH Sub-Model : TGARCH
## Mean Model : ARFIMA(0,0,0)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000088 0.000022 3.9199 8.9e-05
## alpha1 0.044731 0.004021 11.1251 0.0e+00
## beta1 0.957682 0.002623 365.0429 0.0e+00
## eta11 0.847113 0.131896 6.4226 0.0e+00
## shape 10.600068 2.169103 4.8868 1.0e-06
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.000088 0.000022 4.0894 4.3e-05
## alpha1 0.044731 0.003102 14.4196 0.0e+00
## beta1 0.957682 0.001062 901.9899 0.0e+00
## eta11 0.847113 0.123180 6.8770 0.0e+00
## shape 10.600068 2.052390 5.1647 0.0e+00
##
## LogLikelihood : 6997.551
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.3947
## Bayes -6.3816
## Shibata -6.3947
## Hannan-Quinn -6.3899
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.797 0.09444
## Lag[2*(p+q)+(p+q)-1][2] 5.974 0.02236
## Lag[4*(p+q)+(p+q)-1][5] 8.062 0.02845
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1844 0.6676
## Lag[2*(p+q)+(p+q)-1][5] 2.5384 0.4978
## Lag[4*(p+q)+(p+q)-1][9] 4.5788 0.4938
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.222 0.500 2.000 0.6376
## ARCH Lag[5] 3.321 1.440 1.667 0.2465
## ARCH Lag[7] 3.683 2.315 1.543 0.3946
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.4511
## Individual Statistics:
## omega 0.9334
## alpha1 0.9300
## beta1 0.9712
## eta11 0.2573
## shape 0.0309
##
## 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 0.5962 0.5511
## Negative Sign Bias 0.5260 0.5989
## Positive Sign Bias 0.9925 0.3211
## Joint Effect 5.0397 0.1689
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 20.49 0.3657
## 2 30 26.87 0.5788
## 3 40 33.56 0.7160
## 4 50 40.23 0.8097
##
##
## Elapsed time : 0.9067209W tym przypadku również wskazane jest stosowanie modelu T-Garch zakładający pochodzenie rozkładu stopy zwrotu z rozkładu t-studenta. Wszytskie parametry sa istotne, kryteria niższe niż w przypadku poprzedniego modelu.Omawiany model jest ostatecznym modelem klasy T-GARCH. Ostatnim elementem analizy modelu T-Garch jest prześledzenie oszacowania warunkowej wariancji:
Podobnie jak w przypadku EGARCH(1,1) większa warunkowa wariancja wystepuje, gdy na rynku wystepuje niewpewność.
Podsumowując, przekrojowa analiza modelami klasy GARCH umożlwiła stworzenie modeli, w które są odporn na leptokurtyczność rozkładu, grupowanie wariancji oraz efek dźwigni. Ostatecznie najlepszym modelem okazał się być EGARCH(1,1) zakładający zmienna objaśnianą pochodzącą z rozkładu t-student, ponieważ ma minimalnie mniejsze kryteria informacyjne niż T-GARCH z takim samym załozeniem.
Porównanie warunkowej wariancji i wartości narażonej na ryzyko VaR w okresie in-sample
Po oszacowaniu parametrów modeli klasy GARCH, EGARCH i TGARCH zostaną zaprezentowane wykresy warunkowej wariancji odpowiednich modeli.
Na powyżej zamieszczonych wykresach możemy zaobserwować, że warunkowa wariancja jest bardzo podobna dla wszystkich 3 modeli. Dramatyczny wzrost warunkowej wariancji występuje w okolicach 2009 roku, czyli w momencie globalnego kryzysu finansowego. Jedynie model GARCH odróżnia się na tle pozostałych modeli, gdyż jego wariancja warunkowa jest stabilniejsza na przestrzeni całego analizowanego okresu i obserwowany przebieg funkcji jest bardziej płaski, co szczególnie dobrze widać w okresie od 2012 do 2015 roku.
Przejdzmy zatem do porównania wartości narażonej na ryzyko (ang. Value at Risk - VaR) w okresie in-sample
W celu obliczenia wartości narażonej na ryzyko konieczne jest oszacowanie empirycznego 1% kwantyla standaryzowanego rozkładu stopy zwrotów oraz warunkowej wariancji na podstawie modelu.
# szacowanie modelu 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 = "std")
insample.garchs11 <- ugarchfit(spec = spec,
data = na.omit(insample$r))
#standaryzacja zwrotów
insample$rstd <- (insample$r - mean(insample$r, na.rm=T)) /
sd(insample$r ,na.rm = T)
# kwantyl empiryczny
q01 <- quantile(insample$rstd, 0.01, na.rm = T)
q01## 1%
## -2.993809Po obliczeniu empirycznego 1% kwantyla, wartość narażoną na ryzyko obliczona została jako iloczyn 1% kwantyla oraz oszacowania warunkowego odchylenia standardowego na następną sesję.
# obliczanie wartości narażonej na ryzyko (VaR)
insample$VaR <- q01 * insample.garchs11@fit$sigma
head(insample$VaR)## [1] -0.03446014 -0.03464309 -0.03390682 -0.03306366 -0.03213962 -0.03125087tail(insample$VaR)## [1] -0.02479537 -0.02435567 -0.02480102 -0.02512847 -0.02452332 -0.02438595mean_VaR<-mean(insample$VaR)
mean_VaR## [1] -0.03201242Średnia wartość narażona na ryzyko wynosi -0.0320209. Jednak punktowa analiza tej wartości nie wiele wnosi do analizy. O wiele lepsze właściwości poznawcze ma wykres oszacowań wartości narażonej na ryzyko oraz ciągłej stopy zwrotu z portfela.
# wyres zwrotóW vs. VaR
plot(insample$Date, insample$r, col = "red", lwd = 1, type = 'l',
ylim = c(-0.1, 0.1),
main = "Wartość narażona na ryzyko i ciągła stopa zwrotu z aktywów")
abline(h = 0, lty = 2)
lines(insample$Date, insample$VaR, type = 'l', col = "green")
par(xpd=TRUE) #extend the area of plotting
par(new=TRUE) #to add new graph "layers"
plot(0:1, 0:1, type="n", xlab="",ylab="", axes=FALSE) #redo the x/y limits
#first legend
legend(0.9, 1.2, leg=c("r","VaR"), fill=c("red","green"))
Krzywa wartości dopasowanych Value at Risk poprawnie zmienia trend w w warunkach zwiększonej niepewności. Jednak w około 1% przypadków powinna być niedoszacowana. Aby sprawdzić stopień restrykcyjnoci modelu obliczono częstość niedoszacowania wartości narażonej na ryzyko:
freq_VaR<-sum(insample$r < insample$VaR) / length(insample$VaR)
freq_VaR## [1] 0.005944216Value at Risk jest wyższy niż stopa zwrotu w jedynie 0.6% przypadków. Obserwacja ta wskazuje, że model jest zbyt restrykcyjny. Kolejnym etapem analizy jest przeprowadzanie bliźniaczej analizy dla dwóch wybranych rozszerzeń modelu klasy GARCH.
W przypadku EGARCH, modelu najlepiej dospasowanego do danych, również oszacowano wartość narażoną na ryzyko i przedstawiono wykres:
## [1] -0.0317127Value at Risk z modelu EGARCH średnio wynosi -0.03171779. W porównaniu do oszacowań z modelu GARCH zawuażyć można, że omawiany model jest bardziej restrykcyjny w warunkach zwiększonej niepewności oraz mniej w stabilnych interwałach. Wskazują na to wyższe co do wartości bezwzględnej oszacowani wartości narażonej na ryzyko w początkach kryzysu a szczególnie mniejsze oszacowania w momentach stabliniejszej koniunktury (widoczne dla małych wahań r, że zielona linia jest poniżej niebieskiej). Wskazan cecha jest wynikiem uwzględniania w modelu efektu dźwigni, czyli asymetrycznej rekacji na szoki.Warto sprawdzić, czy model jest mniej restrykycjny:
freq_VaR2<-sum(insample$r < insample$VaR2) / length(insample$VaR2)
freq_VaR2## [1] 0.005029721rbind(freq_VaR,freq_VaR2)## [,1]
## freq_VaR 0.005944216
## freq_VaR2 0.005029721Okazuje się, że model EGARCH niedoszacowuje Value at Risk tylko w 0.503% przypadków. Jest bardziej restrykcyjny niż GARCH. Kolejnym etapem analizy oszacowań Value at Risk dla klasy modeli GARCH jest interprtacja modelu TGARCH:
## [,1]
## mean_VaR -0.03201242
## mean_VaR2 -0.03171270
## mean_VaR3 -0.03167261Średnie oszacowania Value at Risk jest najniższe w przypadku omawianego modelu i wynosi -0.03167232. Z wykresu wynika że model TGARCH generuje prognozy bardzo zbliżone do EGARCH. Podbnie poprawa w stosunku do GARCH wynika z uwzględnienia efektu dźwigni. Warto jeszcze sprawdzić stopień restrykcyjności modelu:
freq_VaR3<-sum(insample$r < insample$VaR3) / length(insample$VaR3)
freq_VaR3## [1] 0.004572474rbind(freq_VaR,freq_VaR2,freq_VaR3)## [,1]
## freq_VaR 0.005944216
## freq_VaR2 0.005029721
## freq_VaR3 0.004572474Co ciekawe, model TGARCH jest bardziej restrykcyjny jak EGARCH. Niższy odsetek warunkowej wariancji został niedoszacowany przy niższej średniej prognozie. Ostateczna decyzja, który model ma lepsze własności prognostyczne jest trudna do podjęcia, ponieważ modele TGARCH i EGARCH generują podobne wyniki.
Warunkowa wariancja i VaR w okresie out-of-sample
Następnym krokiem przeprowadzonym w raporcie jest wygenerowanie jednookresowej prognozy warunkowej wariancji dla modeli klasy GARCH, EGARCH, TGARCH na okres out-of-sample. W tym celu użyto odpowiednij pętli która generuje prognozę na jeden okres w przód dla każdej iteracji. Kod pętli prezentuje się następująco:
for (k in start:finish) {
tmp.data <- portfel[portfel$obs <= (k - 1), ]
tmp.data <- tmp.data[as.Date("2007-06-04") <= tmp.data$Date, ]
spec <- ugarchspec(variance.model = list(model = "eGARCH",
garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0),
include.mean = F),
distribution.model = "norm")
tmp.egarch11 <- ugarchfit(spec = spec, data = na.omit(tmp.data$r))
sigma.forecast <- ugarchforecast(tmp.egarch11, n.ahead = 1)
sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
warvar[k - start + 1] <- sigma.forecast2
}Zaczniemy od prognozy dla modelu GARCH(1,1) bez stałej z założeniem o rozkładzie t-studenta.
Wykres jednokresowych prognoz warunkowej wariancji prezentuje się następująco:
Na wykresie można zaobserwować, że prognoza warunkowej wariancji rośnie w momencie gdy obserwowane jest zjawisko grupowania się wariancji zwrotów w out-of-sample.
Teraz model EGARCH(1,1) bez stałej z założeniem o rozkładzie t-studenta.
Wykres jednokresowych prognoz warunkowej wariancji prezentuje się następująco:
Tutaj także prognoza warunkowej wariancji rośnie w momencie występowania zjawiska grupowania wariancji zwrotów portfela.
Przejdźmy do modelu TGARCH(1,1) bez stałej z założeniem o rozkładzie t-studenta
Wykres jednokresowych prognoz warunkowej wariancji prezentuje się następująco:
Wnioski są identyczne - warunkowa wariancja rośnie w momencie grupowania wariancji zwrotów z portfela. Porównajmy wykresy prognoz warunkowej wariancji dla wszystkich modeli:
Tutaj także wykresy zbliżone do siebie, jednak prognoza warunkowej wariancji dla modelu GARCH różni się od pozostałych. Dobrze to widać w okresie od listopada do stycznia, w którym to warunkowa wariancja modelu GARCH osiągą największy wzrost na tle modeli EGARCH i TGARCH. W tym wypadku prognoza dla modelu GARCH(1,1) cechuje się większą amplitudą wachań wariancji w analizowanym okresie. Wykresy EGARCH i TGARCH są bardziej płaskie szczególnie w okresie od stycznia do maja. Efekt jest odwrotny do warunkowej wariancji w okresie in-sample gdzie model GARCH cechował się stabilniejszą warunkową wariancją na tle pozostałych modeli.
Następnie zostaną oszaowane jednookresowe prognozy wartości narażonej na ryzyko (VaR) w okresie out-of-sample.
Zaczniemy od GARCH(1,1)
Trend VaR reaguje wzrostem w warunkach zwiększonej niepewności, gdy obserwowany jest wzrost wariancji zwrotów portfela. Zwroty przekraczają oszacowaną wartość VaR tylko w jednym przypadku. Aby zbadać restrykcyjność modelu obliczona została częstość niedoszacowania wartości narażonej na ryzyko:
## [1] 0.004081633Zwroty przekraczają wartość narażoną na ryzyko tylko w około 0,4% przypadków. Przejdźmy teraz do modelu EGARCH(1,1) bez stałej zakładającego rozkład t-studenta
EGARCH
Trend wartości narażonej na ryzyko dla modelu EGARCH bardzo podobny do modelu GARCH jednak występują minimalne różnice
W ilu przypadkach zwroty przekraczają VaR?
## [1] 0.004081633W 0,4% przypadków. Wartość tej statystyki jest identyczna jak w modelu GARCH.
Dla TGARCH:
W ilu przypadkach przekracza?
## [1] 0.004081633Tutaj także w 0,4%.
Wykresy wartości narażonej na ryzyko są bardzo podobne w modelach EGARCH i TGARCH jedynie VaR dla modelu GARCH delikatnie się różni. Zwroty z portfela przekraczają Value at Risk w takiej samej ilości przypadków dla wszystkich modeli. Na podstawie tego kreyterium nie można zatem wybrać najlepszego modelu.
Analiza wrażliwości przy zmianie długości okna do estymacji parametrów
Warto jednak sprawdzić, zachowanie prognoz w ostatnim roku,w tym celu przeprowadzona została podbna anliza, tylko na okresie od 2015-06-01 do 2016-06-01. Interpretacja graficzna pozwoli na zdecydowanie, który model wydaje się być odpowiednim do prognozowania wartości narazonej na ryzyko dla któszego horyznotu.
Próba zostaje ograniczona do ostatniego roku w celu pominięcia największych wahań stóp zwrotu
## [,1]
## mean2_VaR -0.02778383
## mean2_VaR2 -0.03061045
## mean2_VaR3 -0.02754326## [,1]
## freq2_VaR 0.012244898
## freq2_VaR2 0.008163265
## freq2_VaR3 0.012244898Na podstawie wykresu można zauważyć, że model TGARCH częściej prognozuje znacznie niższą wartość narażoną na ryzyko przy 50% większej częstości niedoszacowania co EGARCH. Ponadto, na badanej próbie nawet GARCH charakteryzuje się niższym niedoszacowaniem VaR. Wybór najlepszego modelu jest w tym przypadku trudny. Model GARCH cechował się najgorszymi własnościami na całej analizowanej próbie oraz nie uwzględnia efektu dźwigni, więc wybór należy ograniczyć do TGARCH i EGARCH. Pierwszy z modeli ma najwyższą częstość niedoszacowania VaR i nie mieści się ona w przyjmowanym poziomie 1%. z kolei EGARCH ma niższą częstość ale prognozuje wyższą VaR. na portfelu 100000 zł różnica średnich VaR wynosi 330.4 zł. Oszacowanie wartości narażonej na ryzyko EGARCH jest o około 11% niższe niż TGARCH, z kolei prawdopodbieństwo utraty 330.4 zł jest 50% wyższe. Warto sprawdzić to na tym konkretnym przykładzie,czy model EGARCH i TGARCH niedoszacowały VaR dla tej samej obserwacji
##
## m_Var3 nm_Var3
## m_Var2 2 0
## nm_Var2 1 242W badananym rocznym okresie model TGARCH niedoszacował VaR trzy razy, z kolei EGARCH raz, z tym, że dwie pomyłki były wspólne. W tym przypadku, ze względu na awersje do ryzyka, inwestorzy wybrali model EGARCH za odpowiedni. Kolejnym etapem analizy jest wygenerowanie prognoz dla próby out-of-sample. Wygenerowane zostaną wykresy prognoz jednookresowych out-of-sample dla VaR dla modeli GARCH(1,1), EGARCH(1,1) i TGARCH(1,1) dla 2 różnych okresów in-sample: 1 letni i 9 letni.
W wspomnianych modelach zwroty przekraczają VaR:
GARCH dla 1 roku:
## [1] 0.008163265GARCH dla 9 lat:
## [1] 0.004081633EGARCH dla 1 roku:
## [1] 0.008163265EGARCH dla 9 lat:
## [1] 0.004081633TGARCH dla 1 roku:
## [1] 0.0122449TGARCH dla 9 lat:
## [1] 0.004081633Jak można się było spodziewać we wszystkich trzech modelach skrócenie okna estymacji wpływa negatywnie na przekroczenia wartości narażonej na ryzyko. Funkcja prezentująca wartości VaR dla estymacji 1 rocznej znajduje się powyżej VaR 9 letniej co oznacza, że częściej będzie dochodzić do przekroczeń. Wyżej przedstawione wartości potwierdzają wnioskowanie gdyż w przypadku estymacji 9 letniej zwroty z portfela rzadziej przekraczają VaR , są to więc modele bardziej restrykcyjne. W przypadku okna estymacji wynoszącego 1 rok wszystkie trzy modele są równie restrykcyjne. Wybór najlepszego jest zatem trudny.
Podsumowanie
W pracy została przeprowodzona analiza porównawcza oszacowań warunkowej wariancji i wartości narażonej na ryzyko tj. Value at Risk dla 3 modeli klasy GARCH. Podczas poszukiwania modeli najlepiej dopasowanych do danych i ich parametrów wytypowane zostały trzy modele: GARCH(1,1) i dwa jego rozszerzenia EGARCH(1,1) i TGARCH(1,1). Wszystkie wspomniane modele zakładały, że zwroty z portfela pochodzą z rozkładu t-studenta, gdyż założenie to znacznie poprawiało wartości kryteriów informacyjnych. Modele EGARCH i TGARCH zawsze prezentowały podobne wyniki zarówno w oszacowaniach i prognozach warunkowej wariancji jak i w przypadku VaR. Model GARCH prezentował drobne lecz nieistotne różnice. Z tego też powodu decyzja o wyborze najlepszego modelu była trudna do podjęcia. Po skróceniu długości okna do estymacji do jednego roku najlepszy okazał się model EGARCH gdyż w jego przypadku zwroty z portfela w najmniejszym stopniu przekraczały wartość VaR. Skrócenie okna estymacji parametrów zwiększa liczbę przekroczeń VaR.