Введение

Сегодня обсудим две темы: моделирование хвостов и непараметрическое моделирование распределений. Снова будут финансовые штуки, немного статистики и пара новых R-овских фишек. А начнем мы с моделирования хвостов:

Распределение максимумов

Так как наша сверхзадача состоит в моделировании VaR и прочих мер риска, то ключевую роль для нас играют именно экстремальные значения, хвосты распределения. А особенно левый хвост. И мысль здесь следующая: нам не нужно уметь хорошо моделировать все распределение, нам нужно уметь хорошо моделировать только экстремальные события, особенно экстремально плохие события.

Выберем данные для нашего анализа:

library("quantmod")
getSymbols("LNKD", from="2011-01-01", to="2015-09-01")
## [1] "LNKD"
prices <- as.numeric(Ad(LNKD))
plot(prices, type="l", main="Акции LinkedIn")

ret <- diff(prices)/head(prices, -1)
hist(ret, breaks = 80)

Нам нужно будет делить все на блоки, поэтому выкинем последние точки, чтобы лучше делилось. Мы можем делить на блоки произвольной длины, я выберу 20 (в этом будет смысл). А значит я хочу, чтобы массив делился на блоки по 20 нацело.

length(ret)
## [1] 1078
k <- 20
length(ret)/k
## [1] 53.9
M <- floor(length(ret)/k)
M
## [1] 53
ret <- ret[1:(M*k)]

Максимумы убытков

Моделировать мы будем убытки, поэтому перейдем к ним:

loss <- -ret

Теперь пройдем по блокам и в каждом блоке найдем максимальный убыток:

maxs <- rep(0, M)
for (i in 1:M){
  maxs[i] <- max( loss[( (i-1)*k+1 ):( i*k )] )
}
maxs
##  [1] 0.08509791 0.06847321 0.17392732 0.07681260 0.05795244 0.06027495
##  [7] 0.04882692 0.05344365 0.03940885 0.03328412 0.02417045 0.05744122
## [13] 0.07524893 0.04776271 0.05401431 0.06829033 0.01495410 0.04153599
## [19] 0.04801898 0.01356498 0.01889913 0.02099055 0.02736225 0.05404824
## [25] 0.12932019 0.03993273 0.03155222 0.02590548 0.03402026 0.06104300
## [31] 0.09318604 0.04216867 0.04238345 0.05598233 0.06202730 0.04361530
## [37] 0.07820968 0.08367451 0.02904612 0.06251109 0.03550640 0.07614640
## [43] 0.03841779 0.08493059 0.04180848 0.03516329 0.02523469 0.02740041
## [49] 0.03070115 0.18609446 0.02096988 0.05733335 0.10517279
hist(maxs, breaks = 20)

В каком-то смысле это распределение VaR-ов, так как это распределения худших дней из блоков по 20. То есть мы взяли блок, выделили в нем максимальную потерю (то есть 5% худших данных) и построили их распределение.

median(maxs)
## [1] 0.04801898
quantile(loss, 1-0.05)
##        95% 
## 0.04426566

GEV

Большая наука тут в том, что какое бы ни было исходное распределение loss, его хвосты (maxs) распределены как GEV – Generalized extreme value distribution (Обобщенное распределение экстремальных значений). Распределение опирается на три параметра: среднее \(\mu\), дисперсию \(\sigma\) и параметр формы \(\xi\). PDF GEV’а выглядит не так плохо, как ОГР:

\[f(x) = \frac{1}{\sigma}t(x)^{\xi+1}e^{-t(x)}\]

и

\[ t(x) = \left(1+\left(\frac{x-\mu}{\sigma}\right)\xi\right)^{-1/\xi},\quad \xi \ne 0\] \[t(x) = e^{-(x-\mu)/\sigma},\quad \xi=0\]

Лучше один раз увидеть, не так ли? Новая фишка: curve рисует график функции от from до to в n Точках.

library(evd)
curve(dgev(x, loc = 0, scale = 1, shape=0), from = -4, to=4, n = 500, col=1, ylim=c(0,1))
curve(dgev(x, loc = 0, scale = 1, shape=0.5), from = -4, to=4, n = 500, col=2, add=TRUE)
curve(dgev(x, loc = 0, scale = 1, shape=1), from = -4, to=4, n = 500, col=3, add=TRUE)
curve(dgev(x, loc = 0, scale = 1, shape=-0.5), from = -4, to=4, n = 500, col=4, add=TRUE)
curve(dgev(x, loc = 0, scale = 1, shape=-1), from = -4, to=4, n = 500, col=6, add=TRUE)

Подгоняем под наше распределение

maxs.fit <- fgev(maxs)
maxs.fit
## 
## Call: fgev(x = maxs) 
## Deviance: -234.6923 
## 
## Estimates
##     loc    scale    shape  
## 0.03943  0.02023  0.17900  
## 
## Standard Errors
##      loc     scale     shape  
## 0.003179  0.002495  0.115179  
## 
## Optimization Information
##   Convergence: successful 
##   Function Evaluations: 80 
##   Gradient Evaluations: 9
plot(maxs.fit, which=1) # Можно не указывать which, а просто наживать enter

plot(maxs.fit, which=2)

plot(maxs.fit, which=3)

plot(maxs.fit, which=4)

Новые меры риска

Так как мы не работает непосредственно с нашим распределение доходностей, то и говорить прямо о VaR мы не можем. Но можем ввести две новые очень естественные меры, они обе есть на последнем графике. Но сначала извлечем параметры:

