В прошлом посте я рассказывала про свои мучения с моделькой 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
