ПОИСК ЗАКОНОМЕРНОСТЕЙ В ЧИСЛОВЫХ РЯДАХ PYTHON

Кросс-валидация на временных рядах, подбор параметров

Перед тем, как построить модель, поговорим, наконец, о не ручной оценке параметров для моделей.

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

Небольшая загвоздка возникает только в кросс-валидации. Проблема состоит в том, что временной ряд имеет, как ни парадоксально, временную структуру, и случайно перемешивать в фолдах значения всего ряда без сохранения этой структуры нельзя, иначе в процессе потеряются все взаимосвязи наблюдений друг с другом. Поэтому придется использовать чуть более хитрый способ для оптимизации параметров, официального названия которому я так и не нашел, но на сайте CrossValidated, где можно найти ответы на всё, кроме главного вопроса Жизни, Вселенной и Всего Остального, предлагают название «cross-validation on a rolling basis», что не дословно можно перевести как кросс-валидация на скользящем окне.

Суть достаточно проста — начинаем обучать модель на небольшом отрезке временного ряда, от начала до некоторого

, делаем прогноз на

шагов вперед и считаем ошибку. Далее расширяем обучающую выборку до

значения и прогнозируем с

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

Код для кросс-валидации на временном ряду

Значение длины сезона 24*7 возникло не случайно — в исходном ряде отчетливо видна дневная сезонность, (отсюда 24), и недельная — по будням ниже, на выходных — выше, (отсюда 7), суммарно сезонных компонент получится 24*7.

В модели Хольта-Винтерса, как и в остальных моделях экспоненциального сглаживания, есть ограничение на величину сглаживающих параметров — каждый из них может принимать значения от 0 до 1, поэтому для минимизации функции потерь нужно выбирать алгоритм, поддерживающий ограничения на параметры, в данном случае — Truncated Newton conjugate gradient.

Out: (0.0066342670643441681, 0.0, 0.046765204289672901)

Передадим полученные оптимальные значения коэффициентов

и построим прогноз на 5 дней вперёд (128 часов)

Код для отрисовки графика

Судя по графику, модель неплохо описала исходный временной ряд, уловив недельную и дневную сезонность, и даже смогла поймать аномальные снижения, вышедшие за пределы доверительных интервалов. Если посмотреть на смоделированное отклонение, хорошо видно, что модель достаточно резко регирует на значительные изменения в структуре ряда, но при этом быстро возвращает дисперсию к обычным значениям, «забывая» прошлое. Такая особенность позволяет неплохо и без значительных затрат на подготовку-обучение модели настроить систему по детектированию аномалий даже в достаточно шумных рядах.

Экспоненциальное сглаживание, модель Хольта-Винтерса

А теперь посмотрим, что произойдёт, если вместо взвешивания последних

Здесь модельное значение представляет собой средневзвешенную между текущим истинным и предыдущим модельным значениями. Вес

называется сглаживающим фактором. Он определяет, как быстро мы будем «забывать» последнее доступное истинное наблюдение. Чем меньше

Экспоненциальность скрывается в рекурсивности функции — каждый раз мы умножаем

на предыдущее модельное значение, которое, в свою очередь, также содержало в себе

, и так до самого начала.

Двойное экспоненциальное сглаживание

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

В этом нам поможет разбиение ряда на две составляющие — уровень (level, intercept)

В результате получаем набор функций. Первая описывает уровень — он, как и прежде, зависит от текущего значения ряда, а второе слагаемое теперь разбивается на предыдущее значение уровня и тренда. Вторая отвечает за тренд — он зависит от изменения уровня на текущем шаге, и от предыдущего значения тренда. Здесь в роли веса в экспоненциальном сглаживании выступает коэффициент

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

Теперь настраивать пришлось уже два параметра —

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

Тройное экспоненциальное сглаживание a. Holt-Winters

Итак, успешно добрались до следующего варианта экспоненциального сглаживания, на сей раз тройного.

Идея этого метода заключается в добавлении еще одной, третьей, компоненты — сезонности. Соответственно, метод применим только в случае, если ряд этой сезонностью не обделён, что в нашем случае верно. Сезонная компонента в модели будет объяснять повторяющиеся колебания вокруг уровня и тренда, а характеризоваться она будет длиной сезона — периодом, после которого начинаются повторения колебаний. Для каждого наблюдения в сезоне формируется своя компонента, например, если длина сезона составляет 7 (например, недельная сезонность), то получим 7 сезонных компонент, по штуке на каждый из дней недели.

Получаем новую систему:

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

шагов вперёд, что не может не радовать.

Ниже приведен код для построения модели тройного экспоненциального сглаживания, также известного по фамилиям её создателей — Чарльза Хольта и его студента Питера Винтерса.
Дополнительно в модель включен метод Брутлага для построения доверительных интервалов:

— длина сезона,

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

Код для модели Хольта-Винтерса

Домашнее задание №9

В демо-версии домашнего задания вы будете предсказывать просмотры wiki-страницы «Machine Learning». Веб-форма для ответов, там же найдете и решение.

