Продолжаю свою серию заметок про трюки ggplot2, предыдущая заметка:

• Трюки ggplot2 - легенда 🤹🤹🤹
• Неопределенность и бизнес 🏌🏾⛳️ – прелюдия к этой заметке

Гистограммы 📊

Ранее я писал о важности использования гистограмм при анализе процессов в бизнесе и теперь пришло время поделиться несколькими трюками о том как можно сделать такие диаграммы более информативными. Такие трюки – являются итогом моего неустанного поиска более совершенных методов визуализации. Фрагменты этой статьи долгое время хранились в виде черновиков и теперь наступило время их собрать вместе и опубликовать 🖋

Баффет уже не торт 🍰

В качестве игрового примера данных я хочу вновь обратиться к котировкам американского фондового рынка. В этот раз я планирую исследовать феномен оракула из Омахи – Уоррена Баффета. Если кто-то не знает – это легендарный миллиардер, меценат, который сумел сделать состояние на инвестициях и заслужил место в тизере к этой статье. Уоррена Баффет стал примером для целого поколения инвесторов, создав ауру настоящего провидца, но ведущего скромный образ жизни и занимающейся благотворительностью. В общем очень милый, добрый дедушка, обладающий состоянием размером в годовой бюджет Польши или Израиля 👴🏻

Предметом моего небольшого исследования станет попытка выяснить чем котировки компании Баффета Berkshire Hathaway были похожи на остальной рынок и соответственно чем выделялись в лучшую или худшую стороны. Капитализацию Berkshire Hathaway я буду воспринимать как объективную характеристику успеха, полагая, что прибыль компании в той или иной степени трансформируется в капитализацию с поправкой на общерыночный уровень ликвидности. В качестве ориентира для сравнения буду использовать котировки индексов рынка США: S&P500 и Nasdaq 📈

Небольшая подготовка в виде загрузки библиотек и прочих рюшей по заведенной традиции:

library(thematic) # пакет для автоматической установки стилей графиков 
library(quantmod) # пакет для загрузки информации о биржевых котировках
library(tidyverse) # набор пакетов по принципу "все включено", в который включен ggplot2
library(memoise) # пакет кеширования результатов вывода функции
library(ggpp) # расширение для ggplot2
library(DT) # пакет создания интерактивных таблиц
library(patchwork) # создание композиции диаграмм
library(ggforce) # фокусировка на определенных местах диаграмм

# Активируем тему для блога
thematic_rmd(bg = "#1D1E20", accent = "cyan", fg = "grey90", 
             font = font_spec("Roboto"), sequential = firatheme::firaPalette(100), 
             qualitative = palette.colors(palette = "Tableau")) 

# Сохраняем палитру в отдельную переменную
my_pal <- palette.colors(palette = "Tableau") %>% unname() 

Далее для загрузки данных я буду использовать функцию, которую написал для статьи Трюки ggplot2 - легенда 🤹🤹🤹. Функция оказалась очень удачной и удобной, но я решил сделать небольшой апгрейд в виде добавления в параметры функции троеточия ..., что позволяет мне прокидывать параметры внешней функции get_symbols() во внутреннюю функцию getSymbols() без явного указания этих параметров. В данном случае мне нужно указать стартовую дату 01 января 1980 года, с которой я хотел бы получить данные. Такой вот удобный и полезный трюк 🤹

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

# Функция загрузки котировок
get_symbols <- function(sym_vec, ...){
  # Общие имена колонок, которые будут использованы для всех тикеров
  nms <- c("Time", "Open", "High", "Low", "Close", "Volume", "Adjusted")
  
  map(sym_vec, function(x)getSymbols(x, auto.assign = FALSE, ...)  %>%  
        as_tibble(rownames ="Time") %>% set_names(nms)) %>% 
    set_names(sym_vec)  %>%  
    bind_rows(.id = "Ticker")  %>%  
    mutate(Time = as.Date(Time))
}

# Команда кеширования результатов вывода функции 
get_symbols_cache <- memoize(get_symbols)

# Вектор тикеров вместе с вектором краивых имен
ticks <- c("^GSPC", "BRK-A", "^IXIC")
tick_names <- c("S&P500", "Berkshire", "Nasdaq") %>% set_names(ticks) 

