Начинаю серию статей о сравнительно молодой библиотеки ETNA для прогнозирования временных рядов от команды Тинькофф Банк. В работе над статьями использую очень полезную книгу Груздев А. В. "Прогнозирование временных рядов с помощью Facebook Prophet, ETNA,sktime и LinkedIn Greykite".Начинаю с самой простой модели : модель наивного прогноза
import pandas as pd
import numpy as np
import re
import datetime
Модуль re (регулярные выражения) в Python предоставляет мощные инструменты для работы с текстом. Регулярные выражения позволяют искать, заменять и делить текст на основе паттернов или шаблонов. Модуль datetime - работа с временными метками.
Далее по порядку импортируем функции и классы библиотеки ETNA.
# импортируем модели наивного прогноза
from etna.models import NaiveModel
# импортируем классы для вычисления метрик
from etna.metrics import MAE, MSE, SMAPE, MAPE
# импортируем функцию plot_forecast() для визуализации прогнозов
from etna.analysis import plot_forecast
Удобство библиотеки ETNA заключается в том, что мы можем применить выбранную модель прогноза сразу к нескольким временным рядам и получить прогнозы для каждого из них. Будем работать с набором данных, представляющим месячные выручки десяти розничных магазинов одежды в период январь 2015 декабрь 2019.
Загружаем данные из excel файла :
df = pd.read_excel('data/shops_10_rev_month.xlsx')
df.head()
month | shop | rev | |
---|---|---|---|
0 | 2015-01-01 | Sh0 | 3040820 |
1 | 2015-01-01 | Sh1 | 3990770 |
2 | 2015-01-01 | Sh2 | 4356330 |
3 | 2015-01-01 | Sh3 | 4630860 |
4 | 2015-01-01 | Sh4 | 2450160 |
ETNA требует определенного формата данных: столбец, который мы хотим
спрогнозировать, должен называться target, столбец с временными метками
должен называться timestamp. Поскольку библиотека может работать с несколькими
временными рядами, столбец segment является меткой ряда.
# задаем обязательные столбцы target, timestamp и segment
df=df.rename(columns={"month":"timestamp", "shop":"segment", "rev":"target"})
df['timestamp'] = pd.to_datetime(df['timestamp'])
df.head()
timestamp | segment | target | |
---|---|---|---|
0 | 2015-01-01 | Sh0 | 3040820 |
1 | 2015-01-01 | Sh1 | 3990770 |
2 | 2015-01-01 | Sh2 | 4356330 |
3 | 2015-01-01 | Sh3 | 4630860 |
4 | 2015-01-01 | Sh4 | 2450160 |
Библиотека работает со специальной структурой данных TSDataset. Класс
TSDataset – основной класс библиотеки ETNA для работы с временными рядами.
Его главные методы:
.to_dataset() – превращает датафрейм c плоским индексом в датафрейм
c мультииндексом;
.to_pandas() – превращает объект TSDataset в датафрейм c плоским индексом
(flatten=True) или датафрейм c мультииндексом (flatten=False);
.fit_transform() – вычисляет и применяет преобразования и конструирование
признаков (например, логарифмирует зависимую переменную,
прогнозирует тренд и удаляет тренд из зависимой переменной, вычисляет
скользящее среднее с шириной окна 20 дней и создает столбец с
ним), требует списка экземпляров классов-трансформеров;
.inverse_transform() – применяет обратные преобразования (например,
выполняет экспоненцирование прологарифмированной зависимой переменной,
добавляет тренд к зависимой переменной c удаленным трендом),
требует списка экземпляров классов-трансформеров;
.make_future() – создает тестовый набор / набор новых данных, увеличив
обучающий/исторический набор данных на длину горизонта прогнозирования
и применив преобразования и конструирование признаков,
требует списка экземпляров классов-трансформеров;
.train_test_split() – разбивает на обучающую и тестовую выборку c
учетом временной структуры;
.describe() – вывод описательных статистик;
.plot() – визуализирует временной ряд в переменной target (в зависимости
от количества сегментов может быть визуализирован один
или несколько временных рядов).
А сейчас с помощью класса TSDataset превращаем наш датафрейм в объект
TSDataset, указав частоту временного ряда.
ts = TSDataset(df, freq='MS')
ts.head()
segment | Sh0 | Sh1 | Sh2 | Sh3 | Sh4 | Sh5 | Sh6 | Sh7 | Sh8 | Sh9 |
---|---|---|---|---|---|---|---|---|---|---|
feature | target | target | target | target | target | target | target | target | target | target |
timestamp | ||||||||||
2015-01-01 | 3040820 | 3990770 | 4356330 | 4630860 | 2450160 | 3858160 | 4566790 | 4626270 | 6468900 | 4189940 |
2015-02-01 | 2409980 | 3384060 | 3683360 | 3310420 | 1931430 | 2959990 | 3766150 | 3814750 | 4761640 | 3838840 |
2015-03-01 | 2604970 | 3856440 | 4419200 | 4434280 | 2325200 | 3457290 | 4352330 | 4308780 | 5433960 | 4684460 |
2015-04-01 | 3001640 | 4095820 | 5134330 | 5254790 | 2574390 | 3861470 | 4510820 | 4924240 | 6658410 | 5339020 |
2015-05-01 | 3065120 | 3954920 | 5248450 | 5156630 | 2380390 | 3848130 | 4216470 | 4949780 | 6724790 | 4774140 |
С помощью метода .describe() можем посмотреть базовые характеристики
объекта TSDataset:
start_timestamp – стартовая дата временного ряда;
end_timestamp – конечная дата временного ряда;
length – длина временного ряда;
num_missing – количество пропусков;
num_segments – количество сегментов (рядов);
num_exogs – количество экзогенных переменных;
num_regressors – количество регрессоров;
num_known_future – количество регрессоров, добавленных из датафрейма
df_exog;
freq – частота временного ряда.
ts.describe().head().T
segments | Sh0 | Sh1 | Sh2 | Sh3 | Sh4 |
---|---|---|---|---|---|
start_timestamp | 2015-01-01 00:00:00 | 2015-01-01 00:00:00 | 2015-01-01 00:00:00 | 2015-01-01 00:00:00 | 2015-01-01 00:00:00 |
end_timestamp | 2019-12-01 00:00:00 | 2019-12-01 00:00:00 | 2019-12-01 00:00:00 | 2019-12-01 00:00:00 | 2019-12-01 00:00:00 |
length | 60 | 60 | 60 | 60 | 60 |
num_missing | 0 | 0 | 0 | 0 | 0 |
num_segments | 10 | 10 | 10 | 10 | 10 |
num_exogs | 0 | 0 | 0 | 0 | 0 |
num_regressors | 0 | 0 | 0 | 0 | 0 |
num_known_future | 0 | 0 | 0 | 0 | 0 |
freq | MS | MS | MS | MS | MS |
#Давайте визуализируем временной ряд.
# визуализируем временной ряд
ts.plot()
Прогнозирование временных рядов нужно начинать с самых простых моделей.
Прогнозирование начнем с самой простой модели - модели наивного прогноза, т. е. в качестве прогноза берем последнее известное значение. Для этого воспользуемся классом NaiveModel –моделью наивного прогноза и обучим ее.
Горизонт прогнозирования зададим равным 12 месяцам.
# задаем горизонт прогнозирования
HORIZON = 12
Разбиваем набор
(наш объект TSDataset) на обучающую и тестовую выборки с учетом временной
структуры, так чтобы последние 12 месяцев ушли в тестовую выборку.
train_ts, test_ts = ts.train_test_split(
train_start='2015-01-01',
train_end='2018-12-01',
test_start='2019-01-01',
test_end='2019-12-01')
Количество лагов в модели надо задавать равным горизонту прогнозирования и выше его, чтобы избежать повторения одних и тех же прогнозов.
Создаем модель наивного прогноза, задав 12 лагов, т. е. по длине горизонта.
Теперь прогнозами будут последние 12 известных значений.
# создаем экземпляр класса NaiveModel, в котором
# реализована модель наивного прогноза,
# задав 12 лагов
model = NaiveModel(lag=12)
# обучаем модель наивного прогноза
model.fit(train_ts)
С помощью метода .make_future() формируем набор, для которого нужно
получить прогнозы (по сути, тестовую выборку). Параметр future_steps задает
количество моментов времени, которые мы хотим предсказать в будущем
(горизонт прогнозирования). Параметр tail_steps задает количество моментов
времени, предшествующих горизонту прогнозирования, которые мы хотим
использовать для прогнозов (контекст модели). Этот параметр требуется
для простых моделей типа наивного прогноза, сезонного наивного прогноза,
скользящего среднего, сезонного скользящего среднего.
future_ts = train_ts.make_future(future_steps=HORIZON,
tail_steps=model.context_size)
С помощью метода .forecast() получаем наши прогнозы. В метод .forecast()
мы, помимо датафрейма, передаем параметр prediction_size, потому что
прогнозы определяются контекстом модели (количеством моментов времени,
предшествующих горизонту прогнозирования, которые мы хотим использовать
для прогнозов).
# получаем прогнозы
forecast_ts = model.forecast(future_ts,
prediction_size=HORIZON)
forecast_ts
segment | Sh0 | Sh1 | Sh2 | Sh3 | Sh4 | Sh5 | Sh6 | Sh7 | Sh8 | Sh9 |
---|---|---|---|---|---|---|---|---|---|---|
feature | target | target | target | target | target | target | target | target | target | target |
timestamp | ||||||||||
2019-01-01 | 3390640.0 | 2388210.0 | 5470110.0 | 3944450.0 | 3501260.0 | 4130470.0 | 3674330.0 | 3485820.0 | 6132930.0 | 3601470.0 |
2019-02-01 | 2669090.0 | 2027720.0 | 4399290.0 | 2872760.0 | 2747370.0 | 3117610.0 | 2978520.0 | 2791070.0 | 4419260.0 | 3253810.0 |
2019-03-01 | 2865420.0 | 2308770.0 | 5043260.0 | 3898760.0 | 3285700.0 | 3581490.0 | 3377230.0 | 3065430.0 | 4946960.0 | 3913890.0 |
2019-04-01 | 3269360.0 | 2453530.0 | 5564900.0 | 4667040.0 | 3600160.0 | 3915510.0 | 3419200.0 | 3495010.0 | 5908400.0 | 4485810.0 |
2019-05-01 | 3305000.0 | 2373230.0 | 5422920.0 | 4632980.0 | 3300240.0 | 3826670.0 | 3127490.0 | 3504710.0 | 5820050.0 | 4034560.0 |
2019-06-01 | 2982130.0 | 2216290.0 | 4795610.0 | 4308310.0 | 2992590.0 | 3523600.0 | 2742860.0 | 3099660.0 | 5032720.0 | 3483950.0 |
2019-07-01 | 2820750.0 | 1885370.0 | 3799020.0 | 3630070.0 | 2806560.0 | 2936160.0 | 2334890.0 | 2714280.0 | 4224710.0 | 3352260.0 |
2019-08-01 | 3138740.0 | 1915560.0 | 3954220.0 | 3991890.0 | 3020100.0 | 3173640.0 | 2612630.0 | 3022230.0 | 4493300.0 | 3752480.0 |
2019-09-01 | 3868920.0 | 2427400.0 | 5586130.0 | 5218370.0 | 3916480.0 | 4296840.0 | 3525440.0 | 4044120.0 | 6217680.0 | 5196860.0 |
2019-10-01 | 4880810.0 | 2878200.0 | 6980350.0 | 5968900.0 | 5265910.0 | 5155930.0 | 4458310.0 | 5258060.0 | 7904450.0 | 6638270.0 |
2019-11-01 | 5525890.0 | 3280960.0 | 7617220.0 | 7008110.0 | 6389720.0 | 5806650.0 | 5154010.0 | 6053910.0 | 8480420.0 | 6654180.0 |
2019-12-01 | 4830620.0 | 3216540.0 | 6905310.0 | 6946570.0 | 5821860.0 | 5487170.0 | 4596640.0 | 5446120.0 | 7458580.0 | 5474160.0 |
# вычисляем метрику MAPE
mape = MAPE()
mape_=mape(y_true=test_ts, y_pred=forecast_ts)
mape_ = sorted(mape_.items(), key=lambda item: item[1], reverse=True)
mape_
[('Sh8', 28.17434156973828), ('Sh6', 19.517295260293512), ('Sh2', 13.288054471502697), ('Sh7', 13.25201785225388), ('Sh9', 10.819537264180356), ('Sh5', 10.509905475866542), ('Sh3', 8.591273575294613), ('Sh4', 8.255704643405032), ('Sh0', 5.609373332921791), ('Sh1', 2.419442811335598)]
Выведем графики
plot_forecast(forecast_ts, test_ts,
train_ts, n_train_samples=12)
Таким образом в качестве прогноза мы получили прошлогодние месячные выручки, а в качестве теста - фактические значения выручки в прогнозируемом году. А точность прогноза по каждому магазину показывает изменение годовой выручки в этом году относительно прошлого и этот показатель можно использовать в работе. Так как это абсолютная ошибка, нельзя понять лучше стал работать магазин или хуже. Чтобы ответить на этот вопрос лучше вместо MAPE использовать показатель MPE (
mean forecast error). Для того, чтобы его рассчитать проведем некоторые преобразования : помощью значения параметра flatten=True метода .to_pandas()
превращаем объект TSDataset в датафрейм с плоским индексом
test_df = test_ts.to_pandas(flatten=True)
test_df.head()
imestamp | segment | target | |
---|---|---|---|
0 | 2019-01-01 | Sh0 | 3337850 |
1 | 2019-02-01 | Sh0 | 2603750 |
2 | 2019-03-01 | Sh0 | 2774420 |
3 | 2019-04-01 | Sh0 | 3146650 |
4 | 2019-05-01 | Sh0 | 3164340 |
forecast_df = forecast_ts.to_pandas(flatten=True)
forecast_df.head()
timestamp | segment | target | |
---|---|---|---|
0 | 2019-01-01 | Sh0 | 3390640.0 |
1 | 2019-02-01 | Sh0 | 2669090.0 |
2 | 2019-03-01 | Sh0 | 2865420.0 |
3 | 2019-04-01 | Sh0 | 3269360.0 |
4 | 2019-05-01 | Sh0 | 3305000.0 |
Объединим два датафрейма и переименуем колонки для ясности
df_test_fc=pd.merge(test_df, forecast_df, how='inner', on=['timestamp','segment'])
df_test_fc=df_test_fc.rename(columns={"timestamp":"month", "segment":"shop", "target_x":"rev_test", "target_y":"rev_fc"})
df_test_fc.head()
month | shop | rev_test | rev_fc | |
---|---|---|---|---|
0 | 2019-01-01 | Sh0 | 3337850 | 3390640.0 |
1 | 2019-02-01 | Sh0 | 2603750 | 2669090.0 |
2 | 2019-03-01 | Sh0 | 2774420 | 2865420.0 |
3 | 2019-04-01 | Sh0 | 3146650 | 3269360.0 |
4 | 2019-05-01 | Sh0 | 3164340 | 3305000.0 |
А теперь рассчитаем обе метрики
df_test_fc['MPE']=(df_test_fc['rev_test']-df_test_fc['rev_fc'])/df_test_fc['rev_test']
df_test_fc['MAPE']=abs(df_test_fc['rev_test']-df_test_fc['rev_fc'])/df_test_fc['rev_test']
df_test_fc
month | shop | rev_test | rev_fc | MPE | MAPE | |
---|---|---|---|---|---|---|
0 | 2019-01-01 | Sh0 | 3337850 | 3390640.0 | -0.015816 | 0.015816 |
1 | 2019-02-01 | Sh0 | 2603750 | 2669090.0 | -0.025095 | 0.025095 |
2 | 2019-03-01 | Sh0 | 2774420 | 2865420.0 | -0.032800 | 0.032800 |
3 | 2019-04-01 | Sh0 | 3146650 | 3269360.0 | -0.038997 | 0.038997 |
4 | 2019-05-01 | Sh0 | 3164340 | 3305000.0 | -0.044452 | 0.044452 |
... | ... | ... | ... | ... | ... | ... |
115 | 2019-08-01 | Sh9 | 3414480 | 3752480.0 | -0.098990 | 0.098990 |
116 | 2019-09-01 | Sh9 | 4558770 | 5196860.0 | -0.139970 | 0.139970 |
117 | 2019-10-01 | Sh9 | 5633360 | 6638270.0 | -0.178386 | 0.178386 |
118 | 2019-11-01 | Sh9 | 5456970 | 6654180.0 | -0.219391 | 0.219391 |
119 | 2019-12-01 | Sh9 | 4338710 | 5474160.0 | -0.261702 | 0.261702 |
120 rows × 6 columns
Таким образом мы видим как соотносятся по магазинам текущие выручки м прошлогодние. Если сгруппируем, то получим среднегодовой показатель
df_er=df_test_fc.groupby('shop').agg(
avg_mpe=('MPE', 'mean'),
avg_mape=('MAPE', 'mean')).reset_index()
df_er
shop | avg_mpe | avg_mape | |
---|---|---|---|
0 | Sh0 | -0.056094 | 0.056094 |
1 | Sh1 | 0.024194 | 0.024194 |
2 | Sh2 | -0.132881 | 0.132881 |
3 | Sh3 | -0.048343 | 0.085913 |
4 | Sh4 | 0.082557 | 0.082557 |
5 | Sh5 | -0.105099 | 0.105099 |
6 | Sh6 | -0.195173 | 0.195173 |
7 | Sh7 | 0.132520 | 0.132520 |
8 | Sh8 | -0.281743 | 0.281743 |
9 | Sh9 | -0.055003 | 0.108195 |
В заключении можно сказать что и наивный прогноз можно с пользой использовать.
Комментариев нет:
Отправить комментарий