Streszczenie

Dana praca opiera się na analizie portfela równoważonego. Analizę tą wykonano przy pomocy modeli GARCH. Wybrano EGARCH i TGARCH, które posłużyły do oszacowania modeli w okresie IN-SAMPLE i OUT-OF-SAMPLE. Ważnym krokiem początkowym był odpowiedni dobór danych i wstępna ich analiza. Ostatecznie wyniki IN-SAMPLE potwierdziły poprawność badanych modeli, zaś wynki dla okresu OUT-OF-SAMPLE poprawnie pokazały prognozę.

Wstęp

Praca ta prezentuje analizę portfela równoważonego, czyli takiego, do którym każdy z instrumentów wchodzi z tą samą wagą (1/5). Pozwala to na rozważenie wszystkich instrumentów w tym samym stopniu. Związku z tym celem tej analizy było przeprowadzenie badania dla takiego właśnie portfela przy pomocy funkcji wariancji warunkowej w modelach klasy GARCH. Dla dwóch rozszeżenień tej kalsy modeli należało wykonać oszacowania warunkowej wariancji oraz wartości narazonej na ryzyko (VaR) w okresie in-sample, jak również w okresie out-of-sample. Badanie należało rozpocząć od badania danych w okresie in-sample, ponieważ pokaże nam to czy modele zostały dobrze dopasowane, natomiast okres out-of-sample sluży do oceny otrzymanej prognozy.

Uzasadanienie danych i metodologia

Do badania użyto 5 instrumentów finansowych:

  1. Ceny akcji spółki KGHM

  2. Wartość ceny surowca: miedzi

  3. Arbitraż walutowy: dolar kanadyjski- dolar amerykański

  4. Arbitraż walutowy: funt brytyjski- dolar kanadyjski

  5. Arbitraż walutowy: dolar amerykański- polski złoty

Wszystkie dane wykorzystano ze strony https://stooq.pl. Zdecydowano się na analizę rynku miedzi ze względu na dużą stablizację. Polska spółka KGHM ma dobry kurs i wiele firm chce w niż inwestować. Według Banku Światowego ceny surowców mają rosnąć. Bardzo ważnym aspektem dla polskiej spółki są kursy walutowe, wysoki kurs dolara w stosunku do polskiego złotego powoduje, że nawet przy spadku cen surowca Polska firma radzi sobie dobrze. Dzięki czemu spółka jest płynna i inwestorzy liczą na jej dalszy wzrost. Patrząc na ceny akcji zarówno spółki jak i wartości surowca zdecydowano, że najlepszym okresem do badnia będą dane od 26 czerwca 2012 roku do 26 czerwca 2017 roku. Od 2012 roku obserwowano spadek wartości akcji miedzi, gdzie najniższą w tym czasie wartośc odnotowano na początku stycznia 2016 roku. Notowania w 2016 roku oscylowały na niskim poziomie, zaobserwowano jednak spore wahania. Jednak przełomową datą dla akcji miedzi był październik 2016 roku, od tego momentu notowania poszły w górę. Ważnym aspektem przy wybraniu tego okresu czasowego jest również fakt, że z roku na rok rośnie ilość wydobywanej miedzi, ale w tym samym czasie więcej miedzi powstaje z powtórnej produkcji. Jakoś próbki po przerobieniu jej przy pierwszym podejściu i po recyklingu jest niemal nie do odróżnienia. Co jest zaletą tego surowca. Polepszanie się technologi i wzrost znaczenia miedzi na światowym rynku sprawia, że wykorzystywana jest w coraz większej ilości przemysłów. Jednym z takich przemysłów jest przemysł energetyczny. Do napędzania hybrydowych samochodów zaczęto wykorzystywać miedź, co jest jednym z powodów wzrostu znaczenia i cen akcji. Sugerując się prognozami Banku Światowego zakładamy, że druga połowa 2017 roku odnotuje wzrost cen akcji.

Przygotowanie danych do analizy

Analizę należy rozpocząć od wczytania danych:

url1 <- "https://stooq.pl/q/d/l/?s=usdpln&i=d" # usdpln
url2 <- "https://stooq.pl/q/d/l/?s=gbpcad&i=d" #gbpcad
url3 <- "https://stooq.pl/q/d/l/?s=cadusd&i=d" #cadusd 
url4 <- "https://stooq.pl/q/d/l/?s=kgh&i=d" # KGHM Polska Miedź SA 
url5 <- "https://stooq.pl/q/d/l/?s=hg.f&i=d" #miedź

