Агрегирование наблюдений при работе с временными рядами

Обычной задачей при работе с временными рядами является расчет различных сводных показателей для тех или иных интересующих исследователя интервалов. В этой статье мы рассмотрим различные способы подобных вычислений по данным, хранящимся в таблицах класса tsibble. Этот формат ранее был подробно описан в одном из предыдущих сообщений, а также в книге "Анализ временных рядов с помощью R".

Агрегирование по избранным наблюдениям

Предположим, что перед нами стоит задача вычислить некоторый сводный показатель (например, среднее значение) по некоторым интересующим нам наблюдениям временного ряда. Разберемся, как это можно сделать в случае с данными, хранящимися в таблицах класса tsibble.

Такие таблицы отличаются от обычных таблиц классов data.frame или tibble только наличием дополнительных переменных - обязательной индексной и необязательной группирующей, которые описывают контекст каждого наблюдения. Поэтому к tsibble-таблицам применимы многие функции из пакета dplyr, в частности filter().

В приведенном ниже примере мы загружаем таблицу bitcoin с данными по стоимости биткоина (подробнее об этих данных см. здесь), а затем выбираем наблюдения только за 1--3 января 2019 г.:

require(readr)
require(dplyr)
require(tsibble)

# Ссылка для загрузки файл с данными:
url <- "https://raw.githubusercontent.com/ranalytics/tsa-r/master/data/bitcoin_price.csv"

# Загрузка данных и преобразование в формат `tsibble`
bitcoin <- read_csv(url) %>% 
  as_tsibble(., key = NULL, index = ds)

# Выбор данных за 1-3 января 2019 г.:
bitcoin %>% 
  filter(ds >= "2019-01-01", ds <= "2019-01-03")

#> # A tsibble: 3 x 2 [1D]
#>       y ds        
#> 1 <dbl> <date>    
#> 2 3844. 2019-01-01
#> 3 3943. 2019-01-02
#> 4 3837. 2019-01-03

Однако filter-команда оказалась относительно громоздкой. Чтобы упростить сделать синтаксис таких команд, в пакет tsibble была специально добавлена функция filter_index(). Интересующий пользователя сегмент временного ряда задается в этой функции с помощью формулы вида start ~ end, где start и end - это текстовые значения, соответствующие начальной и конечной точкам сегмента (включительно). В зависимости от стоящей задачи, элементы start или end можно опустить, что даст один из следующих возможных вариантов:

• ~ end или . ~ end: сегмент с начала временного ряда до отметки end;
• start ~ end: сегмент с отметки start до отметки end;
• start ~ .: сегмент с отметки start до конца временного ряда.

Используя функцию filter_index(), приведенный выше пример можно переписать более компактно:

bitcoin %>% 
  filter_index("2019-01-01" ~ "2019-01-03")

При необходимости на функцию filer_index() можно подать также несколько формул, задающих интересующие пользователя сегменты:

bitcoin %>% 
  filter_index("2019-01-01" ~ "2019-01-03",
               "2019-02-01" ~ "2019-02-03")

Для расчета сводных показателей по выбранным наблюдениям далее логично было бы воспользоваться функцией dplyr::summarise(). Давайте посмотрим, что получится, если применить этот подход к расчету средней стоимости биткоина за первые три дня 2019 г.:

bitcoin %>% 
  filter_index("2019-01-01" ~ "2019-01-03") %>% 
  summarise(avg_y = mean(y))

#> # A tsibble: 3 x 2 [1D]
#>   ds         avg_y
#>   <date>     <dbl>
#> 1 2019-01-01 3844.
#> 2 2019-01-02 3943.
#> 3 2019-01-03 3837.

Как видим, желаемого результата достичь не удалось: по сути, мы получили те же исходные данные за первые три дня 2019 г., хотя название второго столбца и изменилось с y на avg_y. Это произошло потому, что агрегирование наблюдений в таблицах класса tsibble допускается только в пределах уникальных значений индексной переменной. В нашем случае такой переменной является ds и она задает дневную периодичность учета наблюдений. В итоге мы получили средние значения по каждому дню, которые, естественно, оказались равны исходным наблюдениям. Для того чтобы рассчитать среднее значение по всем дням, перед вызовом функции summarise() необходимо изменить класс таблицы с tsibble на data.frame или tibble:

bitcoin %>% 
  filter_index("2019-01-01" ~ "2019-01-03") %>%
  as_tibble() %>% 
  summarise(avg_y = mean(y))

