Лекция 3 - Разговоры о волатильности
Илья Езепов
18 сентября 2015
Нюансы с прошлых лекций
Буквально пара моментов, которые не получилось обсудить в первых лекциях.
Считаем доходности
Я немного обманул в первой лекции касательно нахождения доходностей. Давайте сейчас раз и навсегда определимся, как считать ее в R быстро и эффективно. Пусть есть искусственный вектор цен a
a <- c(1, 1.5, 3, 2, 1, 2)
a## [1] 1.0 1.5 3.0 2.0 1.0 2.0Тогда есть несколько способов найти доходности:
a[2:length(a)]/a[1:(length(a) - 1)] - 1## [1] 0.5000000 1.0000000 -0.3333333 -0.5000000 1.0000000tail(a, -1)/head(a, -1) - 1## [1] 0.5000000 1.0000000 -0.3333333 -0.5000000 1.0000000diff(a)/head(a, -1)## [1] 0.5000000 1.0000000 -0.3333333 -0.5000000 1.0000000diff(a)/a[-length(a)]## [1] 0.5000000 1.0000000 -0.3333333 -0.5000000 1.0000000Случайна ли случайность?
Откуда компьютер берет случайные числа? Конечно, он их рассчитывает алгоритмически. А значит они не случайны. Собственно они и называются псевдослучайными. Но запуская много раз, я буду получать разные числа:
runif(1)## [1] 0.8401564runif(1)## [1] 0.9078401runif(1)## [1] 0.003306434Но бывают алгоритмы, которые работают на случайности (см. прошлую лекцию, Монте-Карло). Тогда как добиться их воспроизводимости? Хочется, чтобы при каждом запуске алгоритма я получал одинаковый результат. Это легко сделать, установив ядро генерации – отправную точку для ГСЧ (генератора случайных чисел). После установки ядра вся последовательность чисел будет воспроизведена.
runif(4)## [1] 0.7089046 0.8092093 0.2407360 0.5771166runif(4) # Все поменялось ...## [1] 0.9674804 0.4231410 0.1517501 0.3864418set.seed(42) # Ставим любое число, на свой вкус
runif(4)## [1] 0.9148060 0.9370754 0.2861395 0.8304476set.seed(42) # А теперь нет
runif(4)## [1] 0.9148060 0.9370754 0.2861395 0.8304476Введение
Сегодня займемся моделированием волатильности. Прежде всего обсудим, что такое волатильность. Узнаем, что такое стационарность и как ее проверить. И подгоним GARCH/ARCH/APARCH модели под наши данные. К делу: пакетов нужно будет много, лучше их сразу подключить. Для лекции возьмем котировки Старбакс с 2006 года.
#install.packages("fGarch")
#install.packages("FinTS")
#install.packages("tseries")
library("fGarch")
library("FinTS")
library("quantmod")
library("tseries")
getSymbols("SBUX", from="2006-01-01", to="2015-01-01")## [1] "SBUX"p <- as.numeric(Ad(SBUX))
plot(p, type="l", main="SBUX price", xlab="Price")
L <- length(p)
ret <- p[2:L]/p[1:(L-1)] - 1
hist(ret, breaks = 40)
Волатильность
Волатильность (в финансовом контексте) – мера разброса доходностей, дисперсия распределения доходностей (не совсем, но нас это устроит). Мы легко можем ее найти:
var(ret)## [1] 0.0004680171sd(ret)## [1] 0.0216337Но ведь можно вообразить ее функцией от времени:
vol <- c()
for (i in 31:length(ret)) {
vol <- c(vol, sd(ret[(i-30):i]))
}
plot(vol, type="l", main="Плато волатильности")
Смотрите как интересно! Совершенно явные плато, видимые зависимости. Скорее всего и ACF будет значима.
acf(vol)
acf(ret^2)
Квадраты доходностей иногда используют как меру волатильности. Напомню, что сами доходности не автокоррелированы.
acf(ret)
То есть дисперсия распределения меняется во времени, да еще и закономерным образом. А значит можно это все дело замоделировать.
Ну а раз повелась традиция вставлять забавные картинки, то вот картинка связанная с волатильностью. Фото сделано во время falsh crash в 2010:
Стационарность
Beat the market
Чуть-чуть отклонимся от повествования. Мы верим, что цены ведут себя случайным образом. Пусть тогда есть вот такой вектор цен:
p_art <- rnorm(100, mean=10)
plot(p_art, type="l", main = "Искусственная динамика цен")
Давайте придумаем торговый алгоритм: Если цена сегодня меньше 10, то покупать, если выше, то продавать:
trades <- rep(0, length(p_art))
for (i in 2:length(p_art)){
trades[i] <- ifelse(p_art[i] < 10, p_art[i+1] - p_art[i], -p_art[i+1] + p_art[i])
}
eq <- cumsum(trades)
plot(eq, main="Динамика капитала", type="l")
Вот и все. Оказывается, торговать на рынке предельно просто, я вас научил, расходимся.
На самом деле нет. Где я вас обманул?
Проблема была в том, что наша динамика цен была стационарна. Мы априори знали, что цена будет блуждать около 10, да еще и с постоянным стандартным отклонением. А как дело обстоит в реальности? А так, что на этой неделе котировки торгуются со средним 10 и отклонением 2, через неделю – со средним 7 и отклонением 5, через две – среднее 15 и отклонение 0.5. Свойства распределения меняются во времени, ряд не стационарен.
Стационарность в узком смысле: Ряд \(y\) стационарен, если распределение \(y_1, y_2,\dots,y_m\) совпадает с распределением \(y_{1+t}, y_{2+t},\dots,y_{m+t}\) для всяких \(m, t\).
Но это довольно сильное определение, так как мы требуем совпадения распределений. Что это вообще значит? Поэтому:
Стационарность в широком смысле: Среднее, дисперсия и ковариация ряда \(y\) не зависит от \(t\): \[E(y_t) = \mu < \infty\] \[V(y_t) = \gamma_0\] \[Cov(y_t, y_{t-k}) = \gamma_k \]
(Магнус, раздел 11.3).
Вот и выходит, что именно нестационарность реальных цен портит нам жизнь. Так может существуют стационарные цены, надо их только найти? попробуем. Но прежде нам нужен метод, который проверит ряд на стационарность.
Единичный корень
Рассмотрим ряд:
\[ y_t = \sum_{i=1}^p a_i y_t-i + \epsilon_t\]
Вспомним, что есть оператор лага \(L: Ly_t =y_{t-1}\). Тогда можем переписать модель:
\[ (1-\sum_{i=1}^p a_i L^i ) y_t = \epsilon_t\]
И назовем \(a(z) = (1-\sum_{i=1}^p a_i z^i )\) характеристическим полиномом. Если у этого полинома все корни по модулю больше 1, то ряд будет стационарным, если есть хотя бы один корень, равны1 1, то ряд будет не стационарен. Если есть корни меньше 1 по модулю, то процесс будет очень быстро разлетаться к бесконечности, и нам такие вообще не интересны.
Заметим еще, что если есть \(k\) единичных корней, то \(a(z) = (1-z)^k b(z)\), где у \(b(z)\) единичных корней нет. Тогда:
\[(1-\sum_{i=1}^p a_i L^i ) y_t = b(L) (1-L)^k y_t = b(L)\Delta^k y_t = \epsilon_t\]
А раз так, то \(b(L)\Delta^k y_t = \epsilon_t\) – процесс стационарный. То есть мало того, что мы научились определять стационарность процесса, так еще и можем свести любой с \(k\) единичными корнями к стационарному варианту.
Пока мы верим, что цены нестационарны, а доходности может и стационарны.
Тесты на едничный корень
Их очень много.
- Тест Дики-Фуллера: \(\Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \sum_{i=1}^{p-1}\delta_i \Delta y_{t-p+1} + \epsilon_t\). \(H_0: \gamma = 0\).
adf.test(p)##
## Augmented Dickey-Fuller Test
##
## data: p
## Dickey-Fuller = -1.1706, Lag order = 13, p-value = 0.9112
## alternative hypothesis: stationaryadf.test(ret)## Warning in adf.test(ret): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: ret
## Dickey-Fuller = -13.9651, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary- Тест Филлипса-Перрона (PP): \(y_t = \delta y_{t-1} + \epsilon_t\). \(H_0: |\delta| = 1\)
pp.test(p)##
## Phillips-Perron Unit Root Test
##
## data: p
## Dickey-Fuller Z(alpha) = -2.7127, Truncation lag parameter = 8,
## p-value = 0.9472
## alternative hypothesis: stationarypp.test(ret)## Warning in pp.test(ret): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: ret
## Dickey-Fuller Z(alpha) = -2287.525, Truncation lag parameter = 8,
## p-value = 0.01
## alternative hypothesis: stationary- Тест Квятковского — Филлипса — Шмидта — Шина (KPSS): еще более хитрые вещи, все равно никто их не запомнит. Важно, что в нем \(H_0\): стационарность.
kpss.test(p)## Warning in kpss.test(p): p-value smaller than printed p-value##
## KPSS Test for Level Stationarity
##
## data: p
## KPSS Level = 15.1415, Truncation lag parameter = 10, p-value =
## 0.01kpss.test(ret)## Warning in kpss.test(ret): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: ret
## KPSS Level = 0.3409, Truncation lag parameter = 10, p-value = 0.1ARCH и иже с ним
Авторегрессионная условная гетероскедастичность (ARCH — AutoRegressive Conditional Heteroscedasticity) – модель ряда, в котором дисперсия зависит от прошлых значений ряда.
Дальше будем говорить о ряде \(y_t\) и о его дисперсии \(\sigma^2_t\)
ARCH
Пусть дисперсия зависит от прошлых значений: \[\sigma^2_t = \alpha_0 + \sum_{i=1}^q \alpha_i y_{t-i}^2\]
Вот и вся наука, называется это ARCH(q), подгоняется обычным МНК.
GARCH
Обощим? Добавим зависимость не только от ряда, но и от прошлой дисперсии:
\[\sigma^2_t = \alpha_0 + \sum_{i=1}^p \alpha_i y_{t-i}^2 + \sum_{j=1}^p \beta_j \sigma_{t-j}^2\]
Звучит отлично. На практике любят GARCH(1,1). Можно ли использовать МНК для подгонки?
EGARCH
А одинаково ли на рынок влияют плохие и хорошие новости? И сто из этого следует.
Еще чуть-чуть обобщим:
\[\sigma^2_t = \alpha_0 + \sum_{i=1}^p \alpha_i (y_{t-i} - \gamma)^2 + \sum_{j=1}^p \beta_j \sigma_{t-j}^2\]
APGARCH
И совсем обобщим: кто вообще решил, что определяют все квадраты отклонения?
\[\sigma^2_\delta = \alpha_0 + \sum_{i=1}^p \alpha_i (|y_{t-i}| - \gamma y_{t-i})^\delta + \sum_{j=1}^p \beta_j \sigma_{t-j}^\delta\]
И моделей еще много-много-много… На сайте университета Дьюка есть файл С перечнем всех (почти) моеделей семейства. На 31 страницу за вычетом списка литературы. 31, Карл!
Хватит теории
В R много реализаций подгоник GARCH моделей. Два самых интересных пакета: fGarch и rugarch. Второй несколько быстрее и гибче, но с первым легче освоится. Фитировать будем вектор ret, не забываем.
Можно запустить вообще без параметров:
gfit <- garchFit(data=ret)##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 1)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 1
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: norm
## h.start: 2
## llh.start: 1
## Length of Series: 2264
## Recursion Init: mci
## Series Scale: 0.0216337
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -0.32157802 0.321578 0.0321578 TRUE
## omega 0.00000100 100.000000 0.1000000 TRUE
## alpha1 0.00000001 1.000000 0.1000000 TRUE
## gamma1 -0.99999999 1.000000 0.1000000 FALSE
## beta1 0.00000001 1.000000 0.8000000 TRUE
## delta 0.00000000 2.000000 2.0000000 FALSE
## skew 0.10000000 10.000000 1.0000000 FALSE
## shape 1.00000000 10.000000 4.0000000 FALSE
## Index List of Parameters to be Optimized:
## mu omega alpha1 beta1
## 1 2 3 5
## Persistence: 0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 2959.3019: 0.0321578 0.100000 0.100000 0.800000
## 1: 2943.4608: 0.0321584 0.0780867 0.102120 0.792078
## 2: 2933.1948: 0.0321594 0.0760240 0.121405 0.805166
## 3: 2926.4485: 0.0321620 0.0328486 0.139183 0.808272
## 4: 2916.3437: 0.0321733 0.0282469 0.161914 0.848914
## 5: 2899.2262: 0.0322486 0.0127222 0.126101 0.874619
## 6: 2897.7747: 0.0323009 0.0226665 0.0824884 0.888264
## 7: 2896.5686: 0.0323012 0.0233536 0.0857242 0.891574
## 8: 2894.8226: 0.0323029 0.0193562 0.0839212 0.893207
## 9: 2893.1218: 0.0323076 0.0178211 0.0809192 0.901936
## 10: 2891.3598: 0.0323176 0.0113205 0.0705227 0.916076
## 11: 2890.8833: 0.0323391 0.00938556 0.0609869 0.932052
## 12: 2890.4522: 0.0323792 0.00672124 0.0574325 0.935317
## 13: 2889.9587: 0.0324183 0.00828585 0.0610821 0.931482
## 14: 2889.9434: 0.0325292 0.00866676 0.0638896 0.928015
## 15: 2889.8931: 0.0327092 0.00853689 0.0641142 0.928766
## 16: 2889.8726: 0.0328895 0.00806458 0.0638783 0.929261
## 17: 2889.8592: 0.0330712 0.00813097 0.0638903 0.929407
## 18: 2889.7671: 0.0359793 0.00855422 0.0641735 0.928387
## 19: 2889.7211: 0.0388871 0.00752905 0.0641413 0.929836
## 20: 2889.6575: 0.0417954 0.00785071 0.0641363 0.929845
## 21: 2889.6468: 0.0421256 0.00801254 0.0654223 0.928381
## 22: 2889.6465: 0.0427963 0.00804256 0.0641212 0.928977
## 23: 2889.6313: 0.0431319 0.00825005 0.0645454 0.928593
## 24: 2889.6298: 0.0438222 0.00818825 0.0645883 0.928681
## 25: 2889.6298: 0.0440275 0.00828192 0.0649816 0.928202
## 26: 2889.6296: 0.0439246 0.00824527 0.0648440 0.928381
## 27: 2889.6296: 0.0439020 0.00823646 0.0648096 0.928425
## 28: 2889.6296: 0.0439019 0.00823650 0.0648098 0.928424
##
## Final Estimate of the Negative LLH:
## LLH: -5789.421 norm LLH: -2.557165
## mu omega alpha1 beta1
## 9.497606e-04 3.854822e-06 6.480983e-02 9.284245e-01
##
## R-optimhess Difference Approximated Hessian Matrix:
## mu omega alpha1 beta1
## mu -8350500.27 -1.096645e+08 -13804.31 -33309.24
## omega -109664456.84 -3.668464e+12 -639865148.70 -863385829.73
## alpha1 -13804.31 -6.398651e+08 -177402.87 -208366.58
## beta1 -33309.24 -8.633858e+08 -208366.58 -260806.59
## attr(,"time")
## Time difference of 0.03100586 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
## Time difference of 0.131026 secsИ узнать какие параметры есть: ?garchFit
Теперь попросим его говорить поменьше:
gfit <- garchFit(data=ret, trace=FALSE)
gfit##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(data = ret, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x000000000b45bc48>
## [data = ret]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## 9.4976e-04 3.8548e-06 6.4810e-02 9.2842e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 9.498e-04 3.466e-04 2.740 0.00614 **
## omega 3.855e-06 1.311e-06 2.941 0.00328 **
## alpha1 6.481e-02 1.130e-02 5.737 9.62e-09 ***
## beta1 9.284e-01 1.207e-02 76.902 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 5789.421 normalized: 2.557165
##
## Description:
## Fri Oct 02 14:04:46 2015 by user: EzepovISТеперь смысловые замечания. Попробуем более “глубокий” GARCH:
gfit <- garchFit(formula= ~garch(5,5),data=ret, trace=FALSE)## Warning in sqrt(diag(fit$cvar)): созданы NaNgfit##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(5, 5), data = ret, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(5, 5)
## <environment: 0x00000000353cf930>
## [data = ret]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 alpha3 alpha4
## 1.0239e-03 1.3838e-05 7.1360e-02 1.2324e-01 1.3677e-03 1.0000e-08
## alpha5 beta1 beta2 beta3 beta4 beta5
## 3.3063e-02 1.0000e-08 1.0000e-08 4.7697e-01 2.7231e-01 1.0000e-08
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 1.024e-03 3.424e-04 2.991 0.002785 **
## omega 1.384e-05 4.084e-06 3.389 0.000702 ***
## alpha1 7.136e-02 2.626e-02 2.718 0.006574 **
## alpha2 1.232e-01 4.140e-02 2.977 0.002912 **
## alpha3 1.368e-03 NA NA NA
## alpha4 1.000e-08 NA NA NA
## alpha5 3.306e-02 3.376e-02 0.979 0.327445
## beta1 1.000e-08 NA NA NA
## beta2 1.000e-08 NA NA NA
## beta3 4.770e-01 NA NA NA
## beta4 2.723e-01 NA NA NA
## beta5 1.000e-08 1.697e-01 0.000 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 5801.438 normalized: 2.562473
##
## Description:
## Fri Oct 02 14:04:47 2015 by user: EzepovISДобавим ARMA:
gfit <- garchFit(formula= ~arma(1,1) + garch(1,1), data=ret, trace=FALSE)
gfit##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(1, 1) + garch(1, 1), data = ret, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(1, 1) + garch(1, 1)
## <environment: 0x000000001caee498>
## [data = ret]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu ar1 ma1 omega alpha1
## 3.6115e-04 6.1093e-01 -6.4381e-01 3.8961e-06 6.5057e-02
## beta1
## 9.2811e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 3.612e-04 2.705e-04 1.335 0.18187
## ar1 6.109e-01 2.539e-01 2.406 0.01612 *
## ma1 -6.438e-01 2.517e-01 -2.558 0.01053 *
## omega 3.896e-06 1.321e-06 2.949 0.00319 **
## alpha1 6.506e-02 1.137e-02 5.722 1.05e-08 ***
## beta1 9.281e-01 1.216e-02 76.343 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 5791.026 normalized: 2.557874
##
## Description:
## Fri Oct 02 14:04:48 2015 by user: EzepovISА попробуем обобщить до APGARCH
gfit <- garchFit(formula= ~arma(1,1) + garch(1,1), data=ret, trace=FALSE, include.delta=TRUE)
gfit##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(1, 1) + garch(1, 1), data = ret, include.delta = TRUE,
## trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(1, 1) + garch(1, 1)
## <environment: 0x000000000bc33988>
## [data = ret]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu ar1 ma1 omega alpha1
## 0.00024127 0.72692388 -0.76838232 0.00060428 0.07179972
## beta1 delta
## 0.93545755 0.72561261
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 2.413e-04 8.346e-05 2.891 0.00384 **
## ar1 7.269e-01 4.225e-02 17.206 < 2e-16 ***
## ma1 -7.684e-01 4.178e-02 -18.392 < 2e-16 ***
## omega 6.043e-04 1.915e-04 3.155 0.00160 **
## alpha1 7.180e-02 9.687e-03 7.412 1.25e-13 ***
## beta1 9.355e-01 9.045e-03 103.423 < 2e-16 ***
## delta 7.256e-01 1.709e-01 4.247 2.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -3154.19 normalized: -1.393193
##
## Description:
## Fri Oct 02 14:04:49 2015 by user: EzepovISВот на этом этапе вам предоставляется свобода выбора модели для лальнейшей работы. Я возьму ARMA(1,1) + GARCH(1,1), как адекватную. Но призываю в ДЗ попробовать более грубые модели (GARCH(1,0)?) и более глупые (ARMA(5,5) + APGARCH(3,3)?).
И немного об оценке моделей:
plot(gfit, which=13)
plot(gfit, which=10)
NB: Есть еще одна важная штука – моделирование распределения остатков. Мы это не настраиваем, предполагая нормальность остатков. Это не очень хорошо, так как сохраняется автокорреляция остатков. Есть более общее распределение – GED, но так все будет дольше считаться и сложнее находить предсказание VaR. Поэтому остановимся на нормальном распределении, а желающие могут копнуть глубже и разобраться с параметром cond.dist функции garchFit.
Измеряем риски
Предсказываем полатильность
Имея объект fgit ничего не стоит предсказать волатильность на следующий день. Сначала просто картинку:
gfit <- garchFit(formula= ~arma(1,1) + garch(1,1), data=ret, trace=FALSE, include.delta=TRUE)
predict(gfit, 50, plot=TRUE, conf=0.95)
## meanForecast meanError standardDeviation lowerInterval upperInterval
## 1 0.0007525252 0.01224454 0.01224454 -0.02324634 0.02475139
## 2 0.0007882964 0.01262758 0.01261737 -0.02396130 0.02553790
## 3 0.0008142994 0.01301174 0.01299598 -0.02468824 0.02631684
## 4 0.0008332016 0.01339935 0.01338043 -0.02542905 0.02709545
## 5 0.0008469420 0.01379163 0.01377075 -0.02618416 0.02787804
## 6 0.0008569303 0.01418922 0.01416700 -0.02695343 0.02866729
## 7 0.0008641910 0.01459247 0.01456922 -0.02773652 0.02946490
## 8 0.0008694689 0.01500158 0.01497747 -0.02853310 0.03027203
## 9 0.0008733056 0.01541670 0.01539178 -0.02934287 0.03108948
## 10 0.0008760946 0.01583791 0.01581223 -0.03016563 0.03191782
## 11 0.0008781220 0.01626528 0.01623884 -0.03100124 0.03275748
## 12 0.0008795957 0.01669887 0.01667169 -0.03184960 0.03360879
## 13 0.0008806670 0.01713876 0.01711081 -0.03271068 0.03447202
## 14 0.0008814458 0.01758499 0.01755627 -0.03358449 0.03534738
## 15 0.0008820119 0.01803761 0.01800812 -0.03447106 0.03623508
## 16 0.0008824234 0.01849669 0.01846642 -0.03537043 0.03713528
## 17 0.0008827225 0.01896229 0.01893122 -0.03628269 0.03804813
## 18 0.0008829400 0.01943446 0.01940258 -0.03720790 0.03897378
## 19 0.0008830980 0.01991326 0.01988055 -0.03814617 0.03991237
## 20 0.0008832129 0.02039874 0.02036521 -0.03909759 0.04086402
## 21 0.0008832965 0.02089098 0.02085660 -0.04006228 0.04182887
## 22 0.0008833572 0.02139003 0.02135480 -0.04104033 0.04280704
## 23 0.0008834013 0.02189595 0.02185985 -0.04203187 0.04379867
## 24 0.0008834334 0.02240880 0.02237182 -0.04303701 0.04480388
## 25 0.0008834567 0.02292866 0.02289078 -0.04405588 0.04582280
## 26 0.0008834737 0.02345557 0.02341679 -0.04508860 0.04685555
## 27 0.0008834860 0.02398962 0.02394992 -0.04613530 0.04790227
## 28 0.0008834950 0.02453085 0.02449023 -0.04719609 0.04896308
## 29 0.0008835015 0.02507935 0.02503778 -0.04827112 0.05003813
## 30 0.0008835062 0.02563518 0.02559265 -0.04936052 0.05112753
## 31 0.0008835097 0.02619840 0.02615491 -0.05046441 0.05223143
## 32 0.0008835122 0.02676909 0.02672462 -0.05158294 0.05334996
## 33 0.0008835140 0.02734731 0.02730185 -0.05271624 0.05448326
## 34 0.0008835153 0.02793315 0.02788667 -0.05386445 0.05563148
## 35 0.0008835163 0.02852666 0.02847916 -0.05502771 0.05679474
## 36 0.0008835170 0.02912792 0.02907939 -0.05620616 0.05797320
## 37 0.0008835175 0.02973701 0.02968744 -0.05739996 0.05916699
## 38 0.0008835178 0.03035400 0.03030337 -0.05860924 0.06037627
## 39 0.0008835181 0.03097897 0.03092726 -0.05983415 0.06160119
## 40 0.0008835183 0.03161199 0.03155919 -0.06107484 0.06284188
## 41 0.0008835184 0.03225314 0.03219923 -0.06233147 0.06409851
## 42 0.0008835185 0.03290250 0.03284747 -0.06360419 0.06537123
## 43 0.0008835186 0.03356014 0.03350398 -0.06489315 0.06666018
## 44 0.0008835187 0.03422615 0.03416884 -0.06619850 0.06796554
## 45 0.0008835187 0.03490061 0.03484214 -0.06752042 0.06928745
## 46 0.0008835187 0.03558359 0.03552395 -0.06885905 0.07062608
## 47 0.0008835188 0.03627519 0.03621435 -0.07021456 0.07198159
## 48 0.0008835188 0.03697549 0.03691344 -0.07158711 0.07335415
## 49 0.0008835188 0.03768457 0.03762130 -0.07297687 0.07474391
## 50 0.0008835188 0.03840251 0.03833800 -0.07438401 0.07615105Нижняя линия – 95% доверительный интервал, это и есть VaR.
VaR <- as.numeric(predict(gfit, 1)[1] - 1.96*predict(gfit, 1)[3])
VaR## [1] -0.02324678Кривая VaR и Купик
Вам ведь уже знакомо?
N <- 500 # Для больших моделей нужно больше точек
test <- ret[(N+1):length(ret)]
VaR <- rep(0, length(test))
for (i in (N+1):length(ret)){
#cat("\r", i-N, "of", length(ret)-N)
train <- ret[(i-N):(i-1)]
gfit <- garchFit( formula= ~arma(1,1) + garch(1,1), data=train, trace=FALSE, include.delta=TRUE)
VaR[i-N] <- as.numeric(predict(gfit, 1)[1] - 1.96*predict(gfit, 1)[3])
}
plot(test, type="l", main = "Кривая VaR при моделировании ~ARMA(1,1) + GARCH(1,1)")
lines(VaR, col="red")
И тест Купика
L <- length(test)
K <- sum(test < VaR)
a0 <- K/L
a0 # Реальная доля пробитий## [1] 0.04478458a <- 0.05
S <- 2*log( (1-a0)^(L-K) * a0^K ) - 2*log( (1-a)^(L-K) * a^K )
pval <- 1 - pchisq(S, 1)
pval## [1] 0.3065977Мысли поп поводу кривой и результатов теста?
Домашнее задание
Вы уже догадались, не так ли?
- Взять котировки по вкусу (можно из первых ДЗ, лучше дневные, 500-700 торговых дней).
- Выбрать адекватную модель из семейства ARMA+APGARCH.
- построить кривую VaR и провести тест Купика.
- Сделать выводы.