# Запись котировок в переменную с прокидкой параметра from со значением "1980-01-01" внутрь getSymbols
tickers0 <- get_symbols_cache(ticks, from  = "1980-01-01") 

# Некоторые преобразования данных 
tickers1 <- tickers0 %>% 
  mutate(Month = str_sub(Time, 1L, 7L) %>% str_c("-01") %>% as.Date()) %>% 
  group_by(Ticker, Month) %>% 
  # агрегирование по месяцу
  summarise(across(-Time, ~mean(., na.rm = TRUE)),  .groups = "drop_last") %>% 
  # нормирование котировок ценных бумаг на 1
  mutate(across(-Month, ~scales::rescale(., to = c(0, 1))),
         # расчет изменения цен 
         Return  = Close - lag(Close, default = 0)) 

# Переменная для условного форматирования таблицы
brks <- quantile(tickers1[, 6], probs = seq(.05, .95, .05), na.rm = TRUE)

# Табличка с рюшечками и кружавчиками 
slice_sample(tickers1, n = 100) %>% # 100 случайных наблюдений
  datatable(style = 'bootstrap4',  extensions = 'Responsive', 
            options = list(pageLength = 10),
            caption = "Котировки технологических гигантов") %>% 
  formatRound(c(3:9, 8), digits = 4, mark = " ") %>% 
  formatStyle(6, backgroundColor = styleInterval(brks, firatheme::firaPalette(length(brks) + 1)))

Трюк № 1 – переиспользование элемента ggplot2 1️⃣

Лично я не очень люблю много тыкать в клавиатуру, предпочитаю переиспользовать фрагменты кода много раз. Более того, я искренне считаю, что если фрагмент кода нашел свое применение не один раз – это говорит о том, что этот фрагмент хорошо продуман и возможно его даже стоит оформить в виде функции публичного пакета. Далее я запишу функцию ggplot2::annotate() с определенным параметрами в переменную branding для того чтобы можно было эту переменную переиспользовать. Кроме того, я хадействую пакет patchwork, который позволяет создавать композии диаграмм с использованием оператора + в случае размещения диаграмм в горизонтальную линию и / в случае размещения диаграмм в стопку или веритикально.

# Переменная с фрагментом брендирования
branding <- annotate("text_npc", npcx = .5, npcy = .5, alpha = .9, size = 15, 
                     label = "InvestCookies.ru", color = "#1D1E20")

# Совместный график показателей
p1 <- ggplot(tickers1) + 
  branding + 
  geom_line(aes(Month, Close, col = Ticker)) + 
  scale_color_manual(values = my_pal, labels = tick_names)

# Показатели разбиты на фасеты
p2 <- ggplot(tickers1) + 
  branding + 
  geom_line(aes(Month, Close, col = Ticker)) + 
  scale_color_manual(values = my_pal, labels = tick_names) + 
  facet_grid(~Ticker, labeller = labeller(Ticker = tick_names))

p1/p2 + plot_layout(guides='collect') + plot_annotation(title = "Динамика котировок")

Видно, что трюк отлично сработал: брендинговая надпись появилась без особых проблем и вместо двух строчек кода для построения графика получилась только одна, что повысило читаемость кода.
В данном случае я сформировал двойной график: на верхнем можно легко заметить, что траектория всех трех показателей очень похожа. Лично я обратил внимание на кризис доткомов и стремительное падение выскотехнологического индекса Nasdaq в 2001 году после чего наблюдалось явное отставание Nasdaq вплоть до пандемии 2020 года. Также на нижних трех графиках отлично видно, что недавняя просадка Nasdaq вновь вывела этот индекс в аутсайдеры. Собственно тут же можно заметить, что вложения в S&P500 и Berkshire Hathaway принесли бы приблизительно одинаковую прибыль. На этом можно было бы остановиться и разбежаться, но в таком случае совсем не понятны причины феномена Баффета и следует копнуть глубже для чего я буду использовать диаграммы распределений.

Трюк № 2 – переиспользование нескольких элементов ggplot2 🖐

Переходим к композиции функций. К сожалению, простое присвоение работает только с одной функцией ggplot2 т.е. не получится сделать что-то подобное:

branding1 <-  try(annotate("text_npc", npcx = .5, npcy = .5, alpha = .9, size = 15, 
                           label = "InvestCookies.ru", color = "#1D1E20") +
                    scale_color_manual(values = my_pal, labels = tick_names))
