Streszczenie

Praca ma na celu przeprowadzenie analizy ryzyka równoważonego portfela aktywów. Do przeprowadzenia analizy porównawczej wykorzystano model z klasy GARCH i dwa rozszerzenia tego modelu: EGARCH oraz TGARCH. Dokonano wyboru najlepszych modeli oszacowań warunkowej wariancji oraz wartości narażonej na ryzyko w okresach in-sample i out-of-sample. Wszystkie wybory zostały potwierdzone odpowiednimi analizami i testami.

Niezbędne pakiety

library(xts)
library(tseries)
library(car)
library(FinTS) 
library(fGarch) 
library(rugarch) 
library(MSBVAR)
library(lmtest) 
library(fBasics) 
library(urca) 
library(vars)
library(dplyr)

Wybór instrumentów

Do przeprowadzenia analizy porównawczej ryzyka rozumianego jako oszacowanie funkcji warunkowej wariancji w modelach klasy GARCH został zbudowany portfel składający się z pięciu instrumentów finansowych. Dane zostały zaczerpnięte z portalu notowań indeksów stooq.com. Wykorzystano szereg zawierający notowania trzech spółek, w tym dwóch polskich:

-mBank S.A.

-PKO BP S.A.

-J.P. Morgan Chase & Co

oraz dwa notowania kursów walutowych:

-kurs dolara amerykańskiego do euro

-kurs dolara amerykańskiego do funta brytyjskiego

PKO BP S.A. oraz mBank S.A. to spółki notowane na polskiej giełdzie. Obie należą do sektora finansowego i są częścią indeksu WIG20. Spółka J.P. Morgan Chase & Co to również spółka z sektora finansowego, ale jest ona notowana na Giełdzie Papierów Wartościowych w Nowym Jorku i jest on częścią indeksu DJIA. Wszystkie wykorzystane w portfelu spółki charakteryzują się wysoką płynnością.

Wczytanie danych

url1 <- "https://stooq.pl/q/d/l/?s=mbk&i=d" # mbank
url2 <- "https://stooq.pl/q/d/l/?s=usdgbp&i=d" #USDGDP
url3 <- "https://stooq.pl/q/d/l/?s=usdeur&i=d" # USDEUR
url4 <- "https://stooq.pl/q/d/l/?s=pko&i=d"    # pko
url5 <- "https://stooq.pl/q/d/l/?s=jpm.us&i=d"   # jp

Przygotowanie danych do analizy

Jako dolne ograniczenie szeregu przyjęto datę 2 lutego 2012 roku. Górnym ograniczeniem jest natomiast dzień 21 lipca 2017 roku. Górne ograniczenie wynika z daty rozpoczęcia analizy badanego portfela, natomiast dolne zostało wybrane w wyniku analizy szeregów pojedynczych instrumentów. Wybór wcześniejszej daty zaburzałby wyniki oszacowań modeli, przede wszystkim ze wzgledu na zaburzenia wynikające z kryzysu w 2008 roku i jego echa w latach późniejszych.

Dane dla poszczególnych instrumentów wczytano i przygotowano do analizy w ten sam sposób. Dla przykładu kod dla spółki mBank S.A.

#MBANK
MBANK <- read.csv(url1, 
                  header = TRUE,
                  sep = ",",
                  dec = ".",
                  stringsAsFactors = F)
colnames(MBANK)[5]=cbind("MBANK")
MBANK <-subset(MBANK,select=c(Data,MBANK))
MBANK <- MBANK[MBANK$Data >= "2012-01-02",]
MBANK <- MBANK[MBANK$Data <= "2017-07-21",]
MBANK$Data <- as.Date(MBANK$Data)

Budowa portfela

W pierwszym kroku wyliczono logarytmiczne ciągłe stopy zwrotu z zamknięcia dla poszczególnych instrumentów finansowych. Poniżej procedura dla mBanku oraz wykresy dla wszystkich instrumentów:

MBANK$r5 <- diff.xts(log(MBANK$MBANK))
plot(MBANK$Data, MBANK$r5, type = "l", col = "black", main = "Wykres ciągłych stóp zwrotu dla spółki mBank")

