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

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

Есть много способов решить эту проблему, в этом посте я сосредоточусь на разговоре о некоторых методологиях динамических систем, которые мы можем использовать для применения стандартных алгоритмов обучения, таких как SVM или Gradient Boosting, к данным временных рядов. Другими возможностями такого прогнозирования являются:

  • Статистические методы прогнозирования, такие как ARIMA
  • Методологии глубокого обучения, такие как LSTM и RNN

Во-первых, давайте поймем, почему не стоит брать наши данные временных рядов и напрямую применять их к регрессору.

Коды, используемые для генерации результатов для этого поста (и других), доступны на моем GitHub и в Kaggle.

Почему алгоритмы по умолчанию обычно не работают

Статистическая теория обучения является основой для математического доказательства того, что наши модели машинного обучения действительно способны к обучению (также известному как обобщение). Чтобы эти гарантии сохранялись, должны быть установлены некоторые предположения о данных, на которых мы пытаемся моделировать [1], одно из них:

  • Примеры должны быть выбраны i.i.d.

Когда мы имеем дело с временным рядом, примеры нашего набора данных почти не зависят друг от друга (если они независимы, то ваш ряд, вероятно, является случайным рядом).

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

Теорема вложения Такена

Эта теорема была предложена Такенсом [2] в 80-х годах и говорит о том, как встраивание (переход из одного пространства в другое) нашего временного ряда в пространство его временных задержек, которое мы называем фазовым пространством, устраняет временную зависимость между экземпляры серии.

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

Конечно, не все наши временные ряды будут детерминированными, и поэтому аттракторы в фазовом пространстве не будут так четко определены.

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

  • m -> Называемый размерностью встраивания, этот параметр сообщит нам размерность фазового пространства
  • d -> Называется временной задержкой, это говорит нам, сколько временных задержек будет представлять каждая ось фазового пространства.

Звучит немного абстрактно, правда? Давайте посмотрим на это в виде уравнения, чтобы все стало ясно для понимания.

Пусть f(t) будет нашим временным рядом. Учитывая m и d, наше вложение фазового пространства будет следующей матрицей:

Это означает, что каждая из осей m нашего пространства будет представлена ​​как кратное заданному нами временному отставанию.

Реализация встраивания взятого

Давайте реализуем некоторый код для преобразования временного ряда в наше фазовое пространство. Во-первых, давайте импортируем необходимые библиотеки:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from nolitsa import dimension, delay
import plotly.express as px
import plotly.graph_objects as go
from sktime.datasets import load_airline, load_shampoo_sales, load_lynx
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.model_selection import train_test_split
from gtda.time_series import SingleTakensEmbedding, takens_embedding_optimal_parameters
from openml.datasets.functions import get_dataset

Для реализации Taken мы можем использовать только numpy:

def takens(data, m=2, d=1):
    emb = np.array([data[0:len(data) - d*m]])
    for i in range(1, m):
        emb = np.append(emb, [data[i*d:len(data) - d*(m - i)]], axis=0)
        
    return emb.T

Получим функцию Лоренца из пакета Giotto-TDA:

lorenz = get_dataset(42182).get_data(dataset_format='array')[0][0:, 0]
t = [i for i in range(len(lorenz))]
emb = takens(lorenz, m=3, d=5)
fig = go.Figure()
fig.add_trace(go.Scatter3d(
    x=emb[:, 0], y=emb[:,1], z=emb[:, 2], mode='lines'
))
fig.show()

Код сгенерировал следующий график:

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

Конечно, вам может быть интересно, как я узнал, какие параметры следует использовать для m и d? Что ж, функция Лоренца действительно хорошо изучена, и эти значения были найдены много лет назад.

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

Нахождение m и d

Чтобы найти эти параметры, мы будем использовать некоторые реализации из пакета NoLiTSA для Python. Этот пакет реализует несколько функций и алгоритмов для нелинейного анализа временных рядов. Если у вас еще не установлен этот пакет, вы можете установить его с помощью:

pip install git+https://github.com/manu-mannattil/nolitsa.git

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

Следующий фрагмент кода сгенерирует взаимную информацию для синуса и построит его:

t = [i for i in np.arange(0, 20, 0.1)]
y = [np.sin(i) for i in t]
plt.figure(figsize=(9,4))
plt.xlabel('Time delay (d)')
plt.ylabel('Mutual Information')
plt.plot(delay.dmi(y, maxtau=20))

Дайте нам сюжет:

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

def find_optimal_delay(x, maxtau=50):
    mi = delay.dmi(x, maxtau=maxtau)
    diffmi = np.diff(mi)
return np.where(diffmi > 0)[0][0]

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

