Обнаружение структурных изменений во временных рядах

Практически любая промышленная система сегодня снабжена совокупностью сенсоров, датчиков и других измерительных устройств, позволяющих вести учет ее состояния (часто в реальном времени). Количество таких измерительных устройств и учитываемых с их помощью показателей может быть очень большим, что требует автоматизации процесса мониторинга и своевременного оповещения о нестандартных (например, предаварийных и аварийных) ситуациях. Подобные нестандартные ситуации обычно сопряжены с существенными изменениям во временных рядах, описывающих те или иные показатели. Ранее мы рассмотрели один из типов таких изменений - возникновение единичных точек-выбросов, или аномалий. Другим распространенным типом является сдвиг среднего уровня, который может быть как очень резким (рис. 1, вверху), так и постепенным (рис. 1, внизу). Существует несколько методов, позволяющих автоматически обнаруживать подобные структурные изменения во временных рядах. Одним из них является "E-Divisive with Medians" (EDM), разработанный исследователями из Корнеллского университета и компании Twitter (James et al. 2014) и реализованный в пакете для R  BreakoutDetection. Ниже приведено краткое описание основ работы с этим пакетом.

Рис. 1

Инсталляция пакета BreakoutDetection

Пакет BreakoutDetection можно установить только из репозитория Github (т.е. его нет в центральном хранилище CRAN и стандартная команда install.packages("BreakoutDetection") не сработает). Инсталляцию из Github проще всего выполнить с помощью функции install_github() из пакета devtools (который тоже нужно сначала установить, если его у вас еще нет):

install.packages("devtools")
devtools::install_github("twitter/BreakoutDetection")
library(BreakoutDetection)

Данные, используемые в примерах

В приведенных ниже примерах использованы рассмотренные ранее данные по стоимости биткоина (рис. 2). Этот небольшой набор данных можно найти в Github-репозитории ranalytics/breakout_detection, где хранится также весь код, представленный в этом сообщении.

Рис. 2

Метод "E-Divisive with Medians" (EDM)

Подробное описание теории, лежащей в основе EDM, приведено в статье James et al. (2014) и не является целью этого сообщения. Отметим лишь, что метод EDM основан на использовании т.н. E-статистики ("E" от "energy", т.е. "энергия"), которая отражает расстояние между наблюдениями из случайных выборок. Чем больше эта статистика, тем больше оснований считать, что сравниваемые выборки происходят из разных генеральных совокупностей. Во избежание влияния обычных на практике аномальных наблюдений, метод EDM оперирует не средними значениями, а медианами. Более того, этот метод непараметрический, т.е. он не исходит из предположения, что анализируемые данные происходят из каких-либо конкретных распределений случайных величин. Непараметрическая природа EDM, а также использование медианы в качестве меры центральной тенденции, делают этот метод особенно подходящим для обнаружения структурных изменений во временных рядах при решении практических задач.

Функция breakout() и ее аргументы

В пакете BreakoutDetection есть всего лишь одна функция - breakout(), которая имеет следующие аргументы:

• Z - временной ряд, который может быть представлен либо в виде числового вектора, либо в виде таблицы данных со столбцами timestamp ("временная отметка") и count ("количество").
• min.size - минимальное количество наблюдений, которое должно иметь место между точками излома тренда во временном ряду (по умолчанию равно 30).
• method - название применяемого к данным сценария по обнаружению структурных изменений. Принимает два возможных значения: "amoc" (от "at most one change", т.е. "не более одного изменения") и "multi" (несколько точек излома тренда).
• ... - дополнительные аргументы.

Набор дополнительных аргументов (...) зависит от значения параметра method. В обоих случаях (т.е. и "amoc", и "multi") функция примет:

• plot - логический аргумент. Значение TRUE приведет к построению графика временного ряда с отметками позиций обнаруженных точек излома тренда.
• xlab и ylab - подписи осей графика.
• title - заголовок графика.

В случае, когда method = "amoc", функция breakout() примет также:

• alpha - параметр, используемый в EDM для придания разных весов расстояниям между наблюдениями. Принимает любое значение из интервала (0, 2] (по умолчанию alpha = 2).
• exact - логический аргумент. Значение TRUE (принято по умолчанию) приведет к расчету точных значений медиан при нахождении точек излома (в случае с длинными рядами это может сопровождаться долгими вычислениями). FALSE вызовет расчет приблизительных значений медиан, что существенно сократит время вычислений.
• sig.lvl - желаемый уровень статистической значимости теста для обнаруженной точки излома тренда (по умолчанию 0.05).
• nperm - количество итераций для выполнения перестановочного теста. По умолчанию nperm = 0, т.е. перестановочный тест не выполняется.