Jednak aby z danych można było poprawnie korzystać należy je odpowiednio przygotować. Poczynając od wczytania danych z wyżej wymienionych stron internetowych i utworzenie z nich tabel. Z stooq.pl wczytają nam się dane dużo strasze niż są potrzebne, z tego powodu należało wybrać jedynie potrzebny do badania okres czasu. Aby przejść do wstępnej analizy danych, należało na początku policzyć stopy zwrotu dla każdego instrumentu (przedstawione poniżej), by na końcu przygotowywania danych policzyć średnią ze wszystkich zwrotów i na tej średniej wykonać badanie.

Poniżej przedstawione zostały wykresy dla ciągłych stóp zwrotu dla każdego instrumentu. Krok ten jest konieczny do wykonania modeli GARCH.

#pierwsze roznice dla arbitrażu USD-PLN
data1$r <- diff.xts(log(data1$Zusdpln))
# obliczamy ciagle stopy zwrotu
plot(data1$Data, data1$r, type = "l", main = "Wykres ciągłych stóp zwrotu dla arbitrażu USD-PLN")

Powyższy wykres przestawia ciągłych stóp zwrotu dla arbitrażu USD-PLN. Obserwuje się, że ciągłe stopy zwrotu wykazują duże wahania, a nawet zdarzają się obserwacje, które znacząco odbiegają od przeciętnej. Widoczna jest również asymetria spadków i wzrostów. Największą różnicę obserwuje się w połowie 2016 roku.

Powyższy wykres przestawia ciągłych stóp zwrotu dla arbitrażu GBP-CAD. Obserwuje się, że ciągłe stopy zwrotu wykazują duże wahania, a nawet zdarzają się obserwacje, które znacząco odbiegają od przeciętnej. W porównaniu do wykresu dla arbitrażu USD-PLN wahania są duże mniejsze i mniej obserwacji znacząco różni się od wartości przeciętnej. Można powiedzieć, że ten arbitraż jest bardziej stabilny i wiarygodny.

Na wykresie dla ciągłych stóp zwrotu dla arbitrażu CAD-USD widać duże wahania, nie ma natomiast znacząco odstających obserwacji jak w przypadku pierwszych dwóch arbitraży.

Jeśli chodzi o wykres ciągłych stóp zwrotu dla spółki KGHM, ponowanie obserwuje się duże wahania, jak również obserwacje, które znacząco odbiegają od wartości średnich. Największe wahania, czyli takie, dla których stopy zwrotu mają największą bądź najmniejszą wartość były zanotowane pod koniec 2016 roku. Co może pokazać, że już wtedy zaczynała się zmieniać sytuacja na rynku tego surowca.

Wykres dla ciągłych stóp zwrotu dla miedzi prezentuje szybko zmieniającą się sytuację na rynku surowca. Widać duże wahania i obserwacje odstające. Jednakże takich obserwacji, które znacząco odbiegają od pozostałych nie jest dużo. W pewnych okresach(od 2015r do 2017r) czasu wahania mają podobną wartość (oddalają się od wartości przeciętnej o mniewięcej podobną wartość).

Wstępna analiza

Po uzyskaniu ciągłych stóp zwrotu i policzeniu średniej ze wszystkich instrumentów należało przedtsawić dane na korelogramach. Które przedstwiają graficznie i liczbowo funkcję autokorelacji dla kolejnych opóźnień. Interesują nas te opóźnienia, których “siła” korelacji jest największa.

par(mfrow = c(2, 1))
acf(Data$r,  lag.max = 48,
    lwd = 7,
    col = "dark green", na.action = na.pass)
pacf(Data$r, lag.max=48,
     lwd = 7,
     col = "dark green", na.action = na.pass)

par(mfrow = c(1, 1))

Na korelogramie ACF widać, że pierwsze i trzecie opóźnienie przekracza poziom ufności, natomiast na korelogramie PACF żadne z pierwszych 10 opóźnień nie przekracza poziomu ufności.

Ważne jest by sprawdzić statystyki opisowe danych wraz z histogramem, który połączony jest z gęstością rozkładu normalnego. Nie można powiedzieć, że dane mają rozkład normalny.Do poniższego histogramu dodano rozkład gęstości rozkłąu normalnego, który jest potwierdzeniem, że dane nie charakteryzują się tym rozkładem.