#> # A tibble: 1 x 1
#>   avg_y
#>   <dbl>
#> 1 3875.

Если в данных есть одна или нескольких группирующих переменных, то перед вычислением сводного показателя по каждой группе необходимо также применить стандартную функцию group_by() из пакета dplyr.

Агрегирование по календарным периодам

Для вычисления сводных показателей за определенный календарный период (например, средние или суммарные значения за день, неделю, месяц, квартал и т.п.) хорошо подходит связка из двух функций - tsibble::index_by() и dplyr::summarise().

Функция index_by() группирует данные по заданному пользователем временным интервалам. Для создания этих интервалов можно воспользоваться такими функциями из пакета tsibble, как yearweek() (неделя), yearmonth() (месяц) и yearquarter() (квартал). Ниже приведены примеры использования этих функций для расчета средних значений и стандартных отклонений стоимости биткоина по данным из таблицы bitcoin:

# Недельная периодичность:
btc_w <- bitcoin %>% 
  index_by(yw = yearweek(ds)) %>% 
  summarise(avg_y = mean(y), sd_y = sd(y))
glimpse(btc_w)

#> Rows: 191
#> Columns: 3
#> $ yw    <week> 2015 W53, 2016 W01, 2016 W02, 2016 W03, 201~
#> $ avg_y <dbl> 432.5933, 443.0057, 411.5671, 395.8229, 383.~
#> $ sd_y  <dbl> 2.281060, 11.477000, 32.634473, 15.399856, 9~

# Месячная периодичность:
btc_m <- bitcoin %>% 
  index_by(ym = yearmonth(ds)) %>% 
  summarise(avg_y = mean(y), sd_y = sd(y))
glimpse(btc_m)

#> Rows: 44
#> Columns: 3
#> $ ym    <mth> 2016 Jan, 2016 Feb, 2016 Mar, 2016 Apr, 2016~
#> $ avg_y <dbl> 410.8448, 404.4079, 416.5265, 434.3383, 461.~
#> $ sd_y  <dbl> 28.620019, 24.898333, 6.212295, 14.661023, 2~

# Квартальная периодичность:
btc_q <- bitcoin %>% 
  index_by(ym = yearquarter(ds)) %>% 
  summarise(avg_y = mean(y), sd_y = sd(y))
glimpse(btc_q)

#> Rows: 15
#> Columns: 3
#> $ ym    <qtr> 2016 Q1, 2016 Q2, 2016 Q3, 2016 Q4, 2017 Q1,~
#> $ avg_y <dbl> 410.7290, 512.4926, 615.7023, 732.7214, 1034~
#> $ sd_y  <dbl> 22.43261, 102.08236, 36.47290, 90.00748, 124~

При вычислении агрегированных показателей для нескольких временных рядов необходимо воспользоваться либо обычной функцией dplyr::group_by(), либо функцией tsibble::group_by_key(). Результат окажется идентичным, но в случае с group_by() необходимо будет указать имя группирующей переменной, тогда как group_by_key() этого в явном виде не требует.

Помимо описанных выше yearweek(), yearmonth() и yearquarter() для агрегирования по календарным периодам можно использовать также ряд функций из популярного пакета lubridate. В частности, при работе с внутридневными данными удобны функции round_date(), floor_date() и ceiling_date(), округляющие временные отметки до заданного пользователем разрешения (например, 5 секунд, 10 минут, 1 час и т.п.). Функции as_date() и year() позволяют также агрегировать данные по дням и годам соответственно. Вот пример:

# Расчет среднегодовой стоимости биткоина:
bitcoin %>%  
  index_by(year = year(ds)) %>% 
  summarise(avg_y = mean(y))

#> # A tsibble: 4 x 2 [1Y]
#>    year avg_y
#>   <dbl> <dbl>
#> 1  2016  568.
#> 2  2017 4006.
#> 3  2018 7572.
#> 4  2019 6771.

Агрегирование в пределах скользящего окна

Агрегирование в пределах скользящих блоков ("окон") применяется для решения широкого круга задач, связанных с анализом временных рядов. В качестве примеров можно привести вычисление скользящих средних, сумм и накопленных сумм, различных рыночных индикаторов и индексов, а также извлечение признаков для предсказательных моделей. Ранее подобные вычисления можно было выполнять с помощью нескольких функций из пакета tsibble. Однако с недавних пор эти функции были исключены из tsibble в пользу специализированного пакета slider. Рассмотрим основные особенности работы с ним.