Na wykresach można zaobserwować, że stopy zwrotu wykazują duże wahania. Zdarzają się obserwacje znacznie odbiegające od przeciętniej wartości. Widoczna jest również asymetria spadków i wzrostów, co do wartości bezwzględnej spadki są większe niż wzrosty.
Na wykresach zwrotów dla kursów walutowych warto zwrócić uwagę na znaczne wzrosty w pierwszej połowie 2016 roku. Miały one związek z sytuacją polityczną w Stanach Zjednoczonych, a konkretnie z wyborami prezydenckim. Pozostałe instrumenty również znacznie zareagowały na wspomniane wybory prezydenckie. Wynika to przede wszystkim z faktu, ze wszystkie spółki pochodzą z sektora finansowego, dlatego są podatne na wahania kursowe głównych walut.

W kolejnym kroku logarytmiczne ciągłe stopy zwrotu zostały połączone do jednej tablicy danych:

s1 <- merge(PKO, USDGBP, by.x="Data", by.y="Data", all.x=TRUE)
s2 <- merge(s1, JP, by.x="Data", by.y="Data", all.x=TRUE)
s3 <- merge(s2, USDEUR, by.x="Data", by.y="Data", all.x=TRUE)
EWP <- merge(s3, MBANK, by.x="Data", by.y="Data", all.x=TRUE)

Następnie zbudowano portfel równoważony aktywów. Dla portfela równoważonego każdy instrument wchodzi do portfela z tą samą wagą, dlatego w celu zbudowania tego portfela uśredniono ciągłe stopy zwrotu dla danego dnia.

EWP <- mutate(EWP, average_ror=((r1+r2+r3+r4+r5)/5))
EWP <- mutate(EWP, marking=((USDGBP+MBANK+JP+PKO+USDEUR)/5))
EWP <- subset(EWP, select=c(Data,average_ror,marking))
data <- na.omit(EWP)
EWP <- data

Dodano również numer dla każdej obserwacji:

EWP$obs<- 1:length(EWP$average_ror)

Aby nie zaburzyć prognozy pomijano obserwacje charakteryzujące się brakiem danych. Wykres stóp zwrotów dla portfela zrównoważonego prezentuje się następująco:

plot(EWP$Data, EWP$average_ror, type = "l", 
     main = "Wykres stóp zwrotów portfela", col = "red",
     xlab ="Data", ylab="zwroty")
abline(h = 0, col="black")

Statystyki opisowe oraz wstępne testy

Na początek przedstawiono statystyki opisowe dla badanego szeregu.

##             X..EWP.average_ror
## nobs               1347.000000
## NAs                   0.000000
## Minimum              -0.032929
## Maximum               0.037960
## 1. Quartile          -0.004102
## 3. Quartile           0.004648
## Mean                  0.000353
## Median                0.000373
## Sum                   0.475905
## SE Mean               0.000204
## LCL Mean             -0.000047
## UCL Mean              0.000753
## Variance              0.000056
## Stdev                 0.007487
## Skewness              0.112429
## Kurtosis              1.663234

Dane zawierają wartości z przedziału [-0.032929 ; 0.037960]. Współczynnik skośności jest dodatni, co sugeruje prawostronną asymetrię rozkłady. Kurtoza będąca miarą spłaszczeni rozkładu wynosi 1,663234 . Dodatnia kurtoza sugeruje, że rozkład jest leptokurtyczny, czyli wartości cechy są bardziej skoncentrowane niż w przypadku rozkładu normalnego.

W kolejnym kroku narysowano histogram oraz wykres Q-Q dla stóp zwrotu z portfela.

hist(EWP$average_ror, prob = T, breaks = 100, main="Histogram zwrotów")
curve(dnorm(x,
            mean = mean(EWP$average_ror, na.rm = T),
            sd = sd(EWP$average_ror, na.rm = T)),
      col = "darkblue", lwd = 2, add = TRUE)
abline(v=0, col = "red", lwd = 2)

qqnorm(EWP$average_ror, col = "darkblue")
qqline(EWP$average_ror, distribution = qnorm, col = "red", lwd = 1)
abline(h=0, col = "darkgreen", lwd = 2)

Smukłość histogramu wskazuje na rozkład leptokurtyczny. Można zauważyć również grube ogony rozkładu. Wykres Q-Q dla stóp zwrotu nie zawiera się w prostej wyznaczonej dla rozkład normalnego. Odbiega od niej w wielu miejscach i często wykracza poza przedział ufności.

Analiza powyższych wykresów pozwala stwierdzić, że badany rozkład stóp zwrotu dla zrównoważonego portfela nie ma rozkładu normalnego. W dalszej części zostanie przeprowadzona dokładniejsza analiza za pomocą testów.

Test Jarque-Bera