basicStats(Data$r)
##               X..Data.r
## nobs        1249.000000
## NAs            0.000000
## Minimum       -0.034125
## Maximum        0.026158
## 1. Quartile   -0.003359
## 3. Quartile    0.003494
## Mean          -0.000063
## Median         0.000002
## Sum           -0.078460
## SE Mean        0.000168
## LCL Mean      -0.000393
## UCL Mean       0.000268
## Variance       0.000035
## Stdev          0.005954
## Skewness      -0.175854
## Kurtosis       2.167697
hist(Data$r, prob = T, breaks = 100, main="Histogram")
curve(dnorm(x,
            mean = mean(Data$r, na.rm = T),
            sd = sd(Data$r, na.rm = T)),
      col = "darkblue", lwd = 2, add = TRUE)

Należało również policzyć qwantyle, które będą potrzebne do poźniejszej analizy. Wykres Normal Q-Q dla stóp zwrotu nie zawiera się w prostej wyznaczonej dla rozkładu normalnego. Odbiega on od prostej w wielu miejscach. Na wykresie widać również grube ogony, które oznaczają duże prawdopodobieństwa, że w bazie występują obserwacje odstające. Ważnym aspektem jest również brak symetrii ogonów, gdyż lewy ogon jest większy niż prawy, co oznacza, że częsciej występują ujemne stopy zwrotu.

##        1% 
## -2.767766

Bardzo ważnym testem jest test Jarque–Bera test, który sprawdza normalność zmiennej. Na podstawie otrzymanego wyniku stwierdzono, że należy odrzucić hipotezę o braku autokorelacji. Czyli nie obserwuje się rozkładłu normalnego.

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

Drugim ważnym testem jest test sprawdzający czy występują efekty ARCH. Patrząc na uzyskane wyniki odrzucono hipotezę o braku tych efektów, dzięki czemu można przejść do analizy modeli GARCH.

ArchTest(Data$r, lags = 5)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  Data$r
## Chi-squared = 62.451, df = 5, p-value = 3.783e-12

Modele GARCH:

Pierwszym użytym modelem był model EGARCH(1,3), sugerując się korelogramami badano modele o wartościach p i q należących do zbioru (1;3). Ważne było to, by parametr alpha1 był istotny, gdyż to on w przypadku tego modelu wyjaśnia efekt szoku(asymetria), natomiast parametr gamma1 określa jego wielkość. Przedstawiony poniżej wykres jak również wartości liczbowe potwierdzają istnienie asymetrii.

spec <- ugarchspec( variance.model = list(model ="eGARCH",
                        garchOrder = c(1, 3)),
  mean.model = list(armaOrder = c(0, 0), include.mean = F),
  distribution.model = "norm")
k.egarch13a <- ugarchfit( spec = spec,
                         data = na.omit(Data$r))
k.egarch13a
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : eGARCH(1,3)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## omega  -0.223421    0.019329  -11.5588  0.00000
## alpha1 -0.054604    0.014908   -3.6626  0.00025
## beta1   0.999996    0.000394 2534.9940  0.00000
## beta2  -0.521503    0.001689 -308.8109  0.00000
## beta3   0.499546    0.001800  277.4579  0.00000
## gamma1  0.156023    0.028528    5.4690  0.00000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## omega  -0.223421    0.011391   -19.6146 0.000000
## alpha1 -0.054604    0.024953    -2.1883 0.028649
## beta1   0.999996    0.000210  4761.1307 0.000000
## beta2  -0.521503    0.000313 -1665.4861 0.000000
## beta3   0.499546    0.000195  2561.2668 0.000000
## gamma1  0.156023    0.054936     2.8401 0.004510
## 
## LogLikelihood : 4683.782 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -7.5085
## Bayes        -7.4838
## Shibata      -7.5085
## Hannan-Quinn -7.4992
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      14.36 1.514e-04
## Lag[2*(p+q)+(p+q)-1][2]     14.43 1.296e-04
## Lag[4*(p+q)+(p+q)-1][5]     18.65 4.544e-05
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic  p-value
## Lag[1]                       10.28 0.001342
## Lag[2*(p+q)+(p+q)-1][11]     14.39 0.011565
## Lag[4*(p+q)+(p+q)-1][19]     17.26 0.040814
## d.o.f=4
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[5]     2.024 0.500 2.000  0.1548
## ARCH Lag[7]     2.812 1.473 1.746  0.3498
## ARCH Lag[9]     3.191 2.402 1.619  0.5306
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.2805
## Individual Statistics:              
## omega  0.26638
## alpha1 0.41177
## beta1  0.26195
## beta2  0.26200
## beta3  0.26268
## gamma1 0.09449
## 
## 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           1.6847 0.09229   *
## Negative Sign Bias  1.2220 0.22195    
## Positive Sign Bias  0.4721 0.63693    
## Joint Effect        5.5919 0.13324    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     26.84       0.1084
## 2    30     27.19       0.5612
## 3    40     44.02       0.2676
## 4    50     53.52       0.3050
## 
## 
## Elapsed time : 1.297074
plot(k.egarch13a,which=12)