Когда требуется обнаружить несколько точек излома тренда (method = "multi"), на функцию  breakout() можно подать также следующие аргументы, задающие уровень регуляризации в отношении количества и близости взаимного расположения этих точек:

• degree - параметр, определяющий близость расположения точек излома. Принимает значения 0, 1 (принято по умолчанию) и 2. Чем больше этот параметр, тем меньше точек излома будут расположены близко друг к другу.
• percent - задает минимальное относительное (%) увеличение внутреннего критерия качества EDM (здесь не рассматривается), необходимое для классификации того или иного наблюдения как точки излома тренда. Например, значение 0.25 соответствует 25%-ному увеличению критерия. Аргумент percent не имеет значения, принятого по умолчанию.
• beta - принимает любое числовое значение (по умолчанию равен 0.008) и используется для регуляризации если degree и/или percent не заданы. Чем больше этот параметр, тем меньше точек излома будет обнаружено алгоритмом.

Применим функцию breakout() с ее принятыми по умолчанию параметрами к данным по стоимости биткоина (рис. 3):

M0 <- breakout(dat$y, plot = TRUE,
               xlab = "Дата", ylab = "log стоимости ($)")
M0$plot

Рис. 3

Как видно из рис. 3, метод EDM обнаружил одну единственную точку излома тренда (отмечена вертикальной прерывистой линией), что обусловлено принятым по умолчанию значением method = "amoc" (см. выше). Для оценивания статистической значимости этого сдвига во временном ряду следует воспользоваться аргументом nperm (см. выше). Выполним такое оценивание с помощью перестановочного теста с 1000 итераций (внимание: выполнение этой команды займет какое-то время):

M0perm <- breakout(dat$y, nperm = 1000)

Объект M0perm представляет собой список со следующими элементами:

M0perm
$loc
[1] 561

$stat
[1] 107.6948

$time
[1] 36.84

$pval
[1] 0.000999001

где loc  - это положение точки излома (т.е. ее порядковый номер во временном ряду), stat - значение E-критерия, time - время выполнения команды (сек.), а pval - р-значение, полученное с помощью перестановочного теста. Как видим, обнаруженный сдвиг уровня временного ряда оказался статистически значимым (при уровне значимости 0.05).

Запросим теперь нахождение нескольких точек излома (method = "multi"; рис. 4):

M1 <- breakout(dat$y, method = "multi", plot = TRUE,
               xlab = "Дата", ylab = "log стоимости ($)")
M1$plot

Рис. 4

Рассмотрим влияние параметров регуляризации. Сначала отключим регуляризацию с помощью параметра degree (рис. 5):

M2 <- breakout(dat$y, method = "multi", plot = TRUE, degree = 0,
               xlab = "Дата", ylab = "log стоимости ($)")
M2$plot

Рис. 5

Как было отмечено выше, порядковые номера наблюдений, соответствующих обнаруженным точкам излома, хранятся в элементе loc получаемого объекта-списка:

M2$loc
 [1]   43   74  113  158  203  246  302  484  680  755  840
[12]  879  927  957  992 1048 1144 1201

Похожего эффекта (увеличение количества обнаруженных точек излома) можно добиться также с помощью низкого значения параметра percent (рис. 6):

M3 <- breakout(dat$y, method = "multi", plot = TRUE, percent = 0.05,
               xlab = "Дата", ylab = "log стоимости ($)")
M3$plot

Рис. 6

Наконец, снижение параметра beta тоже приведет к увеличению числа обнаруженных алгоритмом точек излома (рис. 7):

M4 <- breakout(dat$y, method = "multi", plot = TRUE, beta = 0.0001,
               xlab = "Дата", ylab = "log стоимости ($)")
M4$plot

Рис. 7

Заключение

Мы рассмотрели принципы работы с пакетом BreakoutDetection, в котором реализован один из наиболее эффективных методов обнаружения структурных изменений во временных рядах - "E-Divisive with Medians". Как отмечено выше, этот метод хорошо подходит для использования в автоматических системах мониторинга параметров технических систем (именно для этого он изначально и был разработан в компании Twitter). Другие примеры, где EDM может оказаться полезным, включают:

• анализ результатов, полученных в ходе A/B-тестирования (аргумент method функции breakout() со значением "amoc" позволит оценить влияние экспериментального фактора на сдвиг интересующей исследователя переменной-отклика);
• обнаружение существенных изменений в переменных, описывающих поведение людей (например, активность в социальных сетях) или животных;
• обнаружение сдвигов во временных рядах экономических показателей и т.п.

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