jarque.bera.test(na.omit(EWP$average_ror))
## 
##  Jarque Bera Test
## 
## data:  na.omit(EWP$average_ror)
## X-squared = 159.4, df = 2, p-value < 2.2e-16

Przeprowadzony test silnie odrzuca \(H_0\) mówiącą o normalności rozkładu zwrotów. Tym samym zwroty dla badanego portfela nie pochodzą z rozkładu normalnego.

Test Dublina-Watsona

durbinWatsonTest(lm(EWP$average_ror ~ 1),
                 max.lag = 5)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.05255135      1.894042   0.046
##    2     -0.09045395      2.179801   0.002
##    3     -0.04273685      2.083366   0.112
##    4     -0.03950977      2.075899   0.148
##    5     -0.02830753      2.053098   0.258
##  Alternative hypothesis: rho[lag] != 0

Przy wykorzystaniu testu Durbina - Watsona stwierdzono występowanie istotnego opóźnienia 2, opóźnienie 1 jest na granicy istotności. Przyjęto poziom istotności 5%.

Ljung-Box test

Box.test(EWP$average_ror, type = "Ljung-Box", lag = 24)
## 
##  Box-Ljung test
## 
## data:  EWP$average_ror
## X-squared = 54.316, df = 24, p-value = 0.0003868

Przeprowadzony test odrzuca \(H_0\) o braku autokorelacji.

Test efektów ARCH

ArchTest(EWP$average_ror, lags = 5)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  EWP$average_ror
## Chi-squared = 26.03, df = 5, p-value = 8.804e-05

P-value jest mniejsze od wartości krytycznej 0,05, więc test odrzuca \(H_0\) o braku występowania efektów ARCH.

ACF zwrotów z portfolio

Na wykresie ACF stóp zwrotów z portfela zaobserwowano istotne opóźnienia: 1, 2, 3, 8. Dla wykresu PACF tóp zwrotów z portfela zaobserwowano istotne opóźnienia: 1, 2, 11.

acf(EWP$average_ror,  lag.max = 36,
    lwd = 7,
    col = "darkblue", na.action = na.pass,
    main="Wykres ACF zwrotów z portfela")

pacf(EWP$average_ror, lag.max=36,
     lwd = 7,
     col = "darkblue", na.action = na.pass,
 main="Wykres PACF zwrotów z portfela")

Wykresy dla kwadratów zwrotów z portfela.

Na wykresach ACF i PACF stóp zwrotów dla portfela widać, że wiele opóźnień zyskało na istotności w stosunku do zwykłych stóp zwrotów.

Wybór modeli z rodziny GARCH

Rozszerzeniem modelu ARCH (q) jest model uogólnionej warunkowej autoregresyjnej heteroskedastyczności GARCH (p,q). Model ten ma następującą postać:
\(y_t = x_t\beta+\epsilon_t\)
\(\epsilon_t=u_t+\sigma_t\)
\(\sigma^2_t=\alpha_1\sigma^2_{t-1}+\alpha_1\sigma^2_{t-2}+\dots+\theta_0+\theta_1\epsilon^2_{t-1}+\dots+\theta_q\epsilon^2_{t-q}\)

gdzie \(\sigma^2_t=E(\epsilon_t^2|\epsilon_{t-1}\dots\epsilon_{t-q})\) oraz \(\alpha_i\geq0\), \(\theta_i\geq0\).

Model GARCH

Analizie poddano kilka modeli GARCH, o różnych wartościach p i q. Wartości kryteriów informacyjnych dla najważniejszych modeli podano w tabeli poniżej.

Ocena modelu GARCH według kryteriów informacyjnych

P Q Akaike Bayes Shibata Hannan-Quinn
1 1 -6,9904 -6,9750 -6,9904 -6,9846
1 2 -6,9895 -6,9701 -6,9895 -6,9822
2 1 -6,9889 -6,9696 -6,9889 -6,9817
2 2 -6,9880 -6,9648 -6,9880 -6,9793

Powyższa analiza sugeruje, że najlepiej dopasowanym modelem pod względem kryteriów informacyjnych jest model GARCH(1,1). Charakteryzował się on również najlepszymi wynikami testów, dlatego zostanie on wykorzystany w dalszej analizie.

Poniżej przedstawiono oszacowania wybranego modelu.

spec <- ugarchspec(variance.model = list(model = "sGARCH",
                                      garchOrder = c(1, 1)),
                                      mean.model = list(armaOrder = c(0, 0),
                                      include.mean = T),
                                      distribution.model = "norm"
                                      )