Ważnym punktem jest test Ljunga-Box’a, gdzie testy: Ljung-Box Test R^2 Q(10), Ljung-Box Test R^2 Q(15), Ljung-Box Test R^2 Q(20) i test LM ARCH TEST R TR^2 by wartości p-Value dla nich były większe niż 0,05. Określają one czy występuje autokorelacja czy nie. Kwadraty wystaryzowanych reszt muszą być białym szumem. W tym przypadku model jest poprawny, gdyż nie występuje autokorelacja.

Można zauważyć, że utworzono 2 różne wartości k.egarch13 i k.egarch13a, powodem jest to, że jedna jest klasy fGarch, a druga uGARCHfit. Drugi wymieniony typ klasy jest potrzebny do funkcji, która oszacowuje model w okresie OUT-OF-SAMPLE. Bez tej czynności nie uda się wykonać operacji przedstawionej poniżej dla okresu OUT-OF-SAMPLE. Taką samą operację wykonano również i dla drugeigo modelu.

k.egarch13 <- garchFit(formula = ~ garch(1, 3),
                       data = na.omit(Data$r),
                       include.mean = F,
                       cond.dist = "norm",
                       trace = F)
summary(k.egarch13)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 3), data = na.omit(Data$r), cond.dist = "norm", 
##     include.mean = F, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 3)
## <environment: 0x0000000017be74b8>
##  [data = na.omit(Data$r)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1       beta2       beta3  
## 1.1697e-06  1.0530e-01  4.0218e-01  8.7413e-03  4.5166e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  1.170e-06   5.254e-07    2.226   0.0260 *  
## alpha1 1.053e-01   2.597e-02    4.056    5e-05 ***
## beta1  4.022e-01   2.591e-01    1.553   0.1205    
## beta2  8.741e-03   2.124e-01    0.041   0.9672    
## beta3  4.517e-01   2.100e-01    2.150   0.0315 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  4678.966    normalized:  3.75519 
## 
## Description:
##  Fri Jun 30 10:53:03 2017 by user: Admin 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  102.1617  0           
##  Shapiro-Wilk Test  R    W      0.9893917 7.628871e-08
##  Ljung-Box Test     R    Q(10)  25.63617  0.004261601 
##  Ljung-Box Test     R    Q(15)  30.53508  0.01013179  
##  Ljung-Box Test     R    Q(20)  34.68074  0.02187209  
##  Ljung-Box Test     R^2  Q(10)  10.13523  0.4287103   
##  Ljung-Box Test     R^2  Q(15)  12.42352  0.6467306   
##  Ljung-Box Test     R^2  Q(20)  18.41042  0.5603914   
##  LM Arch Test       R    TR^2   10.46931  0.5748597   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -7.502353 -7.481777 -7.502385 -7.494617

Celem poniższego kodu jest policzenie bezwarunkowej wariancji i zestawienie jej na wykresie z warunkową wariancją. Można zauważyć, że warunkowa wariancja w dłuższym okresie dąży do bezwarunkowej wariancji.

var_uncond1 <- k.egarch13@fit$par[1] / (1 - k.egarch13@fit$par[2]
                                       - k.egarch13@fit$par[3]
                                       - k.egarch13@fit$par[4]
                                       - k.egarch13@fit$par[5])
names(var_uncond1) <- "unconditional variance"
var_uncond1
## unconditional variance 
##           3.642238e-05
k.fore200 <- predict(k.egarch13, n.ahead = 200)
plot(k.fore200[, 3] ^ 2, type = "l")
abline(h = var_uncond1, col = "red", lty = 2)
title(main = "Warunkowa i bezwarunkowa wariancja zwrotow portfela")