Актуальные и обновляемые версии демо-заданий – на английском на сайте курса, вот первое задание. Также по подписке на Patreon («Bonus Assignments» tier) доступны расширенные домашние задания по каждой теме (только на англ.).

Линейная регрессия vs XGBoost

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

Построение линейной регрессии

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

Также можно провести оценку модели на кросс-валидации, тому же принципу, что был использован ранее. Для этого воспользуемся функцией (с небольшими модификациями), предложенной в посте Pythonic Cross Validation on Time Series

%%time
performTimeSeriesCV(X_train, y_train, 5, lr, mean_absolute_error)

На 5 фолдах получили среднюю абсолютную ошибку в 4.6 K пользователей, достаточно близко к оценке качества, полученной на тестовом датасете.

Код для построения прогноза с XGBoost

XGB_forecast(dataset, test_size=0.2, lag_start=5, lag_end=30)

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

Подготовка к анализу

Прежде чем мы начнем работу, нам нужно убедиться, что у нас есть все необходимые библиотеки для анализа временных рядов:

1. pandas: Pandas — это библиотека для работы с данными, которая предоставляет удобные структуры данных, такие как DataFrame, и множество функций для анализа и манипуляции данными. Она идеально подходит для работы с временными рядами, так как позволяет легко хранить и анализировать временные данные.

2. numpy: NumPy — это библиотека для работы с массивами и матрицами чисел. Она полезна при выполнении вычислений над данными временных рядов.

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

4. statsmodels: Statsmodels — это библиотека для статистического анализа данных. Она содержит множество методов для анализа временных рядов, включая модели ARIMA и SARIMA.

5. scikit-learn: Scikit-learn предоставляет множество инструментов для машинного обучения, включая модели регрессии и классификации, которые можно использовать для анализа временных рядов.

!pip install pandas numpy matplotlib seaborn statsmodels scikit-learn

Теперь, когда у нас есть необходимые библиотеки, давайте перейдем к работе с датами и временем.

Одной из ключевых характеристик временных рядов является наличие временных меток, которые позволяют нам связывать данные с конкретными моментами времени. При работе с временными рядами необходимо уметь обращаться с датами и временем.

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

Результат будет следующим:

Дата Продажи
0 2023-01-01 1000
1 2023-02-01 1200
2 2023-03-01 1300
3 2023-04-01 1100
4 2023-05-01 1400

Теперь у нас есть DataFrame с данными о продажах и столбцом ‘Дата’, который имеет тип datetime64. Это позволяет нам выполнять различные операции, такие как выбор данных по дате или вычисление временных интервалов.

Дата Продажи
2 2023-03-01 1300
3 2023-04-01 1100

Таким образом, работа с датами и временем в pandas позволяет нам легко фильтровать и анализировать данные временных рядов.

Обработка пропущенных значений

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

Давайте рассмотрим, как можно обработать пропущенные значения в DataFrame с помощью pandas. Для примера допустим, что у нас есть следующий dataset:

Дата Продажи
0 2023-01-01 1000.0
1 2023-02-01 NaN
2 2023-03-01 1300.0
3 2023-04-01 1100.0
4 2023-05-01 1400.0

Как видите, у нас есть пропущенное значение (NaN) в столбце ‘Продажи’. Существует несколько способов обработки таких значений:

1. Удаление строк с пропущенными значениями:

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

2. Замена пропущенных значений:

Здесь мы заменяем все пропущенные значения на 0. Вы можете выбрать другое значение для замены вместо 0 в зависимости от контекста.

Интерполяция позволяет заполнить пропущенные значения на основе соседних значений. Это может быть полезно, если данные имеют некоторую структуру.

Обработка пропущенных значений зависит от конкретной задачи и данных временных рядов. Важно выбирать подходящий метод в каждом случае.

Движемся, сглаживаем и оцениваем

Небольшое определение временного ряда:

Временной ряд – это последовательность значений, описывающих протекающий во времени процесс, измеренных в последовательные моменты времени, обычно через равные промежутки

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

Импортируем нужные библиотеки. В основном нам понадобится модуль statsmodels, в котором реализованы многочисленные методы статистического моделирования, в том числе для временных рядов. Для поклонников R, пересевших на питон, он может показаться очень родным, так как поддерживает написание формулировок моделей в стиле ‘Wage ~ Age + Education’.

import sys
import warnings
warnings.filterwarnings(‘ignore’)
from tqdm import tqdm

import pandas as pd
import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error

import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs
from scipy.optimize import minimize

import matplotlib.pyplot as plt

В качестве примера для работы возьмем реальные данные по часовому онлайну игроков в одной из мобильных игрушек:

Избавляемся от нестационарности и строим SARIMA

Попробуем теперь построить ARIMA модель для онлайна игроков, пройдя все круги ада стадии приведения ряда к стационарному виду. Про саму модель уже не раз писали на хабре — Построение модели SARIMA с помощью Python+R, Анализ временных рядов с помощью python, поэтому подробно останавливаться на ней не буду.

Код для отрисовки графиков

Out: Критерий Дики-Фуллера: p=0.190189