k.garch11a <- ugarchfit(spec = spec,
                         data = na.omit(EWP$average_ror))

ewp.xts <- xts(na.omit(EWP$average_ror),
               order.by  = EWP$Data[1:length(EWP$Data)])

k.garch11a <- ugarchfit(spec = spec,
                         data = ewp.xts)
k.garch11a
## 
## *---------------------------------*
## *          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|)
## mu      0.000396    0.000191  2.06867 0.038577
## omega   0.000001    0.000001  0.82984 0.406632
## alpha1  0.038934    0.016237  2.39795 0.016487
## beta1   0.942352    0.018998 49.60215 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000396    0.000175 2.256698 0.024027
## omega   0.000001    0.000016 0.067927 0.945844
## alpha1  0.038934    0.166370 0.234023 0.814967
## beta1   0.942352    0.206306 4.567732 0.000005
## 
## LogLikelihood : 4712.041 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.9904
## Bayes        -6.9750
## Shibata      -6.9904
## Hannan-Quinn -6.9846
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic  p-value
## Lag[1]                       3.91 0.048006
## Lag[2*(p+q)+(p+q)-1][2]      8.12 0.006013
## Lag[4*(p+q)+(p+q)-1][5]     12.86 0.001649
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2401  0.6241
## Lag[2*(p+q)+(p+q)-1][5]    0.4678  0.9625
## Lag[4*(p+q)+(p+q)-1][9]    1.5228  0.9537
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.007396 0.500 2.000  0.9315
## ARCH Lag[5]  0.642436 1.440 1.667  0.8409
## ARCH Lag[7]  1.547473 2.315 1.543  0.8118
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  304.6117
## Individual Statistics:               
## mu      0.09581
## omega  25.82689
## alpha1  0.23014
## beta1   0.29317
## 
## 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.1300 0.8966    
## Negative Sign Bias  0.9976 0.3187    
## Positive Sign Bias  0.3396 0.7342    
## Joint Effect        1.2483 0.7414    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     22.93       0.2403
## 2    30     31.60       0.3378
## 3    40     48.13       0.1499
## 4    50     60.76       0.1209
## 
## 
## Elapsed time : 0.205138

Suma parametrów jest mniejsza od 1, nie wszystkie oszacowania są istotne, jednak usunięcie zmiennych pogarsza jakosć modelu, dlatego zdecydowano się na ich zatrzymanie. Ljung - Box test na wystandaryzowanych resztach pozwala odrzucić hipotezę o braku autokorelacji, jednak ten sam test przeprowadzony na kwadratach reszt nie pozwala na odrzucenie hipotezy zerowej mówiącej o braku autokorelacji dla żadnego z opóźnień. Test efektów ARCH nie pozwala na odrzucenie hipotezy zerowej mówiącej o braku efektu ARCH, co oznacza, że model wychwycił ten efekt, ponieważ występował on we wstępnych testach. Test Pearsona wskazuje na to, że rozkład modelu jest zgodny z założonym.

Wariancja warunkowa i bezwarunkowa

Jak wiemy, wariancja bezwarunkowa jest określona jako:

\(\sigma^2=\frac{\theta_0}{1-\theta_1- \dots-\theta_q}\).

Oszacowano model GARCH(1,1).

k.garch11 <- garchFit(formula = ~ garch(1, 1),
                        data = na.omit(EWP$average_ror),
                        include.mean = F,
                        cond.dist = "norm",
                        trace = F)
summary(k.garch11)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = na.omit(EWP$average_ror), 
##     cond.dist = "norm", include.mean = F, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x00000000131c3780>
##  [data = na.omit(EWP$average_ror)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1  
## 1.2397e-06  4.0627e-02  9.3765e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  1.240e-06   6.771e-07    1.831 0.067104 .  
## alpha1 4.063e-02   1.063e-02    3.821 0.000133 ***
## beta1  9.377e-01   1.950e-02   48.076  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  4709.966    normalized:  3.496634 
## 
## Description:
##  Thu Aug 10 01:45:23 2017 by user: Admin 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  128.8266  0           
##  Shapiro-Wilk Test  R    W      0.9891678 1.943606e-08
##  Ljung-Box Test     R    Q(10)  20.54072  0.02453382  
##  Ljung-Box Test     R    Q(15)  26.94877  0.02915653  
##  Ljung-Box Test     R    Q(20)  30.61335  0.06050625  
##  Ljung-Box Test     R^2  Q(10)  3.64353   0.9619987   
##  Ljung-Box Test     R^2  Q(15)  5.939292  0.9807495   
##  Ljung-Box Test     R^2  Q(20)  7.623024  0.9940825   
##  LM Arch Test       R    TR^2   5.272946  0.9482354   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.988813 -6.977219 -6.988823 -6.984471