Następnie policzono VaR w okresie IN-SAMPLE i zademonstrowano jego wykres. Zielona linia określa wyliczoną wartość VaRa. Sprawdzamy ile razy czerwona linia przekroczyła zielony wykres, “na oko” nie jest to wiarygodne badanie, z tego powodu poniżej została zademonstrowana funkcja, która odnotuje ile takich zdarzeń miało miejsce. Okazało się, że 11 razy straty przekroczyły zakładany poziom VaR.

Data$VaR1 <- q01 * k.egarch13@sigma.t
plot(Data$Data, Data$r, col = "red", lwd = 1, type = 'l',ylim = c(-0.1, 0.1))
abline(h = 0, lty = 2)
lines(Data$Data, Data$VaR1, type = 'l', col = "green")

sum(Data$r < Data$VaR1)
## [1] 11

Następnie policzono VaR w okresie OUT-OF-SAMPLE, jednodniową prognozę warunkowego odchylenia standardowego. Czerwona część wykresu przedstwia prognozę na 1 dzień , jak widać przewiduje ona dalszy spadek.

plot(ugarchforecast(k.egarch13a, n.ahead = 1), which=3)

sigma.Forecast <- ugarchforecast(k.egarch13a, n.ahead = 1)
sigma.Forecast@forecast$sigmaFor
##     1973-05-31 01:00:00
## T+1         0.006163405
sigma.Forecast2 <- sigma.Forecast@forecast$sigmaFor[1, 1]

Przy pomocy poniższej funkcji wyliczono jednodniową wartość VaR, dzięki temu wiadomo, że strata nie może być większa niż -0.01705886.

q01 * sigma.Forecast2
##          1% 
## -0.01705886

W kolejnym kroku oszacowoano jednodniowego VaRa dla całego okresu OUT-OF-SAMPLE. Poniżej przedstawiona pętla “for” wylicza VaRy dla całego okresu, zaś funkcja umiejscowniona zaraz za pętlą dopisuje te wartości do bazy.

Data$obs<-1:length(Data$Data)
start  <- Data$obs[Data$Data == as.Date("2015-06-26")]
finish <- Data$obs[Data$Data == as.Date("2017-06-26")]
DaTa <-Data[start:finish, ]
VaR <- rep(NA, times = finish - start + 1)

for (k in start:finish) {
  tmp.data <- Data[Data$obs <= (k - 1), ]
  tmp.data <- tmp.data[as.Date("2012-06-26") <= tmp.data$Data, ]
  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 = "eGARCH",
                                           garchOrder = c(1, 3)),
                     mean.model = list(armaOrder = c(0, 0),
                                       include.mean = T),
                     distribution.model = "norm")
  tmp.garch13a <- ugarchfit(spec = spec, data = na.omit(tmp.data$r))
  sigma.Forecast  <- ugarchforecast(tmp.garch13a, n.ahead = 1)
  sigma.Forecast2 <- sigma.Forecast@forecast$sigmaFor[1, 1]
  VaR[k - start + 1] <- q01 * sigma.Forecast2
}
DaTa$VaR <- VaR

Na koniec zrobiono wykres zwrotów w porównaniu do VaRów dla okresu OUT-OF-SAMPLE i policzono w ilu przypadkach straty przekroczyły zakładny poziom VaR- wyszło ich 4. Co pokazuje, że okres OUT-OF-SAMPLE jest bardziej restrykcyjny niż okres IN-SAMPLE, gdyż tam 11 obserwacji przekroczyło poziom VaR.

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

sum(DaTa$r < DaTa$VaR)
## [1] 4

Model TGARCH (1,3) Drugim wybranym modelem jest TGARCH, w tym modelu bardzo ważnym parametrem jest parametr etta, który określa czy w modelu występuje asymetria czy nie. Patrząc na wyniki istotności tego parametru jak również na poniższy wykres można stwierdzić, że asymetria występuje w badanym modelu.

spec <- ugarchspec(
  # rownanie warunkowej wariancji
  variance.model = list(model = "fGARCH",
                        garchOrder = c(1, 3),
                        submodel = "TGARCH"),
  mean.model = list(armaOrder = c(0, 0),
                    include.mean = F),
  distribution.model = "norm")