Как и следовало ожидать, исходный ряд стационарным не является, критерий Дики-Фуллера не отверг нулевую гипотезу о наличии единичного корня. Попробуем стабилизировать дисперсию преоразованием Бокса-Кокса.

Out: Критерий Дики-Фуллера: p=0.079760
Оптимальный параметр преобразования Бокса-Кокса: 0.587270

Уже лучше, однако критерий Дики-Фуллера по-прежнему не отвергает гипотезу о нестационарности ряда. А автокорреляционная функция явно намекает на сезонность в получившемся ряде. Возьмём сезонные разности:

Out: Критерий Дики-Фуллера: p=0.002571

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

Out: Критерий Дики-Фуллера: p=0.000000

Наконец, получили стационарный ряд, по автокорреляционной и частной автокорреляционной функции прикинем параметры для SARIMA модели, на забыв, что предварительно уже сделали первые и сезонные разности.

Начальные приближения Q = 1, P = 4, q = 3, p = 4

ps = range(0, 5)
d=1
qs = range(0, 4)
Ps = range(0, 5)
D=1
Qs = range(0, 1)

from itertools import product

parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)

Код для подбора параметров перебором

Лучшие параметры загоняем в модель:

Проверим остатки модели:

Что ж, остатки стационарны, явных автокорреляций нет, построим прогноз по получившейся модели

Код для построения прогноза и отрисовки графика

В финале получаем достаточно адекватный прогноз, в среднем модель ошибалась на 1.3 K пользователей, что очень и очень неплохо, однако суммарные затраты на подготовку данных, приведение к стационарности, определение и перебор параметров могут такой точности и не стоить.

Эконометрический подход

Снова небольшое лирическое отступление. Часто на работе приходится строить модели, руководствуясь одним основополагающим принципом – быстро, качественно, недорого. Поэтому часть моделей могут банально не подойти для «продакшн-решений», так как либо требуют слишком больших затрат по подготовке данных (например, SARIMA), либо сложно настраиваются (хороший пример – SARIMA), либо требуют частого переобучения на новых данных (опять SARIMA), поэтому зачастую гораздо проще бывает выделить несколько признаков из имеющегося временного ряда и построить по ним обычную линейную регрессию или навесить решаюший лес. Дешево и сердито.

Возможно, этот подход не является значительно подкрепленным теорией, нарушает различные предпосылки, например, условия Гаусса-Маркова, особенно пункт про некоррелированность ошибок, однако на практике нередко выручает и достаточно активно используется в соревнованиях по машинному обучению.

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

Добавлю только про еще один вариант кодирования категориальных признаков – кодирование средним. Если не хочется раздувать датасет множеством дамми-переменных, которые могут привести к потере информации о расстоянии, а в вещественном виде возникают противоречивые результаты а-ля «0 часов < 23 часа», то можно закодировать переменную чуть более интерпретируемыми значениями. Естественный вариант – закодировать средним значением целевой переменной. В нашем случае каждый день недели или час дня можно закодировать сооветствующим средним числом игроков, находившихся в этот день недели или час онлайн. При этом важно следить за тем, чтобы расчет среднего значения производился только в рамках тестового датасета (или в рамках текущего наблюдаемого фолда при кросс-валидации), иначе можно ненароком привнести в модель информацию о будущем.

Создадим новый датафрейм и добавим в него час, день недели и выходной в качестве категориальных переменных. Для этого переводим имеющийся в датафрейме индекс в формат datetime, и извлекаем из него hour и weekday.

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

code_mean(data, ‘weekday’, «y»)

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

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

Функция для создания переменных

Примеры на практике

В этом разделе мы рассмотрим несколько практических примеров использования анализа временных рядов для решения разнообразных задач.

Прогнозирование продаж на основе временных рядов продаж электроники

В этом примере мы будем использовать данные о ежемесячных продажах электроники в магазине МВидео для прогнозирования будущих продаж.

Cоздадим небольшой dataset с данными о продажах. В качестве исходных данных предположим следующие продажи за последние два года:

Продажи
Дата
2022-01-31 1355
2022-02-28 4665
2022-03-31 3154
2022-04-30 3490
2022-05-31 3569

Визуализация данных

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

Подготовка данных

ADF статистика: -3.336001912478917
p-значение: 0.013343910713214318
Критические значения:
1%: -3.859073285322359
5%: -3.0420456927297668
10%: -2.6609064197530863

Если p-значение ниже некоторого порогового значения (обычно 0.05), то мы можем считать ряд стационарным.

Выбор и обучение модели

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

Оценка качества прогноза

После обучения модели мы можем оценить ее качество на основе имеющихся данных.

MSE: 1178795.2830408707
MAE: 906.8571433244668

Прогноз на будущее

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

На графике показаны исходные данные о продажах и прогноз продаж на следующие 12 месяцев. Этот пример демонстрирует, как использовать анализ временных рядов для прогнозирования будущих продаж в магазине электроники.

Заключение

В данной статье мы рассмотрели основные концепции и методы работы с временными рядами в Python. Мы изучили, как импортировать данные временных рядов, визуализировать их, а также провели анализ стационарности и сезонности.

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