Oceny parametrów zapisane są w:

k.garch11@fit$par
##        omega       alpha1        beta1 
## 1.239738e-06 4.062723e-02 9.376513e-01

Oszacowania bezwarunkowej wariancji:

var_uncond1 <- k.garch11@fit$par[1] / (1 - k.garch11@fit$par[2]
                                         - k.garch11@fit$par[3])
var_uncond1
##        omega 
## 5.707424e-05

W kolejnym kroku dokonano oszacowań prognozy warunkowej wariancji na kolejne 200 okresów

ewp.fore200 <- predict(k.garch11, n.ahead = 200)
head(ewp.fore200)
##   meanForecast   meanError standardDeviation
## 1            0 0.006959340       0.006959340
## 2            0 0.006972814       0.006972814
## 3            0 0.006985969       0.006985969
## 4            0 0.006998815       0.006998815
## 5            0 0.007011359       0.007011359
## 6            0 0.007023609       0.007023609
plot(ewp.fore200[, 3] ^ 2, type = "l")
abline(h = var_uncond1, col = "red", lty = 2)
title(main = "Warunkowa i bezwarunkowa wariancja zwrotów portfela")

Można zaobserwować, że w długim okresie wariancja warunkowa zbiega do wariancji bezwarunkowej.

Model TGARCH

Ocena modelu TGARCH według kryteriów informacyjnych

P Q Akaike Bayes Shibata Hannan-Quinn
1 1 -6,9997 -6,9842 -6,9997 -6,9939
1 2 -6,9982 -6,9788 -6,9982 -6,9909
2 1 -6,9964 -6,9733 -6,9965 -6,9878
2 2 -6,9994 -6,9723 -6,9994 -6,9892

Wszystkie kryteria zgodnie wskazują na model TGARCH(1,1), jednak wyniki testów dla tego modelu okazły się grosze niż dla modelu zajmujących drugie miejsce, tj. model TGARCH(1,2) oraz model TGARCH(2,2). Wyniki testu Persona dla tych modeli są lepsze niż dla modelu TGARCH(1,1), ale model TGARCH(2,2) ma więcej nieistotnych parametrów, tym samym dalszej analizie poddano model TGARCH(1,2).

spec <- ugarchspec(
  variance.model = list(model = "fGARCH",
                        garchOrder = c(1, 2),
                        submodel = "TGARCH"),
                        mean.model = list(armaOrder = c(0, 0),
                        include.mean = F),
                        distribution.model = "norm")

k.tgarch12a <- ugarchfit(spec = spec, data = na.omit(EWP$average_ror))
k.tgarch12a
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : fGARCH(1,2)
## fGARCH Sub-Model : TGARCH
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## omega   0.000184    0.000081   2.2769 0.022792
## alpha1  0.052246    0.011281   4.6312 0.000004
## beta1   0.725517    0.033041  21.9581 0.000000
## beta2   0.211045    0.030560   6.9059 0.000000
## eta11   0.545727    0.193460   2.8209 0.004789
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## omega   0.000184    0.000100   1.8528 0.063911
## alpha1  0.052246    0.012628   4.1374 0.000035
## beta1   0.725517    0.006568 110.4600 0.000000
## beta2   0.211045    0.012029  17.5446 0.000000
## eta11   0.545727    0.225972   2.4150 0.015734
## 
## LogLikelihood : 4718.265 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.9982
## Bayes        -6.9788
## Shibata      -6.9982
## Hannan-Quinn -6.9909
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic  p-value
## Lag[1]                      3.941 0.047132
## Lag[2*(p+q)+(p+q)-1][2]     7.856 0.007064
## Lag[4*(p+q)+(p+q)-1][5]    12.148 0.002547
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.2188  0.6400
## Lag[2*(p+q)+(p+q)-1][8]     1.7556  0.8946
## Lag[4*(p+q)+(p+q)-1][14]    3.1819  0.9449
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]  0.008857 0.500 2.000  0.9250
## ARCH Lag[6]  1.746365 1.461 1.711  0.5492
## ARCH Lag[8]  2.466787 2.368 1.583  0.6472
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.145
## Individual Statistics:             
## omega  0.4072
## alpha1 0.3102
## beta1  0.3447
## beta2  0.3437
## eta11  0.1657
## 
## 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.2595 0.7953    
## Negative Sign Bias  0.4412 0.6591    
## Positive Sign Bias  0.5857 0.5582    
## Joint Effect        0.7635 0.8582    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     29.94      0.05255
## 2    30     45.90      0.02403
## 3    40     54.19      0.05371
## 4    50     64.25      0.07076
## 
## 
## Elapsed time : 0.7434959