mu <- maxs.fit$estimate[1]
mu
##        loc 
## 0.03943258
sigma <- maxs.fit$estimate[2]
sigma
##      scale 
## 0.02023188
xi <- maxs.fit$estimate[3]
xi
##     shape 
## 0.1789951

Новая мера: уровень, который будет превзойдён в среднем 1 раз за n блоков по k наблюдений

n <- 12
R.kn <- mu + sigma/xi*((-log(1-1/n))^(-xi)-1)
R.kn
##       loc 
## 0.1013892

То есть раз в 12 блоков по 20 наблюдений мы будем ловить убыток в 10.1%.

Вторая мера - среднее время наступления события. Обратная предыдущей в смысловом плане:

BigLoss <- 0.2
n.k <- 1/(1-pgev(BigLoss, loc=mu, scale=sigma, shape=xi))
n.k
##    shape 
## 140.0888

То есть потерю в 20% можно ожидать раз в 140 блоков по 20 наблюдений (то есть где-то раз в 11 лет).

GPD

Вторая интересная идея – моделировать только хвост распределения. Опять же есть классная теорема, что хвост распределения должен хорошо фитироваться обобщенным распределение Парето (GPD, Generalized Pareto distribution). Позвольте мне хоть его формулу не записывать вот ссылка.

Будем моделировать худшие 20% доходностей.

alpha <- 0.2
GPD.fit <- fpot(loss, threshold = quantile(loss, 1-alpha), model = "gpd")
plot(GPD.fit, which = 1)

plot(GPD.fit, which = 2)

plot(GPD.fit, which = 3)

plot(GPD.fit, which = 4)

Можем найти Value-at-Risk:

scale <- GPD.fit$estimate[1]
xi <- GPD.fit$estimate[2]
VaR <- scale/xi*((alpha/0.05)^xi-1) + quantile(loss, 1-alpha)
-VaR
##       scale 
## -0.04429894
quantile(ret, 0.05)
##          5% 
## -0.04426566

Круто, работает.

Непараметрическое моделирование

Не будем навязывать природе (рынкам) наши мечты о формах распределения, а замоделируем все “как есть”. Наивная оценка плотности: просто будем скользить окном вдоль оси и считать, сколько точек попадают в наше окно. Но это будет просто гистограмма.

L <- 1e4 # В стольких точках оценим плотность
h <- 0.01 # Окно
points <- seq(-0.3, 0.3, length.out = L)
naive <- rep(0, L)
for (i in 1:L) naive[i] <- sum( (ret > points[i] - h/2) & (ret < points[i] + h/2) )
plot(points, naive, type="l")

Чтобы сгладить это дело можно использовать ядра. Тут можно все описать строго, можно рисовать картинки. Но лучше я расскажу это все у доски, так будет понятнее.

Разные формы ядер:

Kernels

Разные ядра не так сильно влияют на вид кривой:

hist(ret, breaks = 30)

plot(density(ret, kernel = "gaussian"))

plot(density(ret, kernel = "epanechnikov"))

plot(density(ret, kernel = "rectangular")) # Это и есть наша наивная оценка

plot(density(ret, kernel = "triangular"))

Куда более важный параметр – ширина окна:

hist(ret, breaks = 30)

plot(density(ret, 0.1, kernel = "epanechnikov"))

plot(density(ret, 0.01, kernel = "epanechnikov"))

plot(density(ret, 0.001, kernel = "epanechnikov"))

plot(density(ret, 0.0001, kernel = "epanechnikov"))

Если ничего не указать, то алгоритм сам выберет оптимальное окно:

hist(ret, breaks = 30)

plot(density(ret, kernel = "epanechnikov"))

density(ret)
## 
## Call:
##  density.default(x = ret)
## 
## Data: ret (1060 obs.);   Bandwidth 'bw' = 0.004786
## 
##        x                  y            
##  Min.   :-0.20045   Min.   : 0.000004  
##  1st Qu.:-0.09358   1st Qu.: 0.035824  
##  Median : 0.01329   Median : 0.176167  
##  Mean   : 0.01329   Mean   : 2.336988  
##  3rd Qu.: 0.12016   3rd Qu.: 1.695820  
##  Max.   : 0.22703   Max.   :19.751226

Теперь попросим оценить плотность в более широком диапазоне для будущей генерации случайных чисел. Сначала только узнаем этот диапазон:

range(ret)
## [1] -0.1860945  0.2126682
dens.fit <- density(ret, n=1e4, from = -0.4, to = 0.4) # Лучше брать n=1e5
x <- dens.fit$x
f <- dens.fit$y
f <- f/sum(f)  # Получим функцию плотности
plot(x, f, type="l")

F <- cumsum(f) # А так мы получим функцию распределения
plot(x, F, type="l")

Осталось только найти квантиль:

VaR <- x[which.min(abs(0.05 - F))] # Кто понимает, почему это квантиль?
VaR
## [1] -0.04548455

Задание

  1. Получить 1000+ значений котировок для какой-либо акции.
  2. Прогнать метод с GEV: получить уровень и среднее время наступления события. Приветствуется использования параметров отличных от тех, что в этой лекции.
  3. GPD: Прогнать циклом и построить кривую VaR. Проверить тестом Купика.
  4. Непараметрические оценки: выбрать ядро, выбрать окно. Прогнать кривую VaR, проверить Купиком.

В пункте 4 интереснее брать не очень большую обучающую выборку: 30-100 точек и смотреть на изменение теста Купика. В пункте 3 интересно поварьировать долю данных, которую мы хотим описать (в лекции это 20%).

Все непроверенные домашки не пропадают, они просто ждут своего часа :)