## Error in try(annotate("text_npc", npcx = 0.5, npcy = 0.5, alpha = 0.9,  : 
##   Cannot add ggproto objects together. Did you forget to add this object to a ggplot object?

В данном примере функция try() используется просто для вывода ошибки и демонстрации того, что подход не работает. К счастью существует несколько способов выхода из положения. Первый достаточно прямолинейный и подразумевает написание функции, в которую первым аргументом передается объект ggplot2:

my_distr1 <-  function(gg, ...){
  gg + 
    annotate("text_npc", npcx = .5, npcy = .5, alpha = .9, size = 15, 
             label = "InvestCookies.ru", color = "#1D1E20") +
    geom_histogram(...) +
    scale_fill_manual(values = my_pal, labels = tick_names) + 
    scale_x_continuous(label = scales::label_percent()) 
}

ggplot(tickers1, aes(Return, fill = Ticker)) %>% 
  my_distr1(position = "dodge", alpha = .5, binwidth = .002)

Функция в целом работает, но требует указания операторов %>% или |> вместо привычного +, что несколько сбивает с толку. Есть второй метод, который позволяет использовать + для формирования каскада функций ggplot2. Дело в том, что к объекту ggplot2 можно обращаться как к списку через list():

my_distr2 <-  function(...){
  list(
    annotate("text_npc", npcx = .5, npcy = .5, alpha = .9, size = 15, 
             label = "InvestCookies.ru", color = "#1D1E20"),
    geom_histogram(...),
    scale_fill_manual(values = my_pal, labels = tick_names), 
    scale_x_continuous(label = scales::label_percent()) 
  )
}

ggplot(tickers1, aes(Return, fill = Ticker)) +
  my_distr2(position = "dodge", alpha = .5, binwidth = .002)

Просто перечисляем функции через запятую и получаем результат. В итоге имеет ровно такую же диаграмму, как и предыдущую, но уже с использованием привычного для ggplot2 операторы +.

Трюк № 3 – гистограмма и плотность распределения 🗑

Иногда очень полезно дополнять гистограмму сглаженной кривой распределения, которая может показать то, что не показывает гистограмма. При этом хотелось бы чтобы такая кривая накладывалась на гистограмму, но это затруднительно в силу того, что кривая распредления нормирована определенным образом в то время как гистограмма не нормирована вовсе.

Существует несколько способов аккуратного наложения двух диаграмм, но лично я пользуюсь единственным наиболее прозрачным вариантом, который управляется параметром binwidth т.е. буквально шириной корзины, если следовать аналогии, которую я приводил в статье Неопределенность и бизнес 🏌🏾⛳️.

Для построения графика плотности вероятности можно использовать конструкцию geom_density(aes(y = after_stat(density))), но она дает не вполне нужный эффект:

bw <- .004

filter(tickers1, Ticker == "BRK-A") %>% 
  ggplot(aes(Return)) + 
  geom_histogram(binwidth = bw) + 
  geom_density(aes(y = after_stat(density)), col = my_pal[3], fill = my_pal[3], alpha = .5)

Как можно заметить, кривая распределения не покрывает гистограмму и тем самым искажает представление о том, что происходит в данных.

Для того чтобы кривая распределения проходила по верхушкам столбиков гистограммы необходимо использовать конструкцию geom_density(aes(y = after_stat(density*bw))), где bw соответствует binwidth в гистограмме:

filter(tickers1, Ticker == "BRK-A") %>% 
  ggplot(aes(Return)) + 
  geom_histogram(binwidth = bw) + 
  geom_density(aes(y = after_stat(count*bw)), col = my_pal[3], fill = my_pal[3], alpha = .5)

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

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

• bw – для управления детализацией гистограммы и кривой распределения
• split – для включения/выключения гистограммы
• ... – для прокидки любых параметров внешней функции во внутренние