Ljung - Box test dla wystandaryzowanych reszt wskazuje na odrzucenie hipotezy zerowej mówiącej o braku autokorelacji. Jednak ten sam test dla wystandaryzowanych kwadratów reszt pokazuje, że nie ma podstaw do odrzucenia hipotezy zerowej mówiącej o braku auotkorelacji. Według test efektów ARCH nie ma podstaw do odrzucenia hipotezy zerowej mówiącej o braky efektu ARCH, tym samym model wychwycił wcześniej występujący efekt ARCH. Test stabilności Nybloma wskazuje na to, że model jest stabilny przy 5% poziomie istotnści. Test Pearsona nie daje jednoznacznych wyników, jednak dla większosci przypadków nie pozwala odrzucić hipotezy mówiącej o tym, że rozkład jest zgodny z założonym.

Wariancja warunkowa i bezwarunkowa

Następnie dokonano oszacowania modelu TGARCH(1,2) oraz oszacowano bezwarunkową wariancję dla badanego modelu.

##        omega       alpha1        beta1        beta2 
## 1.604933e-06 5.684766e-02 4.931516e-01 4.220035e-01
## unconditional variance 
##           5.732474e-05

Oszacowano prognozy wariancji warunkowej na kolejne 200 okresów i narysowano wykres dla warunkowej wariancji portfela dla modelu TGARCH(1,2).

Model EGARCH

Ocena modelu EGARCH według kryteriów informacyjnych

P Q Akaike Bayes Shibata Hannan-Quinn
1 1 -7,0012 -6,9857 -7,0012 -6,9954
1 2 -6,9999 -6,9806 -6,9999 -6,9926
2 1 -7,0000 -6,9768 -7,0000 -6,9913
2 2 -7,0026 -6,9756 -7,0027 -6,9925

Kryteria informacyjne w tym wypadku wskazują na dwa modele: EGARCH(1,1) oraz EGARCH(2,2). Testy w obu przypadkach dają podobne oszacowania, ale model EGARCH(2,2) zawiera nieistotne parametry, w zwiazku z tym uznano, że lepszym modelem będzie 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")
k.egarch11a <- ugarchfit( spec = spec, data = na.omit(EWP$average_ror))
k.egarch11a
## 
## *---------------------------------*
## *          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.197566    0.002664  -74.1736    0e+00
## alpha1 -0.050568    0.011855   -4.2656    2e-05
## beta1   0.979428    0.000432 2268.7333    0e+00
## gamma1  0.082873    0.007282   11.3799    0e+00
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## omega  -0.197566    0.006310  -31.3090 0.000000
## alpha1 -0.050568    0.013962   -3.6219 0.000292
## beta1   0.979428    0.000612 1600.4942 0.000000
## gamma1  0.082873    0.013053    6.3489 0.000000
## 
## LogLikelihood : 4719.291 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -7.0012
## Bayes        -6.9857
## Shibata      -7.0012
## Hannan-Quinn -6.9954
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic  p-value
## Lag[1]                       3.85 0.049740
## Lag[2*(p+q)+(p+q)-1][2]      7.85 0.007090
## Lag[4*(p+q)+(p+q)-1][5]     12.32 0.002295
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3997  0.5272
## Lag[2*(p+q)+(p+q)-1][5]    0.7910  0.9049
## Lag[4*(p+q)+(p+q)-1][9]    2.0374  0.9005
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2477 0.500 2.000  0.6187
## ARCH Lag[5]    0.9637 1.440 1.667  0.7438
## ARCH Lag[7]    1.9827 2.315 1.543  0.7210
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.903
## Individual Statistics:             
## omega  0.3904
## alpha1 0.1673
## beta1  0.4034
## gamma1 0.1223
## 
## 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.2040 0.8384    
## Negative Sign Bias  0.5911 0.5545    
## Positive Sign Bias  0.6971 0.4859    
## Joint Effect        0.9573 0.8116    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     27.86      0.08611
## 2    30     44.20      0.03511
## 3    40     55.08      0.04545
## 4    50     59.87      0.13743
## 
## 
## Elapsed time : 0.278182

