В прошлом посте я рассказывала про свои мучения с моделькой ARIMA. Здесь же я расскажу о следующей серии издевательств над временными рядами, SARIMAX и экспоненциальным сглаживанием.
Для начала хочу исправить косяки прошлой статьи и проговорить базовые вещи. Можно смело пропускать этот абзац, если Вы и так шарите в основах временных рядов.
Занудная математика
Гауссов(ский) процесс: если конечные распределения процесса ξt для любых t1,. . ., tn являются нормально распределенными (гауссовскими), то процесс ξt называется гауссовским.
Как справедливо заметили в комментариях к предыдущему посту, очень часто гауссовость случайной части зашита во многие модели, которые к процессу применяются. Не рискну лезть вглубь математики в этой части, но замечание крайне ценное. После этого я стала смотреть не только на поведение модели на тестовой выборке, но и на распределение остатков (residuals). Подробнее можно почитать здесь. За наводку спасибо @adeshere.
Временной ряд: значения меняющихся во времени признаков, полученные в некоторые моменты времени.Временные ряды делятся на стационарные и нестационарные, и работа с ними принципиально отличается. Подробнее про математическую базу можно почитать здесь.
Единственный недостаток Яндексовского учебника по ML — очень мало примеров кода. Мне, как новичку, было крайне тяжело связать математику и код с первого раза.
Прочие косяки прошлой статьи
Также в прошлой статье в комментах заметили, что ACF и PACF не гарант хорошего подбора гиперпараметров. Я в этом убедилась на собственном печальном опыте, но ещё раз это подчеркну.
К сожалению, пока не было времени попробовать инструменты из R, модель ARCH, Pycaret и fbprophet.
Данные
У меня была статистика компании за 2.5 года по дням, для пары экспериментов я пробовала склеить её в недельные данные. К сожалению или к счастью никакой физики за процессом тут не лежит. Требовалась просто сносная модель, которая бы предсказывала данные на следующий месяц и при этом давала результат лучше, чем предсказания аналитиков через конверсии. И была возможность вытащить из данных ещё дополнительную колонку внешних параметров.
SARIMA
На модель SARIMA у меня были большие надежды, всё-таки в отличие от ARIMA она могла оперировать внешними параметрами и с сезонностью у неё дела обстояли сильно лучше. В итоге я остановилась на статистике по дням.
Чтобы работать со стационарным рядом, я дифференцировала изначальный:
diff1 = df.copy(deep=True)
diff1['1'] = df['1'].diff()
diff1['2'] = df['2'].diff()
diff1.dropna(inplace=True)Затем сделала пробный тычок пальцем в небо, чтобы посмотреть на работу модели:
from statsmodels.tsa.statespace.sarimax import SARIMAX
# Ensure the SARIMAX model is initialized properly with your training data and exogenous variables
model = SARIMAX(diff1['1'][:962], exog=diff1['2'][:962], order=(1, 1, 1), seasonal_order=(1, 0, 1, 12))
# Fit the model
result = model.fit(disp=False)
# Forecast using the correct steps and aligned exogenous data
forecast = result.get_forecast(steps=len(df_test), exog=forecast_exog)
predictions = result.forecast(len(df_test), exog=forecast_exog)
# predictions
predictions = pd.Series(predictions.values, index=df_test.index)
# If the forecast returns NaNs, check alignment and data quality
forecast_df = pd.DataFrame({
'predicted_mean': forecast.predicted_mean,
'lower_ci': forecast.conf_int().iloc[:, 0],
'upper_ci': forecast.conf_int().iloc[:, 1]
}, index=df_test.index)Оптимизация с помощью Optuna:
#Generate all different combinations of p, d and q triplets
#Generate all different combinations of p, d, q and s triplets
import itertools
import optuna
p = range(1, 3) # Smaller range for non-seasonal AR component
d = range(0, 2)
q = range(0, 3)
P = range(0, 2) # Smaller range for seasonal AR component
D = range(0, 2)
Q = range(0, 2)
s = 12
pdq = list(itertools.product(p, d, q))
pdqs = [(x[0], x[1], x[2], s) for x in list(itertools.product(P, D, Q))]
#%%
def objective_sarima(trial):
non_seasonal_order = trial.suggest_categorical('non_seasonal_order', pdq)
seasonal_order = trial.suggest_categorical('seasonal_order', pdqs)
trend = trial.suggest_categorical('trend', ['n', 'c', 't', 'ct', None])
# # Generate predictions
# predictions = mdl.forecast(len(df_test))
model = SARIMAX(diff1['1'][:962], exog=diff1['2'][:962], order=non_seasonal_order, seasonal_order=seasonal_order)
# Fit the model
result = model.fit(disp=True)
# Forecast using the correct steps and aligned exogenous data
forecast = result.get_forecast(steps=len(df_test), exog=forecast_exog)
predictions = result.forecast(len(df_test), exog=forecast_exog)
# predictions
predictions = pd.Series(predictions.values, index=df_test.index)
# If the forecast returns NaNs, check alignment and data quality
forecast_df = pd.DataFrame({
'predicted_mean': forecast.predicted_mean,
'lower_ci': forecast.conf_int().iloc[:, 0],
'upper_ci': forecast.conf_int().iloc[:, 1]
}, index=df_test.index)
# Calculate residuals and error metric
residuals = diff1['1'] - predictions
mse = np.sqrt(np.mean(residuals**2))
return mse
study=optuna.create_study(direction="minimize")
study.optimize(objective_sarima,n_trials=10)Затем я на всё это счастье визуализировала и смотрела. Получилось устрашающе:
plt.figure(figsize=(10, 6))
plt.plot(df_test['1'], label='Actual Future', marker='o', color='green')
plt.plot(forecast_df['predicted_mean'], label='Forecasted', marker='o', color='red')
plt.fill_between(forecast_df.index, forecast_df['lower_ci'], forecast_df['upper_ci'], color='red', alpha=0.2)
plt.title('Forecasted vs Actual Counts')
plt.xlabel('Date')
plt.ylabel('Count')
plt.legend()
plt.grid(True)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()
После я сравнила ещё и обратно проинтегрированный ряд с исходными данными, которые мне предоставили аналитики. Меня просили оформить всё с табличкой-сравнением, так что код получился монструозный. Тут важно опять-таки не забывать про константу интегрирования, которая в коде носит имя offset. Иначе это чеховское ружьё выстрелит в ногу.
#integral constant!
offset = int(list(dataset_sends['count'])[-len(df_test)])
print(offset)
inverted_forecast = pd.Series(forecast_df['predicted_mean']).cumsum()
inverted_forecast = pd.DataFrame(forecast_df['predicted_mean'])
inverted_true = pd.DataFrame(list(dataset_sends['count'][-len(df_test):]), index=range(0, len(df_test)))
inverted_forecast = inverted_forecast + offset
inverted_forecast = pd.DataFrame(list(inverted_forecast['predicted_mean'][-len(df_test):]), index=range(0, len(df_test)))
dates = pd.DataFrame(list(dataset_sends['date'][-len(df_test):]), index=range(0, len(df_test)))
df_merged = pd.concat([inverted_forecast, inverted_true, dates], axis=1)
df_merged.columns = ['predicted_mean', 'count', 'date']
df_merged.dropna(inplace=True)
df_merged['1'] = df_merged['1'].apply(int)
df_merged['predicted_mean'] = df_merged['predicted_mean'].apply(int)
# Absolute Error
df_merged['error'] = abs(df_merged['count'] - df_merged['predicted_mean'])
df_merged['date'] = pd.to_datetime(df_merged['date'])
# Set 'date' column as the index
df_merged.set_index('date', inplace=True)
# Group by week and sum the other columns
df_merged = df_merged.resample('W').sum()
# Reset index to make 'date' a column again
df_merged.reset_index(inplace=True)
# Relative Error
df_merged['relative_error, %'] = abs(df_merged['count'] - df_merged['predicted_mean']) / df_merged['count'] * 100
df_merged.dropna(inplace=True)
display(df_merged[-17:-1])
plt.plot(df_merged.index, df_merged['count'], label='count')
# Plot 'predicted_mean'
plt.plot(df_merged.index, df_merged['predicted_mean'], label='predicted_mean')
# Set labels and title
plt.xlabel('Time')
plt.ylabel('Count')
plt.title('Comparison of count and predicted_mean')
plt.legend()
# Show plot
plt.show()Увы, показать табличку я не могу из-за NDA, и так хожу по офигенно тонкому льду:)
Вывод: хотя я и питала большие надежды, результат SARIMA меня не удовлетворил. И я пошла гуглить дальше, что ещё умного делают люди для предсказания временных рядов.
Экспоненциальное сглаживание
Этот подход упоминают в Яндекс-учебнике по ML, но уделяют не так много внимания, как моделям ARIMA/SARIMA. А вот мне он дал неожиданно классные результаты.
Здесь я тоже юзала пакетную историю и работала с данным, лежащими в переменной aust. Также руками поставила период в 40 дней, который подбирала ручками, без всяких оптимизаторов. Не делайте так:)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
import statsmodels.api as sm
from pandas.plotting import register_matplotlib_converters
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error,mean_absolute_error
import warnings
from statsmodels.tsa.holtwinters import SimpleExpSmoothing, Holt
warnings.filterwarnings('ignore')
period = 40
fit1 = ExponentialSmoothing(
aust,
seasonal_periods=period,
trend="add",
seasonal="add",
use_boxcox=True,
initialization_method="estimated",
).fit()
fit2 = ExponentialSmoothing(
aust,
seasonal_periods=period,
trend="add",
seasonal="mul",
use_boxcox=True,
initialization_method="estimated",
).fit()
fit3 = ExponentialSmoothing(
aust,
seasonal_periods=period,
trend="add",
seasonal="add",
damped_trend=True,
use_boxcox=True,
initialization_method="estimated",
).fit()
fit4 = ExponentialSmoothing(
aust,
seasonal_periods=period,
trend="add",
seasonal="mul",
damped_trend=True,
use_boxcox=True,
initialization_method="estimated",
).fit()
results = pd.DataFrame(
index=[r"$\alpha$", r"$\beta$", r"$\phi$", r"$\gamma$", r"$l_0$", "$b_0$", "SSE"]
)
params = [
"smoothing_level",
"smoothing_trend",
"damping_trend",
"smoothing_seasonal",
"initial_level",
"initial_trend",
]
results["Additive"] = [fit1.params[p] for p in params] + [fit1.sse]
results["Multiplicative"] = [fit2.params[p] for p in params] + [fit2.sse]
results["Additive Dam"] = [fit3.params[p] for p in params] + [fit3.sse]
results["Multiplica Dam"] = [fit4.params[p] for p in params] + [fit4.sse]
ax = train_data.plot(marker="o",
color="black",
title="Forecasts from Holt-Winters' multiplicative method",
)
ax.set_xlim(1697587200000000000, 1713139200000000000)
ax.set_ylim(0, 20)
ax.set_ylabel("counts")
ax.set_xlabel("Year")
fit1.fittedvalues.plot(style="--",marker="o", color="red", ax=ax)
fit2.fittedvalues.plot(style="--", marker="o",color="green", ax=ax)
fit3.fittedvalues.plot(style="--", marker="o",color="blue", ax=ax)
fit4.fittedvalues.plot(style="--", marker="o",color="green", ax=ax)
forecast1 = fit1.forecast(121).rename("Holt-Winters (add-add-seasonal)")
forecast2 = fit2.forecast(121).rename("Holt-Winters (add-mul-seasonal)")
forecast3 = fit3.forecast(121).rename("Holt-Winters (add-add-seasonal heuristic)")
forecast4 = fit4.forecast(121).rename("Holt-Winters (add-mul-seasonal heuristic)")
predictions1 = pd.Series(forecast1.values, index=test_data.index)
predictions2 = pd.Series(forecast2.values, index=test_data.index)
predictions3 = pd.Series(forecast3.values, index=test_data.index)
predictions4 = pd.Series(forecast4.values, index=test_data.index)
predictions1.dropna(inplace=True)
predictions2.dropna(inplace=True)
predictions3.dropna(inplace=True)
predictions4.dropna(inplace=True)
print("1mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions1)], predictions1), 5))
print("1mean squared error : ",round(mean_squared_error(test_data[:len(predictions1)], predictions1), 5))
print("1Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions1)], predictions1,squared=False), 5))
print("2mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions2)], predictions2),5))
print("2mean squared error : ",round(mean_squared_error(test_data[:len(predictions2)], predictions2),5))
print("2Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions2)], predictions2,squared=False),5))
print("3mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions3)], predictions3),5))
print("3mean squared error : ",round(mean_squared_error(test_data[:len(predictions3)], predictions3),5))
print("3Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions3)], predictions3,squared=False),5))
# print("4mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions4)], predictions4),5))
# print("4mean squared error : ",round(mean_squared_error(test_data[:len(predictions4)], predictions4),5))
# print("4Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions4)], predictions4,squared=False),5))
forecast1.plot(ax=ax, style="--", marker="o", color="green", legend=True, figsize=(20, 10))
forecast2.plot(ax=ax, style="--", marker="o", color="red", legend=True)
forecast3.plot(ax=ax, style="--", marker="o", color="blue", legend=True)
forecast4.plot(ax=ax, style="--", marker="o", color="yellow", legend=True)
plt.plot(predictions1, marker="o", color = "green")
plt.plot(predictions2, marker="o", color = "red")
plt.plot(predictions3, marker="o", color = "blue")
plt.plot(predictions4, marker="o", color = "yellow")
plt.plot(test_data, marker="o", color = 'black')
plt.plot(train_data, marker="o", color = 'black')
# Show the plot
plt.show()
print("Figure 7.6: Forecasting.")
# resultsКод тоже выглядит монструозно, я надеялась его переписать, но посыпались новые задачи и я забила.
К тому же я не смогла по-человечески вписать сюда историю с адекватным временем. Мне пришлось в самом начале перевести всё в цифирьки и довольствоваться малым:
aust['date'] = pd.to_numeric(pd.to_datetime(aust['date']))Исходная картинка у меня пропала, но картинка была какая-то такая:
Затем я взяла среднее по всем эти предсказаниям и получила вполне сносную картинку
combined_predictions = pd.DataFrame({
'predictions1': predictions1,
'predictions2': predictions2,
'predictions3': predictions3,
'predictions4': predictions4
})
# Calculate the average prediction across all models
average_prediction = combined_predictions.mean(axis=1)
from sklearn.metrics import mean_absolute_percentage_error
plt.figure(figsize=(15, 5))
plt.plot(average_prediction, marker="o", color = "red")
plt.plot(test_data, marker="o", color = 'black')
plt.ylim(0, 20)
plt.show()
print("mean absolute error : ",round(mean_absolute_error(test_data, average_prediction),5))
print("mean squared error : ",round(mean_squared_error(test_data, average_prediction),5))
print("Root mean squared error : ",round(mean_squared_error(test_data, average_prediction,squared=False),5))
print("mean relative error", round(mean_absolute_percentage_error(test_data, average_prediction)), "%")
красные точки — предсказания, чёрные — реальные данные
Вывод:
1. Вероятно, такую хорошую картинку я получила именно из-за сглаживания. Про это писали среди комментов к прошлому посту.
2. Временную шкалу мне пришлось дополнительно обрабатывать, чтобы всё заработало. Можно ли было обойтись без этого — не знаю.
В заключительной части этой "временной" серии расскажу про мой опыт использования пакета FEDOT. И если успею потестить что-то из предложенного ранее в комментах — тоже напишу.
Ермак Марина
Аналитик, SENSE IT