Функция slide() и ее разновидности

Функция slide() - ключевая в пакете slider. Ее основными аргументами являются .x - вектор, над которым необходимо выполнить желаемые вычисления, и .f - существующая функция (например, .f = sum) или формула (например, .f = ~sum(.x)), задающие логику вычислений. Важно подчеркнуть, что таблицы с данными рассматриваются в пакете slider как векторы, элементы которых содержат строки таблиц. Поэтому аргумент .x может принимать также таблицы с данными. Как и при работе с популярным пакетом purrr, поданная на аргумент .f формула создает анонимную функцию (анонимными называют функции, которые объявляются в месте использования и не имеют собственного идентификатора для обращения к ним).

Если указать только аргументы .x и .f, то полученный результат окажется идентичным результату работы функции purrr::map():

require(slider)

slide(.x = 1:4, .f = ~.x) # то же, что и purrr:::map(.x = 1:4, .f = ~.x)

#> [[1]]
#> [1] 1
#> 
#> [[2]]
#> [1] 2
#> 
#> [[3]]
#> [1] 3
#> 
#> [[4]]
#> [1] 4

slide(.x = tibble(group = c("a", "a", "b", "b"), val = 1:4), .f = ~.x)

#> [[1]]
#> # A tibble: 1 x 2
#>   group   val
#>   <chr> <int>
#> 1 a         1
#> 
#> [[2]]
#> # A tibble: 1 x 2
#>   group   val
#>   <chr> <int>
#> 1 a         2
#> 
#> [[3]]
#> # A tibble: 1 x 2
#>   group   val
#>   <chr> <int>
#> 1 b         3
#> 
#> [[4]]
#> # A tibble: 1 x 2
#>   group   val
#>   <chr> <int>
#> 1 b         4

Однако у функции slide() есть и другие аргументы для организации вычислений в пределах скользящих окон:

• .before - количество наблюдений, предшествующих текущему (по умолчанию равно 0);
• .after - количество наблюдений, следующих за текущим (по умолчанию равно 0);
• .step - шаг, на который необходимо сдвинуть окно для выполнения следующей порции вычислений (по умолчанию равен 1).

Например, для формирования окна, охватывающего текущее наблюдение и два предшествующих ему наблюдения, необходимо присвоить аргументу .before значение 2:

slide(.x = 1:4, .f = ~.x, .before = 2)

#> [[1]]
#> [1] 1
#> 
#> [[2]]
#> [1] 1 2
#> 
#> [[3]]
#> [1] 1 2 3
#> 
#> [[4]]
#> [1] 2 3 4

Заметьте, что первые два элемента в полученном списке содержат "неполные окна". Такое поведение принято в функции slide() по умолчанию, т.е. предполагается, что пользователь желает выполнить некоторые вычисления даже по неполным окнам. Если такое поведение нежелательно, то следует воспользоваться аргументом .complete со значением TRUE:

slide(.x = 1:4, .f = ~.x, .before = 2, .complete = TRUE)

#> [[1]]
#> NULL
#> 
#> [[2]]
#> NULL
#> 
#> [[3]]
#> [1] 1 2 3
#> 
#> [[4]]
#> [1] 2 3 4

Та же команда, но со сдвигом на две позиции (.step = 2), даст следующий результат:

slide(.x = 1:4, .f = ~.x, .before = 2, .step = 2, .complete = TRUE)

#> [[1]]
#> NULL
#> 
#> [[2]]
#> NULL
#> 
#> [[3]]
#> [1] 1 2 3
#> 
#> [[4]]
#> NULL

Помимо окон фиксированного размера, функция slide() может также выполнять вычисления в пределах расширяющихся окон. Такие окна охватывают текущее наблюдение и все предшествующие ему наблюдения, что, в частности, позволяет вычислять различные показатели с накоплением. Для выполнения подобных вычислений следует присвоить аргументу .before значение Inf (аргумент .after также может принимать значение Inf, что приводит к эффекту "сужающегося окна"):

slide(.x = 1:4, .f = ~sum(.x), .before = Inf)

#> [[1]]
#> [1] 1
#> 
#> [[2]]
#> [1] 3
#> 
#> [[3]]
#> [1] 6
#> 
#> [[4]]
#> [1] 10