Ljung - Box test dla wystandaryzowanych reszt w każdym z przypadków wskazuje na odrzucenie hipotezy zerowej mówiącej o braku autokorelacji, jednak ten sam test dla kwadratów reszt wskazuje, że nie ma podstaw do odrzucenia hipotezy zerowej w żadnym z przypadków. Test efektów ARCH nie pozwala odrzucić hipotezy o braku wystepowania efektów ARCH, tym samym można wnioskować, że model poradził sobie z wcześniej występującymi efektami ARCH. Test Nybloma wskazuje na stabilność modelu. Zgodnie z testem Pearsona można założyć, że dla większości przypadków rozkład jest zgodny z założonym, aczkolwiek wynik jest dyskusyjny.

Wariancja warunkowa i bezwarunkowa

Następnie dokonano oszacowania modelu EGARCH(1,1) oraz oszacowano warunkową oraz bezwarunkową wariancję dla badanego modelu.

##        omega       alpha1        beta1 
## 1.239738e-06 4.062723e-02 9.376513e-01
## unconditional variance 
##           5.707424e-05

Dokonano oszacowań prognozy wariancji warunkowej na kolejne 200 okresów dla modelu EGARCH(1,1).

ewp.fore200 <- predict(k.egarch22, n.ahead = 200)
head(ewp.fore200)
##   meanForecast   meanError standardDeviation
## 1            0 0.006959340       0.006959340
## 2            0 0.006972814       0.006972814
## 3            0 0.006985969       0.006985969
## 4            0 0.006998815       0.006998815
## 5            0 0.007011359       0.007011359
## 6            0 0.007023609       0.007023609
plot(ewp.fore200[, 3] ^ 2, type = "l")
abline(h = var_uncond3, col = "red", lty = 5)
title(main = "Warunkowa i bezwarunkowa wariancja zwrotów portfela")

Oszacowanie wartości narażonej na ryzyko w okresie in-sample

W kolejnym kroku oszacowano VaR w okresie in-sample. Jest to oszacowanie pokazują możliwe wielkości straty portfela przy określonych prawdopodobieństwach. Okres in sample obejmuje nieco ponad 5 lat natomiast okres out-of sample będzie obejmował około pół roku.

Model GARCH (1,1)

Standaryzacja zwrotów

W pierwszym kroku obcięto szereg stóp zwrotów i dokonano standaryzacji dla portfela oraz wyliczono 1% kwantyl z próby.

EWP_o <- EWP[EWP$Data <= "2017-02-02", ]
EWP_o$rstd <- (EWP_o$average_ror - mean(EWP_o$average_ror, na.rm=T)) / sd(EWP_o$average_ror ,na.rm = T)
q01 <- quantile(EWP_o$rstd, 0.01, na.rm = T)
q01
##        1% 
## -2.434452

1% kwantyl z próby wynosi -2,434452.

Wyznaczono wartość narażoną na ryzyko będący uloczynem wektora oszacowań odchylenia standardowego i 1% kwantyla empirycznego rozkładu. Procedura została powtórzona dla każdego z badanych modeli.

EWP$VaR <- q01 * k.garch11@sigma.t

Wykres zwrotów i wartości narażonej na ryzyko dla portfela

Poniżej zaprezentowano wykresy zwrotów oznaczony czerwoną linią oraz wykres wartości narażonej na ryzyko dla badanego portfela.

plot(EWP$Data, EWP$average_ror, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.1, 0.1))
abline(h = 0, lty = 2)
lines(EWP$Data, -EWP$VaR, type = 'l', col = "green") 
lines(EWP$Data, EWP$VaR, type = 'l', col = "green")

Z wizualnej analizy można wnioskować, że występują przekroczenia. Aby określic dokładny procent przekroczeń użyto równania.

sum(EWP$average_ror < EWP$VaR) / length(EWP$VaR)
## [1] 0.008908686

W przypadku modelu GARCH(1,1) straty przekroczyły zakładany poziom wartości narażonej na ryzyko w 0,89% przypadków.

Model TGARCH (1,2)

Wykres zwrotów i wartości narażonej na ryzyko dla portfela

## [1] 0.009651076

W przypadku modelu TGARCH(1,2) straty przekroczyły zakładany poziom wartości narażonej na ryzyko w 0,965% przypadków.

Model EGARCH (1,1)

Wykres zwrotów i wartości narażonej na ryzyko dla portfela

## [1] 0.008908686