# Собственная функция ggplot2 
my_distr2 <-  function(bw, split = FALSE, ...){
  list(
    annotate("text_npc", npcx = .5, npcy = .5, alpha = .9, size = 15, 
             label = "InvestCookies.ru", color = "#1D1E20"),
    if(split) # правило активации визуализации гистограммы
       geom_histogram(binwidth = bw, alpha = .9, ..., fill = "grey50", col = "grey50"),
    if(split) # правило активации фасетной визуализации 
       facet_grid(Ticker ~ ., labeller = labeller(Ticker = tick_names), scales = "free"),
    geom_density(aes(y = after_stat(count*bw)), ...), # уже известная конструкция
    scale_fill_manual(values = my_pal, labels = tick_names), 
    scale_color_manual(values = my_pal, labels = tick_names), 
    scale_x_continuous(label = scales::label_percent()), 
    labs(x = "Изменения за месяц", y = "Наблюдения")
  )
}

ggplot(tickers1, aes(Return, fill = Ticker, col = Ticker)) +
  my_distr2(position = "dodge", alpha = .2, bw = .005, size = .1) + 
  # Функция пакета ggforce для создания эффекта зума определенной области 
  facet_zoom(xlim =c(-.12, -.05), ylim = c(0, 4), horizontal = FALSE, zoom.size = .5) +
  geom_segment(x = -.12, xend = -.05, y = 20, yend = 20, 
               arrow = arrow(ends = "both", length = unit(0.1, "inches")),
               size = .01, col = my_pal[3]) + 
  annotate(x = -.12, y = 90, geom = "label", label = "Зона\nчерных лебедей", hjust = 0, 
           col = "black", fill = my_pal[3], alpha = .5) 

Отлично! получилась хорошая функция, которая позволила визуализировать распределения в нужной детализации. Как можно заметить, распределения не являются нормальными, более того, после логарифмирования они также не станут нормальными, хотя логарифмировать в лоб такие значения нельзя т.к. они не являются положительно-определенными. В целом, аналитическая идентификация и анализ распределений – это отдельная большая тема, которую оставлю за рамками данной заметки.

Если посмотреть на график, то можно обнаружить:

• Шип в плотности распределения Berkshire Hathaway выше шипов индексов – это означает, что компания Баффета чаще торгуется около нуля чем индексы, т.е. является в каком то смысле более стабильной
• У индексов присутствует горб в диапазоне 1-2% но Berkshire Hathaway имеет больше значений в диапазоне 2-5% т.е. компания Баффета растет более нервно: скорей всего на инфоповодах
• Количество черных лебедей т.е. экстремально отрицательных значений в диапазоне от -12% до -5% приблизительно одинаково у всех показателей, а именно 3 шт. все показатели падали +/- одинаково. Именно тут я искал отличия для Berkshire Hathaway т.к. профессиональные инвестиционные стратегии включают опционы, способные существенно ограничить риски. К сожалению, в данном моменте я также не обнаружил какой-либо магии 🧙‍♂️

Далее, используя написанную функцию, можно построить одновременно гистограмму с наложением диаграммы кривой распределения, что позволит обогатить понимание динамики трех показателей:

ggplot(tickers1, aes(Return, fill = Ticker, col = Ticker)) +
  my_distr2(position = "dodge", alpha = .5, bw = .001, size = .1, split = TRUE) 

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

Для тех, кто не сталкивался с феноменом черных лебедей, но хочет об этом узнать немного больше я бы предложил немного погрузиться в художественное описание того, что происходит в моменты прилета черного лебедя. Обычно, за безжизненными бугорками в распределениях скрываются сломанные судьбы или наоборот возрождение фениксов из пепла.

Краткие итоги

• Для того чтобы явно не указывать аргументы внутренней функции во внешней можно использовать троеточие ..., что позволит лениво прокинуть любые аргументы внутрь 🦥
• Для того чтобы переиспользовать функцию ggplot2 достаточно ее присвоить в переменную 1️⃣
• Для того чтобы переиспользовать каскад функций, либо необходимо обернуть такой каскад в функцию с первым аргументом в виде объекта ggplot2, либо использовать внутри такой функции list() как более органичный способ для ggplot2, который позволяет использовать оператор +🖐
• Построение сглаженной кривой распределения, соответствующей гистограмме возможно с использованием конструкции geom_density(aes(y = after_stat(count*bindwidth))) 🗑
• Баффет уже не торт, да собственно тортом он никогда и не был: вложения в индекс S&P500 за более чем 40 лет принесли бы приблизительно такую же доходность как вложения в Berkshire Hathaway т.е. феномен оракула из Омахи остался неразгаданным 🤯

Простой способ узнать о новых публикациях – подписаться на Telegram-канал: