Данный ноутбук является конспектом по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2019). Автор ноутбука - вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
Вполне естественный вопрос: "Кто нафиг такой этот Максим, и что он забыл в Финансовом мире?!"
Спокойно. Сейчас всё расскажу.
library("ggplot2") # Пакет для красивых графиков
library("grid") # Пакет для субплотов
library("dplyr") # Куда же без пакета под таблички :)
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)
library("quantmod")
getSymbols("YNDX", from="2014-01-01", to="2020-01-01")
df = YNDX # Данные по акциям в таблицу df!
colnames(df) = c('open', 'high', 'low', 'close', 'volume', 'adjusted')
df$time = 1:nrow(df) # Колонка под время для удобства
df %>% head( ) # Все данные дневные
open high low close volume adjusted time 2014-01-02 43.31 43.34 42.36 42.63 1724200 42.63 1 2014-01-03 42.99 43.34 42.63 42.90 1130400 42.90 2 2014-01-06 42.89 43.19 42.61 42.91 1809900 42.91 3 2014-01-07 43.17 43.94 42.06 43.53 2821300 43.53 4 2014-01-08 43.77 44.24 42.22 42.25 3146300 42.25 5 2014-01-09 43.43 45.42 43.17 44.22 5685700 44.22 6
price <- as.numeric(df$close)
R = diff(price)/price[-length(price)]
В начале курса, когда мы обсуждали ЦПТ, мы говорили, что есть две зоны: среднестан и крайнестан.
Среднестан хорош тем, что по нему у нас огромная куча статистики. Случайные величины чаще всего вылезают на нас именно из него, в среднестане мы можем строить хорошие модели
Крайнестан очень непредсказуем и опасен. Данные из хвостов на нас выпрыгивают редко. Сложно набрать достаточное количество статистики, чтобы адекватно оценить с какой вероятностью произойдёт катастрофа. Наши оценки постоянно будут занижены. Для работы с крайнестаном нужны свои статистические методы.
qplot(R, bins=80)
Изучает хвосты распределений, помогает понять какие воздействия выскакивают на нас в экстремальных условиях
Используется в задачах, где важно смотреть на хвосты (риск-менеджмент, климатология, урбанистика и тп)
В ней есть куча теории вероятностей, своих тонкостей и разных теорем, на одну из них мы сейчас посмотрим чуть ближе
Так как наша сверхзадача состоит в моделировании $VaR$ и прочих мер риска, то ключевую роль для нас играют именно экстремальные значения, хвосты распределения. А особенно левый хвост. И мысль здесь следующая: нам не нужно уметь хорошо моделировать все распределение, нам нужно уметь хорошо моделировать только экстремальные события, особенно экстремально плохие события.
Нам нужно будет делить все на блоки, поэтому выкинем последние точки, чтобы лучше делилось. Мы можем делить на блоки произвольной длины, я выберу 20 (в этом будет смысл). А значит я хочу, чтобы массив делился на блоки по 20 нацело.
length(R)
k <- 20
length(R)/k
M <- floor(length(R)/k)
M
R <- R[1:(M*k)]
Моделировать мы будем убытки, поэтому перейдем к ним:
loss <- -R
Теперь пройдем по блокам и в каждом блоке найдем максимальный убыток:
maxs <- rep(0, M)
for (i in 1:M){
maxs[i] <- max(loss[((i-1)*k+1):(i*k)])
}
maxs
qplot(maxs, bins=20)
В каком-то смысле это распределение $VaR$-ов, так как это распределения худших дней из блоков по $20$. То есть мы взяли блок, выделили в нем максимальную потерю (то есть $5\%$ худших данных) и построили их распределение.
median(maxs)
quantile(loss, 1-0.05)
GEV (Generalized extreme value distribution) - обобщённое распределение экстремальных значений.
Распределение опирается на три параметра: среднее $\mu$, дисперсию $\sigma$ и параметр формы $\xi$. Функция распределения у GEV поприятнее, чем у ОГР:
В крайнестане GEV играет роль очень похожую на роль ЦПТ в среднестане.
ЦПТ говорила нам о том, как ведёт себя сумма $X_1 + \ldots + X_n$, при $n \to \infty$ (если ни одна из случаынйх величин не выделяется)
Теорема Фишера-Типпета-Гнеденко говорит нам о том, что максимум из $n$ случайных величин $\cdot \max(X_1, \ldots, X_n)$ при $n \to \infty$ будет вести себя в соответствии с GEV распределением
Важно, чтобы $X_1, \ldots, X_n \sim iid$
Когда мы говорили с вами о сходимостях случайных величин, мы даже доказали эту сходимость для частного случая. Мы посмотрели к чему будет сходится максимум из $n$ равномерных случайных величин.
# install.packages("evd")
library("evd")
maxs.fit <- fgev(maxs)
maxs.fit
Call: fgev(x = maxs) Deviance: -327.4042 Estimates loc scale shape 0.03855 0.01626 0.23726 Standard Errors loc scale shape 0.002285 0.001824 0.113916 Optimization Information Convergence: successful Function Evaluations: 51 Gradient Evaluations: 8
# which это тип картинки
plot(maxs.fit, which=2)
plot(maxs.fit, which=3)
# среднее значение
mu <- maxs.fit$estimate[1]
mu
# разброс
sigma <- maxs.fit$estimate[2]
sigma
# искажение
xi <- maxs.fit$estimate[3]
xi
Мера 1: уровень, который будет превзойдён в среднем $1$ раз за $n$ блоков по $k$ наблюдений (квантиль)
$$ q_{1-\frac{1}{k}} = F^{-1} \left(1-\frac{1}{k} \right) = \mu + \frac{\sigma}{\xi} \cdot \left( \left(-\ln \left(1 - \frac{1}{k} \right) \right)^{-\xi} - 1 \right) $$n <- 12
mu + sigma/xi*((-log(1-1/n))^(-xi)-1)
То есть раз в $12$ блоков по 20 наблюдений мы будем ловить убыток в $9.2\%$.
Мера 2: среднее время наступления события (математическое ожидание геометрической случайной величины).
$$ \frac{1}{1 - F(BigLoss)} $$BigLoss <- 0.2
1/(1-pgev(BigLoss, loc=mu, scale=sigma, shape=xi))
То есть потери в $20\%$ можно ожидать раз в $165$ блоков по $20$ наблюдений, то есть где-то раз в $13$ лет
20*165/252 # 252 торговых дня в году
Вторая интересная идея – моделировать только хвост распределения. Обычно основная проблема заключается в том, что мы не дооцениваем то, что происходит в хвосте из-за нехватки экстремальных данных. Можно попробовать моделировать хвосты с помощью обобщённого распределения Парето, GPD.
Функция распределения будет в этом случае выглядеть как
$$ F_{GPD}(x) = \begin{cases} 1 - \left(1+ \frac{\xi(x-\mu)}{\sigma}\right)^{-1/\xi} & \text{for }\xi \neq 0, \\ 1 - \exp \left(-\frac{x-\mu}{\sigma}\right) & \text{for }\xi = 0. \end{cases} $$Пусть случайная величина $X$ распределена с функцией распределения $F(x)$. Нас интересуют её значения, превышающие некоторый порог $u$, то есть значения из хвоста. Интересно было бы посмотреть как именно будут распределены такие значения. Чтобы найти функцию распределения, нужно воспользоваться формулой условной вероятности:
$$ F_u(x) = P(X - U \le x \mid X > U) = \frac{F(x + u) - F(u)}{1 - F(u)} $$Теорема Пикандса-Балкема-де Хаана говорит, что
$$ \sup_{x \in \mathbb{R}} \mid \hat F_u(x) - F_{GPD}(x) \mid \overset{\text{п.н.}}{\to} 0. $$Иначе говоря, эта теорема разрешает нам моделировать хвосты отдельно при большом числе наблюдений и высоком значении порога.
Будем моделировать худшие $10\%$ доходностей.
# Выставляем порог и находим пороговое значение
th = 0.1
u = quantile(loss, 1-th, names = FALSE)
GPD.fit = fpot(loss, threshold = u, model = "gpd")
GPD.fit
Call: fpot(x = loss, threshold = u, model = "gpd") Deviance: -833.0356 Threshold: 0.0308 Number Above: 136 Proportion Above: 0.1 Estimates scale shape 0.01435 0.18258 Standard Errors scale shape 0.001871 0.101231 Optimization Information Convergence: successful Function Evaluations: 40 Gradient Evaluations: 6
plot(GPD.fit, which = 2)
plot(GPD.fit, which = 3)
Дело за малым, найти VaR. Формулу для его поиска можно по-честному найти из функции распределения. Либо подсмотреть тут. Там же можно найти формулу для ES.
scale <- GPD.fit$estimate[1]
xi <- GPD.fit$estimate[2]
alpha = 0.05
VaR <- scale/xi*((th/alpha)^xi - 1) + u
-VaR
ES <- (VaR + sigma - xi*u)/(1 - xi)
-ES
В основе этой тетрадки (вплоть до копипасты) лежит курс моего друга, Ильи Езепова по рискам и всяким другим прикольным штукам в R. Для того, чтобы успешно его освоить, надо уметь хорошо проверять гипотезы, а также немного подучить эконометрику.
Вы смело можете попробовать ознакомиться с ним.
Пусть ссылки на все эти лекции лежат тут. Когда вы морально будете готовы вникнуть в них, читайте! :) Кстатит говоря, если хотите погрузится в финансовый мир не так резко, как мы это сделали сегодня, прочтите книгу Кванты и книгу Хулиганская экономика. Пацаны очень чётко в них всё расписывают. Такой научпоп сойдёт за курс введения в финансы.