k.tgarch13a <- ugarchfit(spec = spec, data = na.omit(Data$r))
k.tgarch13a
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : fGARCH(1,3)
## fGARCH Sub-Model : TGARCH
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## omega   0.000144    0.000064 2.262460 0.023669
## alpha1  0.087874    0.020655 4.254283 0.000021
## beta1   0.622146    0.081662 7.618537 0.000000
## beta2   0.000001    0.225301 0.000002 0.999998
## beta3   0.285253    0.072229 3.949266 0.000078
## eta11   0.331725    0.125963 2.633522 0.008450
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## omega   0.000144    0.000115 1.252062 0.210547
## alpha1  0.087874    0.047981 1.831449 0.067034
## beta1   0.622146    0.079931 7.783497 0.000000
## beta2   0.000001    0.119722 0.000004 0.999997
## beta3   0.285253    0.056535 5.045617 0.000000
## eta11   0.331725    0.316879 1.046852 0.295168
## 
## LogLikelihood : 4682.221 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -7.5060
## Bayes        -7.4813
## Shibata      -7.5060
## Hannan-Quinn -7.4967
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      13.82 2.010e-04
## Lag[2*(p+q)+(p+q)-1][2]     13.87 1.818e-04
## Lag[4*(p+q)+(p+q)-1][5]     17.96 7.022e-05
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic  p-value
## Lag[1]                       10.19 0.001411
## Lag[2*(p+q)+(p+q)-1][11]     15.63 0.006086
## Lag[4*(p+q)+(p+q)-1][19]     18.66 0.022813
## d.o.f=4
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[5]     1.972 0.500 2.000  0.1602
## ARCH Lag[7]     2.724 1.473 1.746  0.3642
## ARCH Lag[9]     3.225 2.402 1.619  0.5246
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.1249
## Individual Statistics:             
## omega  0.2116
## alpha1 0.2365
## beta1  0.2352
## beta2  0.2358
## beta3  0.2318
## eta11  0.4099
## 
## 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           1.8359 0.06661   *
## Negative Sign Bias  1.2845 0.19920    
## Positive Sign Bias  0.3161 0.75195    
## Joint Effect        5.8338 0.11998    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     24.47       0.1789
## 2    30     24.88       0.6843
## 3    40     36.57       0.5813
## 4    50     56.09       0.2264
## 
## 
## Elapsed time : 0.7070398
plot(k.tgarch13a, which=12)

Ponownie ważnym punktem jest test Ljunga-Box’a, gdzie testy: Ljung-Box Test R^2 Q(10), Ljung-Box Test R^2 Q(15), Ljung-Box Test R^2 Q(20) i test LM ARCH TEST R TR^2 by wartości p-Value dla nich były większe niż 0,05. Określają one czy występuje autokorelacja czy nie. Kwadraty wystaryzowanych reszt muszą być białym szumem. Oszacowany model jest poprawny, gdyż nie występuje autokorelacja.

k.tgarch13 <- garchFit(formula = ~ garch(1, 3),
                       data = na.omit(Data$r),
                       include.mean = F,
                       cond.dist = "norm",
                       trace = F)
k.tgarch13@fit$par
##        omega       alpha1        beta1        beta2        beta3 
## 1.169684e-06 1.053044e-01 4.021824e-01 8.741314e-03 4.516575e-01
summary(k.tgarch13)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 3), data = na.omit(Data$r), cond.dist = "norm", 
##     include.mean = F, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 3)
## <environment: 0x0000000043598bb8>
##  [data = na.omit(Data$r)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1       beta2       beta3  
## 1.1697e-06  1.0530e-01  4.0218e-01  8.7413e-03  4.5166e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  1.170e-06   5.254e-07    2.226   0.0260 *  
## alpha1 1.053e-01   2.597e-02    4.056    5e-05 ***
## beta1  4.022e-01   2.591e-01    1.553   0.1205    
## beta2  8.741e-03   2.124e-01    0.041   0.9672    
## beta3  4.517e-01   2.100e-01    2.150   0.0315 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  4678.966    normalized:  3.75519 
## 
## Description:
##  Fri Jun 30 10:59:03 2017 by user: Admin 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  102.1617  0           
##  Shapiro-Wilk Test  R    W      0.9893917 7.628871e-08
##  Ljung-Box Test     R    Q(10)  25.63617  0.004261601 
##  Ljung-Box Test     R    Q(15)  30.53508  0.01013179  
##  Ljung-Box Test     R    Q(20)  34.68074  0.02187209  
##  Ljung-Box Test     R^2  Q(10)  10.13523  0.4287103   
##  Ljung-Box Test     R^2  Q(15)  12.42352  0.6467306   
##  Ljung-Box Test     R^2  Q(20)  18.41042  0.5603914   
##  LM Arch Test       R    TR^2   10.46931  0.5748597   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -7.502353 -7.481777 -7.502385 -7.494617

