Прогноз сквозных временных рядов с использованием модели ARIMA

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

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

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

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

Я использую данные счетчика шагов, загруженные с моего iPhone Apple Health. Если вы также используете iPhone, вы можете загрузить данные о состоянии здоровья со страницы профиля вашего приложения Apple Health, как показано ниже.

Загруженные данные будут представлять собой файл XML, содержащий несколько секторов Apple Health, включая количество шагов, частоту сердечных сокращений, анализ сна и т. Д. Я использовал сценарий здесь из ссылки, чтобы преобразовать файл XML в формат CSV для анализа. Результатом этого файла будет несколько файлов CSV, и вы можете использовать Python для прямого подключения к файлу CSV. Это также забавный источник данных, где вы можете изучить другие способы анализа или визуализации данных с использованием ваших собственных данных о здоровье.

Исследование данных

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

Для анализа мне понадобятся только два столбца startDate и Value. Первый шаг для меня - преобразовать столбец даты в формат даты и времени, а затем объединить данные в недельную сумму. Это связано с тем, что для многих функций или пакетов Python для временных рядов потребуется один столбец в формате времени, и если я буду прогнозировать данные на ежедневной основе, волатильность будет слишком высокой.

#read data from extracted csv
steps=pd.read_csv('apple_health_export/StepCount.csv')
#convert start date into time format
steps['date']=pd.to_datetime(steps['startDate'].str[:19])
#Aggregate data into weekly sum
sample=steps[['date','value']]
weekly=sample.resample('W', on='date').sum()
#visualize weekly data
weekly.plot(figsize=(15, 6))
plt.show()

Из визуализации выше мы видим, что даже данные о количестве шагов за неделю имеют очень высокую волатильность. Чтобы устранить влияние экстремальных значений, таких как 300 000 шагов в неделю или менее 2000 шагов в неделю, я использую winsorization, чтобы сделать данные более нормальными.

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

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

Я выбираю верхний предел 2% и нижний 2% (это процентное значение может быть определено нами на основе распределения данных). Я использую пакет Python winsorize, в котором вы можете напрямую указать верхний и нижний процентный предел.

#remove the first part of data with no steps
sample=steps[steps['date']>'2015-09-01'][['date','value']]
#aggregate data on daily level
daily=sample.resample('D',on='date').sum()
#Winsorize daily data
daily['winsorized_value']=winsorize(daily['value'], limits=[0.02, 0.02])
#Aggregate daily data into weekly data again
weekly=daily.resample('W').sum()
#visualize new weekly data
weekly.plot(figsize=(15, 6))
plt.show()

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

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

extract=train.set_index('date')
from pylab import rcParams
rcParams['figure.figsize'] = 18, 8
decomposition = sm.tsa.seasonal_decompose(extract, model='additive')
fig = decomposition.plot()
plt.show()

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

Прогнозирование временных рядов

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

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

Подготовка данных. Этап немного отличается от обычного разделения теста на машинное обучение. Это связано с тем, что данные обучения и данные тестирования являются динамическими в зависимости от прогнозируемой недели (в качестве данных обучения мы используем все исторические данные до этой недели). Мы используем количество шагов в 2021 году в качестве ориентира для измерения точности нашей модели.

#determine training and testing group
data=weekly.reset_index()
test=data[data['date']>'2021-01-01'][['date','value']]
train=data[['date','value']]

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

Метод 1: скользящее среднее

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

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

#get rolling moving average for previous 10 weeks
data['SMA_10']=train['value'].rolling(window=10).mean().shift(1)
#Measure the MAE for this measure
test=data[data['date']>'2020-12-31']
from sklearn.metrics import mean_absolute_error
print(mean_absolute_error(test['value'], test['SMA_10']))
#10415.972

Метод 2: экспоненциально взвешенное скользящее среднее

Скользящее среднее по методу 1 - это простое среднее за последние 10 недель. Вы можете сказать, что неделя, ближайшая к дате предсказания, должна иметь более тяжелый вес, чем простое среднее значение. Именно тогда мы можем использовать экспоненциально взвешенное скользящее среднее.

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

#Calculate ewm average value
data['ewma_10']=train['value'].ewm(span=10).mean().shift(1)
test=data[data['date']>'2020-12-31']
from sklearn.metrics import mean_absolute_error
print(mean_absolute_error(test['value'], test['ewma_10']))
#9613.11

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

test_plot_data=test.set_index('date')[['SMA_10','ewma_10']]
train_plot=train[train['date']>'2020-01-01'].set_index('date')
plt.figure(figsize=(15,10))
plt.plot(train_plot,label='Actual Step')
plt.plot(test_plot_data['SMA_10'],label='sma_10')
plt.plot(test_plot_data['ewma_10'],label='ewma_10')
plt.legend(loc='Left corner')
plt.show()

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

Метод 3: Модель SARIMA - Ручная настройка

Модель SARIMA - это модель ARIMA с сезонным трендом. Полное название модели ARIMA - Авторегрессивная интегрированная скользящая средняя. Прежде чем я начну рассматривать метод, давайте кратко рассмотрим, что такое ARIMA и какие факторы мы учитываем при его моделировании.

Краткое резюме:

Авторегрессивный / AR: прогноз результатов зависит от предыдущих наблюдений / наблюдений с задержкой.

Интегрированный / I: данные не являются стационарными и требуют определенного порядка дифференцирования для достижения стационарности.

Скользящее среднее / скользящая средняя. Прогноз результатов зависит от предыдущих условий ошибок / запаздывающих ошибок.