W przypadku modelu EGARCH (1,1) straty przekroczyły zakładany poziom wartości narażonej na ryzyko w 0,89% przypadków.

Analizując częstość przekroczeń strat przy założonym poziomie VaR można zauważyć, że najgorszym modelem jest model TGARCH(1,2), w którym straty przekroczyły założony VaR tylko w 0,965%. W pozostałych dwóch przypadkach przekroczenia stanowią ten sam procent.

Szacowanie wartości narażonej na ryzyko w okresie out-of-sample

W ostatnim kroku dokonano oszacowania wartości narażonej na ryzyko w okresie out-of-sample dla posczególnych modeli. Szereg obcięto do dnia 3.02.2017.

Model GARCH (1,1)

EWP$obs<-1:length(EWP$average_ror)
EWP$obs<-1:length(EWP$Data)
start  <- EWP$obs[EWP$Data == as.Date("2017-02-03")]
finish <- EWP$obs[EWP$Data == as.Date("2017-07-21")]
EWP2 <-EWP[start:finish, ]
VaR2 <- rep(NA, times = finish - start + 1)
for (k in start:finish) {
  tmp.data <- EWP[EWP$obs <= k-1, ]
  tmp.data <- tmp.data[as.Date("2012-01-02") <= tmp.data$Data, ]
  tmp.data$rstd <- (tmp.data$average_ror - mean(tmp.data$average_ror, na.rm = T)) /
    sd(tmp.data$average_ror, 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 = T),
                                      distribution.model = "norm")
  tmp.egarch22 <- ugarchfit(spec = spec, data = na.omit(tmp.data$average_ror))
  
  sigma.forecast  <- ugarchforecast(tmp.egarch22, n.ahead = 1)
  sigma.forecast2 <- sigma.forecast@forecast$sigmaFor[1, 1]
  VaR2[k - start + 1] <- q01 * sigma.forecast2
}

Przypisujemy oszacowania do zbioru

EWP2$VaR2 <- VaR2

Wyres zwrotów wartości narażonej na ryzyko w okresie OUT-OF-SAMPLE

Narysowano wykres zwrotów wartości narażonej na ryzyko w okresie out-of-sample i jest on oznaczony zieloną linią. Czerwona linia pokazuje natomiast rzeczywisty przebieg badanego szeregu.

sum(EWP2$average_ror < EWP2$VaR2) / length(EWP2$VaR2)
## [1] 0

Dla modelu GARCH (1,1) straty ani razu nie przekroczyły założonego poziomu VaR.

TGARCH (1,2)

Wyres zwrotów wartości narażonej na ryzyko w okresie OUT-OF-SAMPLE 

sum(EWP2$average_ror < EWP2$VaR2) / length(EWP2$VaR2)
## [1] 0

Dla modelu TGARCH (1,2) straty ani razu nie przekroczyły założonego poziomu VaR.

EGARCH (1,1)

Wyres zwrotów wartości narażonej na ryzyko w okresie OUT-OF-SAMPLE 

sum(EWP2$average_ror < EWP2$VaR2) / length(EWP2$VaR2)
## [1] 0

Dla modelu EGARCH (1,1) straty ani razu nie przekroczyły założonego poziomu VaR.

Wszystkie modele okazały się równie restrykcyjne jeśli chodzi o oszacowania out-of-sample i w żadnym wypadku straty nie przekroczyły założonego poziomu VaR.

Podsumowanie

Zgodnie z wytycznymi badaniu poddano model GARCH i dwa modele będące rozszerzeniem tego modelu, w tym przypadku były to modele TGARCH oraz EGARCH. Model TGARCH miał za zadanie wychwycenie asymetrii szeregu, natomiast model EGARCH miał za zadanie wychwycić występujący w szeregu zwrotów efekt dźwigni, czyli asymetrycznego wpływu informacji na poziom przyszłej wariancji. Na początku przeprowadzono analizę badanego szeregu. Wstępna analiza pokazała, że badany szereg nie pochodzi z rozkładu normalnego i charakteryzuje się leptokurtycznością. Można było zaobserować autokorelację oraz efekty ARCH. Badane modele w większosci poradziły sobie z występującymi problemami. Dla okresu in sample najgorsze oszacowania dał model TGARCH(1,2), ponieważ wykazał najwięcej przekroczeń strat w stosunku do oszacowanej wartości narażonej na ryzyko. W okresie out-of-sample wszystkie modele okazały sie bardzo restrykcyjne. W żadnym z badanym modeli nie było przekroczeń poziomu VaR.