Dla tego modelu również wylicozno bezwarunkową wariancję i na wykresie pokazano jak warunkowa wariancja zachowuje się w stosunku do niej. Na poniższym wykresie można zauważyć, że warunkowa wariancja w dłuższym okresie(około 200 okresów) dąży do bezwarunkowej wariancji.

## unconditional variance 
##           3.642238e-05

Następnie policzono VaR w okresie IN-SAMPLE i zademonstrowano jego wykres. Zielona linia określa wyliczoną wartość VaRa. Sprawdzamy ile razy czerwona linia przekroczyła zielony wykres, “na oko” nie jest to wiarygodne badanie, z tego powodu poniżej została zademonstrowana funkcja, która odnotuje ile takich zdarzeń miało miejsce. Okazało się, że 11 razy straty przekroczyły zakładany poziom VaR. Zauważa się, że oba modele w podobny sposób oceniają ryzyko, gdyż w obu modelach 11 razy straty przekroczyły wartość VaR.

Data$VaR2 <- q01 * k.tgarch13@sigma.t
plot(Data$Data, Data$r, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.1, 0.1))
abline(h = 0, lty = 2)
lines(Data$Data, Data$VaR2, type = 'l', col = "green")

sum(Data$r < Data$VaR2)
## [1] 11

Następnie policzono VaR w okresie OUT-OF-SAMPLE, jednodniową prognozę warunkowego odchylenia standardowego. Czerwona część wykresu przedstwia prognozę na 1 dzień , jak widać przewiduje ona dalszy spadek.

plot(ugarchforecast(k.tgarch13a, n.ahead = 1), which=3)

sigma.forecast <- ugarchforecast(k.tgarch13a, n.ahead = 1)
sigma.forecast@forecast$sigmaFor
##     1973-05-31 01:00:00
## T+1         0.006229326
Sigma.forecast <- sigma.forecast@forecast$sigmaFor[1, 1]

Obliczamy jednodniową wartość VaR, która wynosi: -0.0172327.

q01 * Sigma.forecast
##         1% 
## -0.0172327

Również dla tego modelu wykonywana jest powyżej przedstawiona pętla “for”, wykonująca ten sam krok co w przypadku modelu EGARCH.

Natępnie zademonstrowano wykres zwrotów w porównaniu do policzonych wartości VaR w okresie OUT-OF-SAMPLE. Ponowanie interesuje nas ile razy czerwony wykres przekroczył wartości zielonego wykresu. By to sprawdzić wykorzystano poniższą funkcję, okazało się, że było 4 takich oserwacji. Co oznacza, że ponownie oba modele w okresie OUT-OF-SAMPLE tak samo oceniają wartość VaR.

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

sum(DaTa$r < DaTa$VaR)
## [1] 4

Podsumowanie

Wstępna analiza danych pokazała, że dane nie mają rozkładu normalnego, ale występują w nich efekty ARCH, co w dalszej perspektywnie było ważne. W dlaszej analizie oszacowano dwa modele typu GARCH w okresie in-sample i out-of-sample. Modele w okresie in-sample potwierdziły ich zgodność, oba modele wykazły asymetrię i reszty okazały się białym szumem. Potwierdza to, że modele zostały poprawnie wybrane. W obu modelach straty 11 razy przekroczyły poziom VaR, patrząc ilość obserwacji w bazie można powiedzieć, że jest to mało i modele są dość restrykcyjne. Natomiast jeśli chodzi o okres OUT-OF-SAMPLE to na wykresach obserwuje się, że przewidziany spadek na kolejny dzień. Czyli początkowe założenie, że obserwować będziemy wzrost nie okazało się prawdą. Jednak oba modele okazały się restrykcyjne i jedynie 4 obserwacji przekroczyło poziom VaR. Cel projektu został osiągnięty, przebadano portfel równoważony dla surowca miedzi i oszacowano go w dwóch okresach: IN-SAMPLE i OUT-OF-SAMPLE. Warto zauważyć, że oba modele dały zbliżone rezultaty, co potwierdza poprawne oszacowanie danych.