Для модели SARIMA необходимо указать 7 различных параметров:

  • p: количество задержек для запроса дополненной реальности.
  • d: количество раз, необходимое для достижения стационарности.
  • q: количество запаздываний для срока скользящей средней.
  • s: сезонные периоды. Относится к количеству периодов времени, в течение которых тот же образец повторяется снова.
  • P, D, Q: то же, что и p, d, q, но это параметры для стороны сезонности.

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

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

# Define the p, d and q parameters for value between 0-2 and iterate for all the value in the range
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
# Generate all different combinations of seasonal p, q and q 
seasonal_pdq = [(x[0], x[1], x[2], 4) for x in list(itertools.product(p, d, q))]
params=[]
seasonal=[]
aic=[]
for param in pdq:
    for param_seasonal in seasonal_pdq:
        mod = sm.tsa.statespace.SARIMAX(train_set,order=param,                                 seasonal_order=param_seasonal,                                    
enforce_stationarity=False,                                          enforce_invertibility=False)
#append all the parameters and result AIC value
results = mod.fit()
        params.append(param)
        seasonal.append(param_seasonal)
        aic.append(results.aic)
parameter_options=pd.DataFrame({'params':params,'seasonal_params':seasonal,'AIC':aic})
#sort the AIC value to find the best fitted model
parameter_options.sort_values(by='AIC')

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

data_updated=data.set_index('date')
train=data[data['date']<'2021-01-01']
prediction=[]
my_order = (0, 1, 1)
my_seasonal_order = (1, 0, 1, 4)
initial=len(train)
initial_train=data_updated.iloc[:initial]
model = sm.tsa.statespace.SARIMAX(initial_train['value'], order=my_order, seasonal_order=my_seasonal_order)
results=model.fit()
#Iteratively update the training data and predict following week
for i in range(initial,len(data)):
    updated_data=data_updated.iloc[i:i+1]['value']
    results=results.append(updated_data,refit=False)
    prediction.append(results.forecast()[0])
from sklearn.metrics import mean_absolute_error
print(mean_absolute_error(test['value'], prediction))
#8469.6746

Модель работает намного лучше, чем результаты двух скользящих средних, с уменьшением значения средней ошибки на 15%.

Метод 4: Auto-ARIMA (автоматическая настройка параметров)

Метод 4 использует ту же модель ARIMA, но с другим пакетом, который обеспечивает функцию автонастройки. Я использую пакет pm.auto_arima, и вы можете узнать больше о функциях пакета по ссылке.

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

train_index=train.set_index('date')
#Spesify the value range for parameters
model = pm.auto_arima(train_index['value'], 
                          start_p=0, start_q=0,d=1, max_p=5, 
                          max_q=5, start_P=0, D=None, start_Q=0, max_P=5, 
                          max_D=5, max_Q=5,stepwise=True,seasonal=True)
prediction=[]
#Recurrently predict following week and add the data into training model once we have predicted that week
for i in test['value']:
    predict=model.predict(n_periods=1)[0]
    prediction.append(predict)
    model.update(i)
from sklearn.metrics import mean_absolute_error
print(mean_absolute_error(test['value'], prediction))
#8863.74

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

Когда я строю прогноз обоих методов SARIMA, я могу сказать, что эти два следуют более точно с реальным трендом подсчета шагов, чем простые или экспоненциально взвешенные методы скользящего среднего. Однако модели SARIMA также не могут очень точно уловить колебания. Частота ошибок составляет 20–30%.

Заключение или рекомендация

  1. В этом примере я периодически рассчитываю следующее недельное значение. Однако весьма вероятно, что вам потребуется спрогнозировать значение следующих N (N ›1) недель, используя исторические данные. На самом деле это намного проще, потому что вы можете указать n_periods на количество time_periods, которое вы хотите предсказать, как только вы получите модель (не нужно итеративно обновлять модель)
  2. В Интернете есть статьи о том, как определить стационарность данных или как определить ценность терминов AR и MA с помощью статистических методов, таких как «Расширенный тест Дики Фуллера». Это может лучше определить распределение данных и предотвратить переобучение.
  3. Подход SARIMA или ARIMA обычно работает лучше, чем подход простого скользящего среднего или экспоненциально взвешенного скользящего среднего, поскольку он учитывает член ошибки
  4. Одна очень важная вещь, на которую следует обратить внимание, заключается в том, что мы всегда должны помнить об использовании данных до даты, которую мы прогнозируем, в качестве обучающих данных, для нас будет сложно использовать перекрестную проверку.

И напоследок окончательный вывод: очень сложно предсказать количество шагов с высокой точностью!

Ссылка или рекомендуемая литература

  1. Https://www.digitalocean.com/community/tutorials/a-guide-to-time-series-forecasting-with-arima-in-python-3 (как выбрать параметры ARIMA с использованием самого низкого AIC)
  2. Https://www.machinelearningplus.com/time-series/arima-model-time-series-forecasting-python/ (как определить значение параметра статистически)
  3. Https://github.com/markwk/qs_ledger/blob/master/apple_health/apple_health_extractor.ipynb. (код для извлечения данных из Apple Health)
  4. Https://towardsdatascience.com/what-are-the-best-metrics-to-evaluate-your-regression-model-418ca481755b (как оценить вашу регрессионную модель)
  5. Https://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html (пакет модели SARIMA, ручная настройка)
  6. Https://alkaline-ml.com/pmdarima/modules/generated/pmdarima.arima.auto_arima.html (пакет модели SARIMA, автонастройка)

Спасибо за прочтение! Если у вас есть какие-либо вопросы или темы, которые вы хотите понять больше, укажите их в поле для комментариев!