Začneme s inštaláciou knižníc potrebných na toto cvičenie. Spustením nasledovného príkazu sa nainštalujú tie z knižníc, ktoré nainštalované nie sú, a zároveň sa aj načítajú.
# názvy balíkov, ktoré budú potrebné pre cvičenie 6
packages <- c("astsa", "urca", "ggplot2", "quantmod", "fGarch")
# inštalácia balíkov, ak nie sú nainštalované
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages],dependencies=TRUE)
}
# načítanie balíkov
invisible(lapply(packages, library, character.only = TRUE))
V prvej časti budeme, tak ako na prednáške, pracovať s logaritmickými výnosmi akcií. Konkrétne sa pozrieme na výnosy akcie AAPL (Apple):
getSymbols("AAPL",
from = "2019-01-01",
to = "2024-12-01",
auto.assign = TRUE)
## [1] "AAPL"
Budeme pracovať s týždennými dátami.
AAPL.tyzden <- to.weekly(AAPL)
head(AAPL.tyzden)
## AAPL.Open AAPL.High AAPL.Low AAPL.Close AAPL.Volume AAPL.Adjusted
## 2019-01-04 38.7225 39.7125 35.500 37.0650 747836000 35.40196
## 2019-01-11 37.1750 38.6325 36.475 38.0725 814824400 36.36425
## 2019-01-18 37.7125 39.4700 37.305 39.2050 621168000 37.44595
## 2019-01-25 39.1025 39.5325 37.925 39.4400 450006400 37.67039
## 2019-02-01 38.9475 42.2500 38.415 41.6300 809187200 39.76213
## 2019-02-08 41.8525 43.8925 41.820 42.6025 605593600 40.86552
Transformácia cien na logaritmické výnosy:
ceny <- AAPL.tyzden$AAPL.Adjusted
vynosy <- diff(log(ceny))
vynosy <- vynosy[-1]
Vykreslenie dát
chartSeries(vynosy, theme=chartTheme('white',up.col="#009E73"))
Skontrolujeme, či v dátch náhodou nie je jednotkový koreň
summary(ur.df(vynosy, lags = 5, selectlags = "BIC", type = "none"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.194704 -0.016046 0.006969 0.027720 0.137304
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.96414 0.08172 -11.798 <2e-16 ***
## z.diff.lag -0.03545 0.05774 -0.614 0.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03928 on 300 degrees of freedom
## Multiple R-squared: 0.4999, Adjusted R-squared: 0.4966
## F-statistic: 149.9 on 2 and 300 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -11.798
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
Hypotézu o jednotkovom koreni zamietame, takže dáta nebude potrebné differencovať.
Vykreslíme si ACF a PACF pre naše dáta
acf2(vynosy)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF -0.02 0.02 0.1 -0.04 -0.01 -0.04 -0.05 -0.08 0.00 -0.03 0.03 -0.07 0.00
## PACF -0.02 0.02 0.1 -0.03 -0.01 -0.05 -0.04 -0.08 0.01 -0.02 0.04 -0.08 -0.01
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF 0.01 0.05 -0.01 0 0.04 -0.12 -0.07 0.06 -0.07 -0.12 0.04 0.08
## PACF 0.00 0.06 -0.02 0 0.02 -0.12 -0.10 0.07 -0.04 -0.11 0.02 0.10
## [,26] [,27] [,28]
## ACF -0.01 -0.02 0.02
## PACF -0.01 -0.06 0.00
Zdá sa, že podľa ACF aj PACF by sme dáta mohli modelovať ako konštanta + biely šum.
vynosy.00 <- sarima(vynosy,0,0,0)
## initial value -3.256473
## iter 1 value -3.256473
## final value -3.256473
## converged
## initial value -3.256473
## iter 1 value -3.256473
## final value -3.256473
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## xmean 0.0062 0.0022 2.8133 0.0052
##
## sigma^2 estimated as 0.001484102 on 307 degrees of freedom
##
## AIC = -3.662081 AICc = -3.662039 BIC = -3.63786
##
vynosy.00
## $fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## xreg = xmean, include.mean = FALSE, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## xmean
## 0.0062
## s.e. 0.0022
##
## sigma^2 estimated as 0.001484: log likelihood = 565.96, aic = -1127.92
##
## $degrees_of_freedom
## [1] 307
##
## $ttable
## Estimate SE t.value p.value
## xmean 0.0062 0.0022 2.8133 0.0052
##
## $ICs
## AIC AICc BIC
## -3.662081 -3.662039 -3.637860
Porovnáme si ACF pre rezídua a pre ich druhé mocniny:
acf1(vynosy.00$fit$residuals)
## [1] -0.02 0.02 0.10 -0.04 -0.01 -0.04 -0.05 -0.08 0.00 -0.03 0.03 -0.07
## [13] 0.00 0.01 0.05 -0.01 0.00 0.04 -0.12 -0.07 0.06 -0.07 -0.12 0.04
## [25] 0.08 -0.01 -0.02 0.02
acf1(vynosy.00$fit$residuals^2)
## [1] 0.12 0.03 0.35 0.05 -0.01 0.08 0.04 -0.05 -0.02 -0.02 -0.03 -0.02
## [13] -0.03 -0.02 -0.05 0.00 -0.04 0.02 0.18 -0.01 -0.02 0.15 0.02 0.02
## [25] 0.10 0.01 -0.03 -0.03
Vidíme, že hoci samotné rezídua nevyzerajú byť korelované, ich druhé mocniny korelované sú. Takúto vlastnosť ale biely šum nemá. Preto budeme \(u\) modelovať ako GARCH(p,q) proces.
Najskôr budeme dáta modelovať iba ako ARCH(p) proces. Začneme základným ARCH(1) procesom.
Otázky:
V R-ku na modelovanie použijeme funkciu garchFit
z
knižnice fGarch
:
~garch(1,0)
- predstavuje formulu, všeobecne v tvare
~arma(p,q)+garch(p,q)
data
- naše dátatrace
- ak chceme vypísať konvergenciu, atď… (FALSE ak
nechceme)vynosy.10 <- garchFit(~garch(1,0), data = vynosy, trace = FALSE)
vynosy.10
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 0), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 0)
## <environment: 0x000001d2a8519790>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1
## 0.0085720 0.0010083 0.3977677
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0085720 0.0020253 4.232 2.31e-05 ***
## omega 0.0010083 0.0001309 7.702 1.33e-14 ***
## alpha1 0.3977677 0.1342793 2.962 0.00305 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 573.3349 normalized: 1.861477
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
Z výstupu môžeme vidieť odhady pre jednotlivé parametre:
mu
- konštanta v rovnici pre strednú hodnotuomega
- konštanta v rovnici pre disperziualpha1
- koeficient pri \(u^2_{t-1}\)Spolu s odhadmi sú dané aj štandardné odchýlky a hodnota testovej štatistiky (H0: daný koeficient je nulový) spolu s p-value.
Všetky tri koeficienty sú signifikantné.
Pozrieme sa ešte na zhodnotenie rezíduí. Dobrý model bude taký, kde
nezostane žiadna autokorelácia v rezíduách ani v ich druhých mocninách a
zároveň hypotéza o homoskedasticite sa nebude zamietať. K týmto hodnotám
pristúpime cez funkciu summary
, ktorá nám ponúka aj hodnotu
informačných kritérií pre daný model. K tým budeme ale pristupovať
neskôr.
summary(vynosy.10)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 0), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 0)
## <environment: 0x000001d2a8519790>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1
## 0.0085720 0.0010083 0.3977677
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0085720 0.0020253 4.232 2.31e-05 ***
## omega 0.0010083 0.0001309 7.702 1.33e-14 ***
## alpha1 0.3977677 0.1342793 2.962 0.00305 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 573.3349 normalized: 1.861477
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 39.9894 2.072105e-09
## Shapiro-Wilk Test R W 0.9823208 0.0007576804
## Ljung-Box Test R Q(10) 6.261832 0.7928064
## Ljung-Box Test R Q(15) 8.93209 0.8810414
## Ljung-Box Test R Q(20) 16.15914 0.7067037
## Ljung-Box Test R^2 Q(10) 33.62831 0.0002134665
## Ljung-Box Test R^2 Q(15) 35.54305 0.002056662
## Ljung-Box Test R^2 Q(20) 40.44098 0.004391653
## LM Arch Test R TR^2 31.46994 0.001669588
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -3.703474 -3.667141 -3.703661 -3.688946
Vidíme, že stále ostáva istá autokorelácia medzi druhými mocninami rezíduí, aj keď samotné rezíduá nie sú korelované. Rovnako sa zamieta aj LM Arch test, teda zamietame hypotézu o homoskedasticite (konštantnosti disperzie v rezíduach).
Úloha: Zobrazte ACF pre štandardizované rezíduá a pre ich druhé mocniny.
Z ACF štvorcov štandardizovaných rezíduí vidíme, že ARCH(1) nie je dobrý model. Skúsime teda zvýšiť rád ARCH procesu na ARCH(2):
vynosy.20 <- garchFit(~garch(2,0), data = vynosy, trace = FALSE)
vynosy.20
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(2, 0), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(2, 0)
## <environment: 0x000001d2b050ddd0>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2
## 0.00866803 0.00093114 0.38156106 0.07229835
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0086680 0.0020605 4.207 2.59e-05 ***
## omega 0.0009311 0.0001562 5.961 2.51e-09 ***
## alpha1 0.3815611 0.1338190 2.851 0.00435 **
## alpha2 0.0722983 0.0935662 0.773 0.43970
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 573.5904 normalized: 1.862306
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
summary(vynosy.20)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(2, 0), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(2, 0)
## <environment: 0x000001d2b050ddd0>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2
## 0.00866803 0.00093114 0.38156106 0.07229835
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0086680 0.0020605 4.207 2.59e-05 ***
## omega 0.0009311 0.0001562 5.961 2.51e-09 ***
## alpha1 0.3815611 0.1338190 2.851 0.00435 **
## alpha2 0.0722983 0.0935662 0.773 0.43970
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 573.5904 normalized: 1.862306
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 38.80513 3.746048e-09
## Shapiro-Wilk Test R W 0.9824344 0.0007981876
## Ljung-Box Test R Q(10) 6.457223 0.775499
## Ljung-Box Test R Q(15) 9.205291 0.866541
## Ljung-Box Test R Q(20) 16.66369 0.674691
## Ljung-Box Test R^2 Q(10) 33.21241 0.0002508757
## Ljung-Box Test R^2 Q(15) 35.08405 0.002392165
## Ljung-Box Test R^2 Q(20) 40.54755 0.004256481
## LM Arch Test R TR^2 30.39293 0.00243602
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -3.698639 -3.650196 -3.698971 -3.679269
Vidíme, že koeficient alpha2
vyšiel nesignifikantný,
zároveň sa však zachovali nedostatky rezíduí z ARCH(1) procesu (v
niektorých prípadoch dokonca horšie p hodnoty).
Skúsme ARCH(3):
vynosy.30 <- garchFit(~garch(3,0), data = vynosy, trace = FALSE)
vynosy.30
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(3, 0), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(3, 0)
## <environment: 0x000001d2b2939450>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 alpha3
## 0.00729292 0.00090668 0.21858537 0.00000001 0.15916432
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 7.293e-03 2.096e-03 3.479 0.000504 ***
## omega 9.067e-04 1.414e-04 6.411 1.45e-10 ***
## alpha1 2.186e-01 9.999e-02 2.186 0.028808 *
## alpha2 1.000e-08 6.392e-02 0.000 1.000000
## alpha3 1.592e-01 6.470e-02 2.460 0.013891 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 581.4303 normalized: 1.887761
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
summary(vynosy.30)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(3, 0), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(3, 0)
## <environment: 0x000001d2b2939450>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 alpha3
## 0.00729292 0.00090668 0.21858537 0.00000001 0.15916432
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 7.293e-03 2.096e-03 3.479 0.000504 ***
## omega 9.067e-04 1.414e-04 6.411 1.45e-10 ***
## alpha1 2.186e-01 9.999e-02 2.186 0.028808 *
## alpha2 1.000e-08 6.392e-02 0.000 1.000000
## alpha3 1.592e-01 6.470e-02 2.460 0.013891 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 581.4303 normalized: 1.887761
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 15.44972 0.0004417083
## Shapiro-Wilk Test R W 0.9885061 0.01539655
## Ljung-Box Test R Q(10) 6.680395 0.7552337
## Ljung-Box Test R Q(15) 10.51955 0.7858498
## Ljung-Box Test R Q(20) 17.02062 0.6516341
## Ljung-Box Test R^2 Q(10) 4.958119 0.8939581
## Ljung-Box Test R^2 Q(15) 6.577033 0.9683427
## Ljung-Box Test R^2 Q(20) 8.100123 0.9911844
## LM Arch Test R TR^2 4.664811 0.9682484
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -3.743054 -3.682501 -3.743570 -3.718842
Opäť máme nesignifikantný koeficient alpha2
, no druhé
mocniny rezíduí sú už bez výraznej autokorelácie, LM Arch test sa tiež
nezamietol, teda môžeme predpokladať konštatnú disperziu v
rezíduách.
ARCH(3) proces vyšiel ako vhodný model pre naše dáta, avšak vyšiel
nám v ňom koeficient alpha2
nesignifikantný. Preto sa
pozrieme aj na zovšeobecnený ARCH proces, GARCH. Skúsime začať
základným, ktorý sa vo veľa prípadoch ukazuje ako dostatočný na
modelovanie volatility. Je ním GARCH(1,1) proces.
Otázky:
Čo znamená, ak \(u\) nemodelujeme ako biely šum, ale ako GARCH(p,q) proces? Ako vtedy modelujeme jeho disperziu?
Aké ohraničenia na parametre musia byť splnené, aby bol GARCH(p,q) proces stacionárny?
Modelujme v R-ku:
vynosy.11 <- garchFit(~garch(1,1), data = vynosy, trace = FALSE)
vynosy.11
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x000001d2b3d76fb8>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## 0.00829601 0.00028659 0.22748633 0.59388094
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0082960 0.0020551 4.037 5.42e-05 ***
## omega 0.0002866 0.0001199 2.390 0.01686 *
## alpha1 0.2274863 0.0869563 2.616 0.00889 **
## beta1 0.5938809 0.1234115 4.812 1.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 579.4268 normalized: 1.881256
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
summary(vynosy.11)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x000001d2b3d76fb8>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## 0.00829601 0.00028659 0.22748633 0.59388094
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0082960 0.0020551 4.037 5.42e-05 ***
## omega 0.0002866 0.0001199 2.390 0.01686 *
## alpha1 0.2274863 0.0869563 2.616 0.00889 **
## beta1 0.5938809 0.1234115 4.812 1.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 579.4268 normalized: 1.881256
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 24.1597 5.672676e-06
## Shapiro-Wilk Test R W 0.9861115 0.004599973
## Ljung-Box Test R Q(10) 6.941273 0.7309772
## Ljung-Box Test R Q(15) 10.87408 0.7614685
## Ljung-Box Test R Q(20) 18.26838 0.5697341
## Ljung-Box Test R^2 Q(10) 12.60088 0.2468506
## Ljung-Box Test R^2 Q(15) 14.65495 0.4765452
## Ljung-Box Test R^2 Q(20) 18.14546 0.5778269
## LM Arch Test R TR^2 11.22305 0.5099066
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -3.736538 -3.688095 -3.736869 -3.717168
Všetky koeficienty vyšli signifikatné. Testy nezamietli nulovú koreláciu medzi rezíduami, ani medzi ich druhými mocninami, nezamietol sa ani LM Arch test pre homoskedasticitu v rezíduách. Teda, tak ako ARCH(3), aj GARCH(1,1) je pre naše dáta dobrým modelom.
Na ukážku ešte vyskúšame GARCH proces vyššieho rádu napríklad GARCH(1,2).
vynosy.12 <- garchFit(~garch(1,2), data = vynosy, trace = FALSE)
vynosy.12
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 2), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 2)
## <environment: 0x000001d2ac9cc138>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1 beta2
## 0.0077846 0.0003941 0.2733337 0.0158758 0.4529177
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0077846 0.0019951 3.902 9.55e-05 ***
## omega 0.0003941 0.0001900 2.074 0.03806 *
## alpha1 0.2733337 0.0900403 3.036 0.00240 **
## beta1 0.0158758 0.1304451 0.122 0.90313
## beta2 0.4529177 0.1546620 2.928 0.00341 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 581.7533 normalized: 1.88881
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
summary(vynosy.12)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 2), data = vynosy, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 2)
## <environment: 0x000001d2ac9cc138>
## [data = vynosy]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1 beta2
## 0.0077846 0.0003941 0.2733337 0.0158758 0.4529177
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0077846 0.0019951 3.902 9.55e-05 ***
## omega 0.0003941 0.0001900 2.074 0.03806 *
## alpha1 0.2733337 0.0900403 3.036 0.00240 **
## beta1 0.0158758 0.1304451 0.122 0.90313
## beta2 0.4529177 0.1546620 2.928 0.00341 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 581.7533 normalized: 1.88881
##
## Description:
## Wed Dec 11 16:20:05 2024 by user: Anna Hlubinová
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 15.46707 0.0004378928
## Shapiro-Wilk Test R W 0.9884898 0.01526751
## Ljung-Box Test R Q(10) 6.774191 0.7465769
## Ljung-Box Test R Q(15) 10.63116 0.7782761
## Ljung-Box Test R Q(20) 17.37188 0.6287104
## Ljung-Box Test R^2 Q(10) 6.623651 0.7604324
## Ljung-Box Test R^2 Q(15) 8.748086 0.8903165
## Ljung-Box Test R^2 Q(20) 10.53113 0.9574879
## LM Arch Test R TR^2 5.773293 0.927076
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -3.745151 -3.684598 -3.745667 -3.720939
Vhodný model pre naše dáta môžeme zvoliť napríklad podľa informačných kritérií, ku ktorým vieme pristúpiť nasledovne.
vynosy.30@fit$ics # ARCH(3)
## AIC BIC SIC HQIC
## -3.743054 -3.682501 -3.743570 -3.718842
vynosy.11@fit$ics # GARCH(1,1)
## AIC BIC SIC HQIC
## -3.736538 -3.688095 -3.736869 -3.717168
vynosy.12@fit$ics # GARCH(1,2)
## AIC BIC SIC HQIC
## -3.745151 -3.684598 -3.745667 -3.720939
Vidíme, že spomedzi troch modelov, ktoré boli pre naše dáta vyhovujúce, má najmenšiu hodnotu BIC práve GARCH(1,1). Ten ďalej použijeme aj na predikcie z modelu.
Predtým však overme stacionaritu nájdeného modelu.
GARCH(p,q) model je stacionárny, ak platí \[ (\alpha_1 + \alpha_2 + \dots + \alpha_p) + ( \beta_1 + \beta_2 + \dots + \beta_q) <1 \] V našom prípade sú odhadnuté koeficienty nasledovné
vynosy.11@fit$coef
## mu omega alpha1 beta1
## 0.0082960129 0.0002865938 0.2274863335 0.5938809441
Teda
sum(vynosy.11@fit$coef[3:4]) < 1
## [1] TRUE
Odhadnutý model je stacionárny.
Podobne ako pri Holt-Wintersovej metóde, predikcie budeme robiť
pomocou funkcie predict
.
vynosy.pred11 <- predict(vynosy.11, n.ahead = 50, plot = TRUE)
Keďže dáta modelujeme ako konštantu plus biely šum, predikcie do budúcnosti bude tiež iba samotná konštanta. Na druhú stranu, môžeme vidieť, že podmienená disperzia postupne narastá a konverguje k nepodmienenej disperzii časového radu. Prejaví sa to tak, že sa nám postupne rozšírujú intervaly spoľahlivosti.
plot(ts(vynosy.pred11[3]))
Zhrnutie:
Ak sa pozeráme len na korelácie, konštanta plus biely šum vyzerá ako dobrý model pre tieto dáta.
Ak sa však pozrieme na druhé mocniny rezíduí z takéhoto modelu, mali by byť nekorelované, ale nie sú. To nasvedčuje tomu, že disperzia nie je konštantná.
Ak modelujeme disperziu ARCH procesmi, potrebujeme vyšší počet parametrov na to, aby mal model dobré rezíduá, ale tieto parametre nie sú signifikantné.
Ak modelujeme disperziu GARCH procesmi, stačí nám menej parametrov a výjdu signifikantné.
Samozrejme, nie je to univerzálny postup a pri iných dátach môže byť vhodnejší ARCH model bez GARCH členov, alebo pri vyhovujúcom modeli zostanú nejaké parametre nesignifikantné.
Pozrime sa teraz na dáta o výmenných kurzoch. Dáta sú za obdobie 3 rokov (od 13.12.2021 do 10.12.2024), voľne dostupné na stránke https://www.ecb.europa.eu/stats/policy_and_exchange_rates/euro_reference_exchange_rates/html/index.en.html
Stiahnuté dáta upravené na logaritmické výnosy si do R-ka môžete načítať príkazom
data <- read.csv("https://anna-hlubinova.github.io/eurofx_log_returns.csv",
header = TRUE, sep=";")
head(data)
## Date USD JPY BGN CYP CZK DKK EEK
## 1 13.12.2021 0.000443439 -0.000078000 0 NA 0.001418272 0.00e+00 NA
## 2 14.12.2021 0.002744944 0.002026185 0 NA -0.002641176 0.00e+00 NA
## 3 15.12.2021 -0.004164642 -0.001558240 0 NA -0.002489874 0.00e+00 NA
## 4 16.12.2021 0.006549276 0.008695032 0 NA 0.000435187 0.00e+00 NA
## 5 17.12.2021 -0.000529427 -0.008461141 0 NA -0.001068524 2.69e-05 NA
## 6 20.12.2021 -0.005043589 -0.002341373 0 NA -0.000831864 -1.34e-05 NA
## GBP HUF LTL LVL MTL PLN ROL RON
## 1 -0.002310675 0.004067432 NA NA NA 0.002122499 NA 0.000000000
## 2 0.002193511 -0.000599520 NA NA NA 0.001707721 NA 0.000000000
## 3 -0.004285931 0.006087472 NA NA NA -0.000604935 NA -0.000020200
## 4 -0.001707741 0.000406322 NA NA NA 0.000129660 NA 0.000101019
## 5 0.004387133 -0.004533062 NA NA NA 0.001252511 NA 0.000060600
## 6 -0.000035200 -0.002396972 NA NA NA 0.000021600 NA -0.000303070
## SEK SIT SKK CHF ISK NOK HRK
## 1 -0.001494338 NA NA -0.000575760 0.001352265 0.002591990 -0.000531703
## 2 0.004875205 NA NA -0.001344732 -0.004062294 0.008595016 0.000039900
## 3 -0.002385881 NA NA 0.001152738 -0.004078864 -0.001924270 -0.000066500
## 4 -0.000975467 NA NA 0.003928525 0.000000000 -0.008049676 -0.000492090
## 5 0.002806773 NA NA -0.004600794 -0.004095569 0.003198176 0.000691517
## 6 0.002498021 NA NA -0.000096100 0.002732242 0.001296036 -0.000970919
## RUB TRL TRY AUD BRL CAD
## 1 0.000258413 NA 0.022790047 0.001901141 0.004234304 0.003341691
## 2 0.004487416 NA 0.009714382 0.004485021 0.008399504 0.007064230
## 3 -0.001353139 NA 0.021483628 -0.006513033 0.009853527 0.000413993
## 4 0.002430731 NA 0.059835792 -0.003113387 0.000372868 -0.000828157
## 5 0.004786452 NA 0.083601123 0.006406822 0.004185209 0.002964395
## 6 -0.001143024 NA 0.047400398 0.002526051 -0.006658052 0.002886997
## CNY HKD IDR ILS INR KRW
## 1 -0.000515353 0.000557182 -0.002914870 0.003115579 0.001311990 0.001688587
## 2 0.002490731 0.002951294 0.001495581 0.005748287 0.005329797 0.003525572
## 3 -0.003689604 -0.004157018 -0.000863800 0.003059839 0.001809050 -0.002176752
## 4 0.006839344 0.006477875 0.006174856 -0.002889602 0.004245084 0.004468062
## 5 0.000706288 -0.000497681 0.000175198 0.006983051 -0.003902257 0.001028262
## 6 -0.004884964 -0.005047336 -0.001511401 0.003178490 -0.006545590 0.000476517
## MXN MYR NZD PHP SGD THB
## 1 -0.003139267 0.003993028 0.000961365 -0.000264371 0.001557633 -0.008814594
## 2 0.010503211 0.003810007 0.004314743 0.002658148 0.003172653 0.002761259
## 3 0.004266886 -0.004649320 -0.002035198 -0.005429431 -0.004145889 -0.002681501
## 4 -0.004826244 0.000986121 -0.004744896 0.001889542 0.003240444 0.008049399
## 5 -0.009438359 0.002555189 0.007676660 -0.001059116 0.000258782 -0.005633073
## 6 -0.004537336 -0.003394379 0.002386921 -0.006467239 -0.002850297 0.004841600
## ZAR
## 1 -0.003998806
## 2 0.008450840
## 3 0.004168566
## 4 -0.004709494
## 5 -0.010216487
## 6 -0.005643965
Budeme modelovať výmenný kurz PLN/EUR, teda štrnásty stĺpec.
vynosy <- as.ts(data[,14])
plot(vynosy)
Otestujeme prítomnosť jednotkového koreňa
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0169828 -0.0023207 -0.0002984 0.0017972 0.0232015
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.94833 0.06203 -15.288 < 2e-16 ***
## z.diff.lag1 0.03217 0.04902 0.656 0.51188
## z.diff.lag2 -0.10189 0.03607 -2.824 0.00486 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003986 on 761 degrees of freedom
## Multiple R-squared: 0.4824, Adjusted R-squared: 0.4804
## F-statistic: 236.4 on 3 and 761 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -15.2875
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
Test zamieta hypotézu o jednotkovom koreni, teda dáta nebude potrebné diferencovať.
Ako vidno z vykreslenia logaritmických výnosov, dáta pozostávajú z viacerých období väčšej a menšej volatility.
Zobrazte ACF a PACF pre naše dáta
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0.06 -0.12 0.08 0.09 0.01 -0.06 0.02 -0.02 -0.03 -0.06 -0.02 -0.04 0.03
## PACF 0.06 -0.13 0.10 0.07 0.02 -0.05 0.02 -0.05 -0.01 -0.06 -0.01 -0.05 0.05
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.04 -0.02 0.00 0 0.01 -0.06 -0.02 -0.03 -0.04 -0.08 0.04 0.06
## PACF -0.05 0.01 -0.02 0 0.00 -0.06 -0.02 -0.04 -0.04 -0.07 0.05 0.04
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF -0.04 -0.06 0.00 0.14 0.00 -0.06 -0.03 0.02 0.01 0.00 0.03 -0.09
## PACF -0.02 -0.04 -0.02 0.12 -0.01 -0.04 -0.06 -0.02 0.00 0.02 0.03 -0.11
## [,38]
## ACF -0.01
## PACF 0.00
Skúsme dáta modelovať ako biely šum.
## [1] "initial value -5.515437 "
## [2] "iter 1 value -5.515437"
## [3] "final value -5.515437 "
## [4] "converged"
## [5] "initial value -5.515437 "
## [6] "iter 1 value -5.515437"
## [7] "final value -5.515437 "
## [8] "converged"
## [9] "<><><><><><><><><><><><><><>"
## [10] " "
## [11] "Coefficients: "
## [12] " Estimate SE t.value p.value"
## [13] "xmean -1e-04 1e-04 -0.6747 0.5001"
## [14] ""
## [15] "sigma^2 estimated as 1.619393e-05 on 769 degrees of freedom "
## [16] " "
## [17] "AIC = -8.187802 AICc = -8.187795 BIC = -8.175734 "
## [18] " "
V tomto prípade sa samotný časový rad nedá považovať za biely šum posunutý o konštantu. V ACF aj PACF pozorujeme autokorelácie mimo intervalov spoľahlivosti. Rovnako nevychádzajú ani Ljung-Boxové testy.
Poďme teda hľadať ďalej.
Nájdite pre dáta dobrý ARMA model. Z dobrých modelov vyberte ten s najnižším BIC. Zobrazte rezíduá a ich druhé mocniny.
Výstup pre kontrolu:
vynosy.ar3 <- sarima(vynosy, 3, 0 ,0) # BIC = -8.180398
## initial value -5.513825
## iter 2 value -5.528730
## iter 3 value -5.529124
## iter 4 value -5.529125
## iter 4 value -5.529125
## iter 4 value -5.529125
## final value -5.529125
## converged
## initial value -5.530717
## iter 2 value -5.530717
## iter 2 value -5.530717
## iter 2 value -5.530717
## final value -5.530717
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 0.0833 0.0358 2.3229 0.0204
## ar2 -0.1343 0.0356 -3.7702 0.0002
## ar3 0.1010 0.0358 2.8212 0.0049
## xmean -0.0001 0.0002 -0.6618 0.5083
##
## sigma^2 estimated as 1.570516e-05 on 766 degrees of freedom
##
## AIC = -8.21057 AICc = -8.210502 BIC = -8.180398
##
## [1] -0.01 0.01 0.00 0.05 0.04 -0.06 0.03 -0.04 -0.01 -0.07 0.00 -0.05
## [13] 0.04 -0.04 -0.01 0.00 -0.01 0.02 -0.07 -0.01 -0.05 -0.03 -0.07 0.04
## [25] 0.05 -0.05 -0.03 -0.01 0.14 0.00 -0.04 -0.04 0.01 0.02 -0.01 0.03
## [37] -0.09 0.00
## [1] 0.17 0.22 0.33 0.17 0.15 0.23 0.22 0.06 0.18 0.00 0.13 0.06
## [13] 0.01 0.11 -0.01 0.09 0.07 0.01 0.01 0.05 0.00 0.02 0.11 -0.01
## [25] -0.02 -0.01 0.04 0.02 0.00 -0.01 0.00 -0.04 0.00 -0.01 0.00 -0.01
## [37] -0.03 -0.04
Vidíme, že zatiaľ, čo autokorelačná funkcia pre rezíduá je v poriadku (výrazne vytŕčajú dve hodnoty, čo nemusí byť veľký problém), druhé mocniny rezíduí sú evidentne korelované.
Nájdite pre dáta vhodné GARCH modely.
Poznámka: Ak sa nejedná o čisto GARCH model, nie je dobré najskôr samostantne určiť rád členov ARMA procesu a potom dourčiť rád GARCH procesu, a to z dôvodu, že tieto dve rovnice sú navzájom prepojené, keďže v oboch by vystupovali nejaké rovnaké predošlé hodnot \(u_{t-s}\).
Ak vám vyšlo viacero dobrých modelov, vyberte najlepší na základe BIC.
Rovnako ako pri výnosoch akcií, spravíme aj teraz predikcie pre výnosy výmenných kurzov 12 dní dopredu.
Value at Risk (VaR) predstavuje mieru rizika, ktorá nám hovorí, aké najväčšie straty môžeme očakávať s pravdepodobnosťou \(1- \alpha\) v nasledujúcom období, kde \(\alpha\) je zvolená hladina pravdepodobnosti (napr. 0.05 alebo 0.01). Ide teda o kvantifikáciu rizika vo forme určitej maximálnej straty za predpokladu, že nenastane situácia v „chvoste“ rozdelenia.
Označme \(X\) hodnotu portfólia v nejakom čase. Potom \(\text{VaR}_\alpha\) je definovaná ako taká hodnota, pre ktorú platí \[ P(X \le \text{VaR}_\alpha) = \alpha. \]
Teda, \(\text{VaR}_\alpha\) predstavuje \(1-\alpha\) kvantil rozdelenia možných hodnôt X. Ak predpokladáme normálne rozdelenie pre naše dáta, tak to nie je nič iné ako \(1-\alpha\) kvantil normálneho rozdelenia.
Poznámka: VaR má svoje nedostatky, ako bolo spomenuté na prednáške. Jedným z nedostatkov je fakt, že VaR hovorí iba o maximálnej očakávanej strate na úrovni pravdepodobnosti \(1-\alpha\), ale neinformuje nás o priemernej veľkosti strát, ktoré sú väčšie ako \(\text{VaR}_\alpha\). Na tento účel slúži metrika Conditional VaR (CVaR) alebo Expected Shortfall.
Výpočet VaR si ukážeme na cenách akcie Apple, pre ktoré sme vyššie odhadli GARCH(1,1) model.
ceny <- AAPL.tyzden$AAPL.Adjusted
vynosy <- diff(log(ceny))
vynosy <- vynosy[-1]
VaR budeme počítať iteračne. Začneme po \(N=50\) hodnotách výnosov. Každý deň odhadneme GARCH(1,1) model pre naše dáta. Následne spravíme predikciu pre štandardnú odchýlku a pomocou nej zostrojíme VaR pre výnos na nasledujúci deň. Odhadovať budeme vždy z posledných 60 pozorovaní.
N <- 60
var.aaple <- rep(0, length(vynosy) - N + 1)
for(i in 1:length(var.aaple))
{
data.aapl <- vynosy[i:(i+N-1)]
garch.aapl <- garchFit(~garch(1,1), data = data.aapl, trace = FALSE)
aapl.pred <- predict(garch.aapl, n.ahead = 1, plot = FALSE)
var.aaple[i] <- qnorm(0.05, mean = aapl.pred[1,1], sd = aapl.pred[1,3])
}
Vytvoríme z výnosov aj z VaR časové rady. VaR začína až od N+1 dňa.
vynosy <- ts(vynosy)
var.aaple <- ts(var.aaple, start = N + 1)
Vykreslíme do jedného grafu
ts.plot(vynosy, var.aaple, col=c("black", "red"), lty = c(1, 1))
legend("bottomright", legend = c("Výnosy AAPL", "95% VaR"), col = c("black", "red"),lty = c(1, 1), bty = "n")
Vypočítajte počet prekročení VaR:
## [1] 18
Percentuálne:
## [1] 7.258065
Vidíme, že frekvencia prekorčení je viac ako 5%, teda náš model nemusí byť najvhodnejší na výpočet Value at risk . Rovnako by ešte trebalo testovať, či sú jednotlivé prekročenia od seba nezávislé.
Kupiecov test (angl. Kupiec’s test), často známy aj ako proportion of failures test, je štatistická metóda používaná na overenie správnosti odhadu Value at Risk (VaR) modelu. Tento test overuje, či je frekvencia porušení VaR (t. j. prípadov, keď reálna strata prekročí odhadovaný VaR) konzistentná s očakávanou pravdepodobnosťou porušenia, ktorá je \(\alpha\) (signifikantná hladina, napr. 5 % pre 95 % VaR).
Definujme funkciu \(I(\alpha) = \sum_{t=0}^n I_t(\alpha)\), kde \(I_t(\alpha)\) je rovné 1 ak \(r_t < \text{VaR}_t(\alpha)\), inak je \(I_t(\alpha)\) rovné 0. Pozorovanú frekvenciu prekročení vieme potom spočítať ako \(\hat{\alpha}= \frac{I(\alpha)}{n}\). Ďalej predpokladáme, že prekročenia sú nezávislé z alternativného rozdelenia s parametrom \(\alpha\). Potom počet prekročení má binomické rozdelenie.
Myšlienka Kupiecovho testu: \[ H_0: \alpha = \hat{\alpha} \quad \text{vs.} \quad H_1: \alpha \ne \hat{\alpha}. \]
Ak platí nulová hypotéza, potom likelihood-ratio štatistika má hodnotu \[ \text{LR}_{pof} = 2 \ln \Big[ \Big( \frac{1- \hat{\alpha}}{1- \alpha} \Big) ^ {n-I(\alpha)} \Big( \frac{\hat{\alpha}}{\alpha}\Big) ^{I(\alpha)} \Big] \sim \chi_1^2 \]
Testujte pomocou Kupiecovho testu prekročenia Value at Risk:
## [1] 0.1252441
freq <- ...
celkovy.pocet <- ...
LRuc <- 2*log(... )
pchisq(..., df = ..., lower.tail = FALSE)
P-value testu vyšla výrazne viac ako 5%, a preto nulovú hypotézu nezamietame. Teda aj keď pozorovaná frekvencia vyšla vyššia ako 5%, je to pravdepodobne iba štatistická chyba a môžeme vyhlásiť, že frekvencia prekročení je zhodná s predefinovanou hladinou VaR.
Aj sériovo závislé prekročenia nám robia problém. Zhlukovanie prekročení môže spôsobiť, že dva dni po sebe budeme pozorovať straty, na ktoré nebudeme kapitálovo pripravení. Preto je dôležité, aby platila nezávislosť prekročení.
Majme proces prekročení \(I_t\), ktorý môžeme modelovať ako Markovský reťazec: \[ \pi_{ij} = P(I_t(\alpha) = j \mid I_{t-1}(\alpha) = i). \]
Nech \(N_{ij}\) je počet pozorovaní v stave \(j\), ak v predošlej perióde bolo v stave \(i\) a označme \(N_0 = N_{00} + N_{01}\) a \(N_1 = N_{11} + N_{10}.\)
Ďalej označme odhad pravdepodobností prechodu zo stavov ako \(\pi_{01}\) a \(\pi_{11}\), ktoré vieme spočítať nasledovne \[ \hat{\pi}_{01} = \frac{N_{01}}{N_{01} + N_{00}}, \quad \hat{\pi}_{11} = \frac{N_{11}}{N_{11} + N_{10}}. \] Hypotézu o nezávislosti môžeme sformulovať ako:
\[ H_0 : \pi_{01} = \pi_{11} = \pi \quad \text{vs} \quad H_1 : \pi_{01} \neq \pi_{11}. \]
Za platnosti nulovej hypotézy dostávame likelihood ratio štatistiku:
\[ LR_{ind} = -2 \ln \left( \frac{(1 - \pi)^{N_{0}} \pi^{N_{1}}}{(1 - \pi_{01})^{N_{00}} \pi_{01}^{N_{01}} (1 - \pi_{11})^{N_{10}} \pi_{11}^{N_{11}}} \right) \sim \chi^2_1, \]
kde \(\hat{\pi}\) môžeme odhadnúť ako:
\[ \hat{\pi} = \frac{N_{01} + N_{11}}{N_{00} + N_{10} + N_{01} + N_{11}}. \]
## [1] 0.7606419
prekrocenia <- (vynosy[(N + 1):length(vynosy)] < var.aaple[-length(var.aaple)])*1
# 2 x 0 za sebou
N00 <- sum( ? )
N01 <- sum( ? )
N10 <- sum( ? )
N11 <- sum( ? )
N0 <- ?
N1 <- ?
p01 <- ?
p11 <- ?
p <- ?
LRind <- -2*log( ? )
LRind
pchisq(?, df = ?, lower.tail = FALSE)
P-hodnota testu vyšla aj v tomto prípade vysoko nad 5 percent, teda nulovú hypotézu opäť nezamietame. Prekročenia sú sériovo nezávislé.