Как видно из приведенных выше примеров, slide() всегда возвращает результат с тем же количеством элементов, сколько их было во входных данных (даже если некоторые результаты вычислений оказываются равными NULL или NA). Кроме того, slide() всегда возвращает список с результатами вычислений. Однако у этой функции есть и другие разновидности, возвращающие объекты определенного класса. Названия этих функций состоят из приставки slide_ и таких окончаний, как lgl (логические значения), int (целые числа), dbl (действительные числа), chr (строковые значения), dfr (таблицы, созданные путем объединения строк) и dfc (таблицы, созданные путем объединения столбцов). Так, чтобы получить результат в виде вектора с целыми числами, следует применить функцию slide_int():

slide_int(.x = 1:4, .f = ~sum(.x), .before = Inf)

#> [1]  1  3  6 10

Функция slide_index() и ее разновидности

Во ременных рядах помимо последовательности из нескольких чисел всегда имеется также индексная переменная, "привязывающая" отдельные наблюдения к соответствующим временным меткам. Рассмотренные в предыдущем подразделе функции не учитывают временные метки в явном виде, поскольку они основаны на предположении, что анализируемые наблюдения учтены с одинаковым временным интервалом и в них нет пропущенных значений. Однако эти условия выполняются далеко не всегда. Так, в данных по стоимости акций обычно есть множество пропущенных значений, приходящихся на выходные и праздничные дни, когда биржевые торги не ведутся.

При наличии пропущенных или нерегулярно учтенных наблюдений агрегирование с помощью функций семейства slide следует выполнять с осторожностью, поскольку результаты вычислений могут отличаться от интуитивно ожидаемых. В таких случаях более надежным будет использование функции slide_index() и ее разновидностей (slide_index_lgl(), slide_index_dbl(), slide_index_chr() и т.д.), которые учитывают регулярность учета наблюдений. Все функции семейства slide_index имеют те же аргументы, что и slide-функции, за исключением одного дополнительного аргумента - .i, на который подается индексная переменная временного ряда.

В качестве примера воспользуемся несколькими наблюдениями из временного ряда стоимости акций Amazon (подробнее см. здесь). Заметьте, что в этом фрагменте временного ряда отсутствуют данные за 9-е и 10-е января:

amzn <- tibble(
  ds = as.Date(c("2016-01-04", "2016-01-05",
                 "2016-01-06", "2016-01-07",
                 "2016-01-08", "2016-01-11")),
  price = c(636.99, 633.79, 632.65, 607.94, 607.05, 617.74)
) %>% 
  as_tsibble(index = ds)

amzn

#> # A tsibble: 6 x 2 [1D]
#>   ds         price
#>   <date>     <dbl>
#> 1 2016-01-04  637.
#> 2 2016-01-05  634.
#> 3 2016-01-06  633.
#> 4 2016-01-07  608.
#> 5 2016-01-08  607.
#> 6 2016-01-11  618.

Применим теперь функции slide_dbl() и slide_index_dbl() для расчета скользящего среднего в пределах окна шириной в три наблюдения и сравним полученные результаты:

amzn %>% 
  mutate(
    mv_avg_slide = slide_dbl(
      .x = price, .f = ~mean(.x), .before = 2,
    ),
    mv_avg_slide_idx = slide_index_dbl(
      .x = price, .i = ds, .f = ~mean(.x),
      .before = 2,
    )
  )

#> # A tsibble: 6 x 4 [1D]
#>   ds         price mv_avg_slide mv_avg_slide_idx
#>   <date>     <dbl>        <dbl>            <dbl>
#> 1 2016-01-04  637.         637.             637.
#> 2 2016-01-05  634.         635.             635.
#> 3 2016-01-06  633.         634.             634.
#> 4 2016-01-07  608.         625.             625.
#> 5 2016-01-08  607.         616.             616.
#> 6 2016-01-11  618.         611.             618.

Различие между этими двумя функциями видно в 6-й строке полученной таблицы. Так, функция slide_dbl() "проигнорировала" пропущенные данные и рассчитала скользящее среднее для 11-го января по наблюдениям за 7-е, 8-е и 11-е января: (618 + 607 + 608) / 3 = 611. В то же время функция slide_index_dbl() учла наличие пропущенных значений и рассчитала скользящее среднее только по значению за 11-е января (618 / 1 = 618).

Функция slide_period() и ее разновидности