Эта метрика сообщает нам, учитывая размерность, сколько соседей наших точек являются «ложными соседями». Ложный сосед — это точка, которая находится по соседству с данной точкой в ​​данном измерении, но ее больше нет после того, как мы развернём пространство в другое измерение.

Мы можем построить пропорцию ложных ближайших соседей, используя следующее:

dim = np.arange(1, 10)
f1, f2, f3 = dimension.fnn(y, tau=4, dim=dim)
plt.figure(figsize=(9,4))
plt.xlabel('Dimension (m)')
plt.ylabel('False Nearest Neighbors')
plt.plot(f1)

А сюжет таков:

Эвристика в этом случае заключается в том, что любое измерение, в котором результат FNN меньше 20%, является хорошим кандидатом. Итак, мы можем создать нашу функцию, чтобы найти его:

def find_optional_dimension(x, tau, max_dim=10):
    dim = np.arange(1, max_dim)
    f1, f2, f3 = dimension.fnn(x, tau=tau, dim=dim)
return np.where(f1 < 0.2)[0][0] + 1

Обратите внимание, что это эмпирические эвристики, поэтому они не гарантируют нахождения наилучших значений. Кроме того, если использовать автоматический поиск из пакета Giotto-TDA (в котором уже реализованы Taken), нет гарантии, что значения совпадут.

Давайте теперь посмотрим, как использовать все это для прогнозирования нашего временного ряда.

Прогнозирование временных рядов с помощью встраивания

Хорошо, идея здесь в том, что после встраивания у нас есть X-матрица с задержками нашего временного ряда в новое пространство. Мы можем использовать первые m-1 записи в этой матрице в качестве признаков для нашего регрессора, чтобы предсказать последний столбец.

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

def forecast_on_phase_space(y, m, d, forecast_horizon):
    emb = takens(y, m, d)
    
    # Divide into train and test
    X = emb[:, :m-1]
    y = emb[:, m-1]
    
    X_train = X[:len(X)-forecast_horizon, :]
    y_train = y[:len(y)-forecast_horizon]
    X_test = X[len(X)-forecast_horizon:, :]
    y_test = y[len(y)-forecast_horizon:]
    
    # Fit the regressor on the training data
    rf = RandomForestRegressor()
    rf.fit(X_train, y_train)
    
    # Predict the test data
    preds = rf.predict(X_test)
    
    print(f'R²: {r2_score(y_test, preds)}')
    print(f'RMSE: {mean_squared_error(y_test, preds, squared=False)}')
    
    # Plot the result
    preds_ = [np.nan for i in range(len(y)-forecast_horizon)] + list(preds)
    t = [i for i in range(len(y))]
    fig = go.Figure()
    fig.add_trace(go.Scatter(
        x=t, y=y, mode='lines'
    ))
    fig.add_trace(go.Scatter(
        x=t, y=preds_, mode='lines', line=dict(color='red')))
    fig.show()

Следующий код сгенерирует следующий прогноз:

Кажется довольно хорошим, не так ли? Что ж, давайте теперь посмотрим, что произойдет, если мы используем реальный набор данных, который не является таким уж детерминированным. Для этого установим библиотеку sktime и воспользуемся данными Airline [3], свободно доступными в библиотеке:

!pip install sktime
airline = load_airline()
forecast_on_phase_space(airline, 3, 1, 50)

Что генерирует:

Не так хорошо, верно? Что ж, мы не проводили никакой предварительной обработки в сериале, мы могли бы многое сделать, чтобы улучшить его, однако это будет рассмотрено в следующем сюжете.

Теперь можно спросить: хорошо, а что, если я хочу предсказать несколько новых значений из неизвестного ряда, не имея для этого задержек? В вашем примере вы использовали задержку, но в реальном сценарии теперь у меня будет X для подачи регрессора на .predict().

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

Заключение

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

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

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

Благодарности

Здесь я должен поблагодарить Родриго Мелло, бывшего профессора Университета Сан-Паулу, который впервые познакомил меня с этими концепциями еще в 2019 году во время занятий в университете.

[1] Палиоза, Лукас и Мелло, Родриго. (2016). Применение функции ядра к данным, зависящим от времени, для обеспечения гарантий обучения с учителем. Экспертные системы с приложениями. 71. 10.1016/j.eswa.2016.11.028.

[2] Takens F. (1981) Обнаружение странных аттракторов в турбулентности. В: Рэнд Д., Янг Л.С. (редакторы) Динамические системы и турбулентность, Уорик, 1980. Конспект лекций по математике, том 898. Springer, Берлин, Гейдельберг. https://doi.org/10.1007/BFb0091924

[3] Box, G.E.P., Jenkins, G.M. и Reinsel, G.C. (1976) Анализ временных рядов, прогнозирование и контроль. Третье издание. Холден-Дэй. Серия Г.