Как было показано выше, slide_index-функции группируют отдельные наблюдения временного ряда по принципу "текущее значение индексной переменной `.i` плюс несколько значений перед ним или после него". Функция slide_period() и ее разновидности (slide_period_lgl(), slide_period_dbl(), slide_period_chr() и т.д.) работают несколько иначе: они сначала разбивают исходный временной ряд на заданные пользователем календарные периоды и уже потом выполняют необходимое агрегирование данных.

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

amzn_two_month <- tibble(
  ds = as.Date(c("2016-01-28", "2016-01-29", 
                 "2016-01-30", "2016-01-31",
                 "2016-02-01", "2016-02-02",
                 "2016-02-03", "2016-02-04")),
  price = c(635.35, 587.00, NA, NA, 574.81, 552.10, 531.07, 536.26)
) %>% 
  as_tsibble(index = ds)

amzn_two_month

# A tsibble: 8 x 2 [1D]
#> ds         price
#> <date>     <dbl>
#> 1 2016-01-28  635.
#> 2 2016-01-29  587 
#> 3 2016-01-30   NA 
#> 4 2016-01-31   NA 
#> 5 2016-02-01  575.
#> 6 2016-02-02  552.
#> 7 2016-02-03  531.
#> 8 2016-02-04  536.

Предположим, что мы хотим вычислить среднюю стоимость акций по имеющимся данным отдельно за январь и февраль. Для этого необходимо воспользоваться аргументом .period со значением "month" (другие значения, которые принимает этот аргумент, включают "year", "quarter", "week", "day", "hour", "minute", "second", "millisecond", а также "yweek" - неделя года, "mweek" - неделя месяца, "yday" - день года и "mday" - день месяца):

slide_period_dfr(
  .x = amzn_two_month,
  .i = amzn_two_month$ds,
  .period = "month",
  .f = ~tibble(
     date = max(.x$ds),
     avg_price = mean(.x$price, na.rm = TRUE)
   )
)

#> # A tibble: 2 x 2
#>   date       avg_price
#>   <date>         <dbl>
#> 1 2016-01-31      611.
#> 2 2016-02-04      549.

Конечно, мы могли бы получить аналогичный результат, воспользовавшись функциями tsibble::index_by() и, например, tsibble::yearmonth():

amzn_two_month %>% 
  index_by(date = yearmonth(ds)) %>% 
  summarise(avg_price = mean(price, na.rm = TRUE))

#> # A tsibble: 2 x 2 [1M]
#>       date avg_price
#>      <mth>     <dbl>
#> 1 2016 Jan      611.
#> 2 2016 Feb      549.

Однако функции семейства slide_period предлагают гораздо больше возможностей. Например, применив аргумент .before = 1, мы можем рассчитать среднее значение стоимости акций в пределах окна, охватывающего данные за текущий календарный период и предшествующий ему период:

slide_period_dfr(
  .x = amzn_two_month,
  .i = amzn_two_month$ds,
  .period = "month",
  .f = ~tibble(
     date = max(.x$ds),
     avg_price = mean(.x$price, na.rm = TRUE)
   ),
  .before = 1
)

#> # A tibble: 2 x 2
#>   date       avg_price
#>   <date>         <dbl>
#> 1 2016-01-31      611.
#> 2 2016-02-04      569.

Еще один полезный аргумент у этих функций - .every. Он предназначен для группирования данных по заданному пользователем количеству календарных периодов. Сравните результат из предыдущего примера с результатом выполнения следующей команды, где .every = 2 (т.е. среднее значение рассчитывается по каждым двум месяцам):

slide_period_dfr(
  .x = amzn_two_month,
  .i = amzn_two_month$ds,
  .period = "month",
  .f = ~tibble(
     date = max(.x$ds),
     avg_price = mean(.x$price, na.rm = TRUE)
   ),
  .every = 2
)

#> # A tibble: 1 x 2
#>   date       avg_price
#>   <date>         <dbl>
#> 1 2016-02-04      569.

***

Эта статья представляет собой выдержку из второго, значительно дополненного издания книги "Анализ временных рядов с помощью R", работу над которым я сейчас завершаю. Как и ранее, электронная версия книги будет бесплатной и доступна онлайн для всех желающих. Кроме того, в московском издательстве ДМК Пресс выйдет бумажная версия книги - подписывайтесь на рассылку, чтобы не пропустить новости о скидках! Если сочтете уместным, вы можете поддержать меня в работе над новым изданием,  переведя любую сумму с помощью приведенной на этой странице формы оплаты.