基于树模型的时间序列预测实战
现在，我们将了解一个与经典ARIMA时间序列建模不同的新领域。在监督学习模型中，仅仅使用单变量时间序列似乎信息有限，预测也比较困难。因此，为了生成足够的特征，我们需要采取一些方法，例如创建大量的滞后变量。此外，关于预测目标值，我们需要用过去的项来预测未来的项，而且需要决定是一步领先还是多步领先。
在这篇文章中，云朵君将和大家一起学习以下内容：
- 从单变量时间序列中创建特征，
- 使用提前一步预测的监督学习框架，
- 建立轻型 GBM 预测模型，并提供模型的可解释性。
以上内容并不仅限于lightGBM，但该过程也适用于其他回归器，如随机森林 (RF)、梯度提升机 (GBM) 或极端梯度提升 (XGB)。感兴趣的伙伴可以自己尝试。
从单变量时间序列中创建特征
在单变量时间序列中，我们只能获得有限的信息。ARIMA 模型使用过去的值来预测未来的值，因此过去的值是重要的候选特征，可以创建许多滞后回归因子。时间指数是一个有价值的领域，我们可以基于此创建特征。由于日历上的事件和年度事件在我们的生活中不断重复，它们为我们的过去留下了印记，为我们的未来提供了教益。因此，我们可以从与时间相关的特征入手。
创建基于时间的特征
创建基于时间的特征，包括日期、星期、季度等各种特征，通过 pandas series 的 "date" 类中提供的一系列函数，我们可以轻松实现这些需求。
def create_date_features(df):
df['month'] = df.date.dt.month
df['day_of_month'] = df.date.dt.day
df['day_of_year'] = df.date.dt.dayofyear
df['week_of_year'] = df.date.dt.weekofyear
df['day_of_week'] = df.date.dt.dayofweek + 1
df['year'] = df.date.dt.year
df['quarter'] = df.date.dt.quarter
df['hour_of_day'] = df.date.dt.hour
df['weekday'] = df.date.dt.weekday
df['is_year_start'] = df.date.dt.is_year_start.astype(int)
df['is_year_end'] = df.date.dt.is_year_end.astype(int)
df['is_month_start'] = df.date.dt.is_month_start.astype(int)
df['is_month_end'] = df.date.dt.is_month_end.astype(int)
df['is_quarter_start'] = df.date.dt.is_quarter_start.astype(int)
df['is_quarter_end'] = df.date.dt.is_quarter_end.astype(int)
df['is_quarter_end'] = df.date.dt.is_quarter_end.astype(int)
return df
df = create_date_features(df)
df.head()%matplotlib inline
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
# Use google colab
from google.colab import drive
drive.mount('/content/gdrive')
# Use the electric consumption dataset
path = '/content/gdrive/My Drive/data/time_series'
data = pd.read_csv(path + '/electric_consumption.csv')
data.tail()
这里我们使用数据集 （获取方法：公众号：数据STUDIO 联系 云朵君 即可）包含 2017 年至 2020 年的每小时用电量记录。除 Date（日期 ）和 Consumption（用电量）字段外，它还包含其他字段：
data.columns
Index(['Date', 'Homestead_maxtempC', 'Homestead_mintempC',
'Homestead_DewPointC', 'Homestead_FeelsLikeC', 'Homestead_HeatIndexC',
'Homestead_WindChillC', 'Homestead_WindGustKmph',
'Homestead_cloudcover', 'Homestead_humidity', 'Homestead_precipMM',
'Homestead_pressure', 'Homestead_tempC', 'Homestead_visibility',
'Homestead_winddirDegree', 'Homestead_windspeedKmph', 'Consumption'],
dtype='object')
我们将只使用 Date（日期 ）和 Consumption（用电量）字段来处理单变量时间序列。
df = data[['Date','Consumption']].copy() # 只复制两个字段
df.columns = ['date','y'] # 将目标重命名为 "y"，以方便编码
df["date"] = pd.to_datetime(df["date"]) # 转换为 pandas 的日期格式
df = df.sort_values(by='date') # 确保按日期排序
df.head()
可视化
plt.figure(figsize=(10,4))
plt.plot(df['date'], df["y"])
plt.xlabel("Date")
plt.ylabel("Consumption")
plt.show()
显而易见的周期性模式。
应用函数来创建日期特征：
df = create_date_features(df)
df.head()
需要进行一个步骤。在我们的模型中，有几个字段不应作为数字特征，而应作为分类特征。我们将把它们转化为虚拟变量
to_dummy = ['weekday', 'month', 'quarter', 'year', 'day_of_month', 'week_of_year', 'day_of_week', 'hour_of_day']
df = pd.get_dummies(df, columns= to_dummy)
我们已经创建了一个特征列表，现在学习如何创建滞后变量。
创建滞后特征和未来特征
在自动回归模型中，回归变量是滞后值。可以使用 .shift(n) 来创建滞后特征。接下来，我将在数据集 ff 中创建三个滞后特征。
ff = df.copy()
ff['y-1'] = ff['y'].shift(1)
ff['y-2'] = ff['y'].shift(2)
ff['y-3'] = ff['y'].shift(3)
ff.head()
编写一个 forloop 来创建多个滞后特征。下面我将在不同的数据集 ff 中创建 5 个滞后变量：
ff = df.copy()
def create_lagged(df, n_vars):
# Use a forloop
for i in range(n_vars):
# The name will be y-1, y-2, etc.
name = ('y-%d' % (i+1))
df[name] = df['y'].shift(i+1)
return df
ff = create_lagged(ff, 5)
ff.head()
显然，我们也可以将数值前移，使其成为未来的目标值，如以下所示。
ff = df.copy()
ff['y+1'] = ff['y'].shift(-1)
ff['y+2'] = ff['y'].shift(-2)
ff['y+3'] = ff['y'].shift(-3)
ff.tail()
正式为建模数据 df 创建 25 个滞后变量。
df = create_lagged(df, 25)
df.columns
数据集中包含了['date', 'y', 'y-1', ..., 'y-25']的数据。在此基础上，可以进行一系列汇总统计，如过去 n 小时、n 天或 n 周的总和或平均值。
创建移动平均值
另外，可以创建6、9、12、18、21和24小时的移动平均值。
def roll_mean_features(df, windows):
df = df.copy()
for window in windows:
df['mv_' + str(window)] = df['y'].transform(
lambda x: x.shift(1).rolling(window=window, min_periods=10, win_type="triang").mean())
return df
df = roll_mean_features(df, [6, 9, 12, 15, 18, 21, 24])
df.tail()
监督学习框架用于提前预测。模型目标是 y，特征包括滞后项 y-1到y-25以及时间相关和移动平均变量。
该模型可以通过yt-1到yt-25产生下一期的yt，即提前一步预测。在现实应用中，多步预测也很常见，传统方法是建立n个模型来预测接下来的n期。
建立 LightGBM 预测模型
LightGBM是微软开发的梯度提升框架，它使用叶向树生长以提高准确性。相比之下，level-wise树会尝试在同一级别的分支上生长，看起来更平衡。由于其能够处理大型数据集和并行化训练，因此比其他提升算法更高效、更快速，同时内存占用更低。此外，它原生支持分类特征，无需进行单次编码。
梯度提升模型是机器学习算法的一种，它将多个较弱的模型组合在一起，从而创建一个强大的预测模型。它的基本思想是迭代训练决策树，每棵树都试图纠正前一棵树所犯的错误。最终的预测结果是所有决策树预测结果的总和。梯度提升模型特别适用于处理复杂的数据集，可以处理大量特征和特征之间的交互，并且对过度拟合也很稳健，同时能够处理缺失值。常用的算法有梯度提升机（GBM）、XGB 和 LightGBM。
首先，我们将数据分为 "实时" 数据和 "非实时" 数据。
划分训练和测试集
下面的代码将时间序列切割成 "实时" 数据作为训练数据，"非实时" 数据作为测试数据。
# Count the days
num_days = (df['date'].max() - df['date'].min()).days
# reserve 20% for out-of-time
oot = num_days*0.2
# Get the cutdate
cutdate = df['date'].max() - timedelta(days = oot)
# Create the training data
train = df.loc[(df['date'] <= cutdate), :]
print("Training data: from", train['date'].min(), "to", train['date'].max() )
# Create the test data
test = df.loc[(df['date'] > cutdate), :]
print("Test data: from", test['date'].min(), "to", test['date'].max() )
LightGBM 建模
LightGBM 有许多超参数可以调整。在 Python 字典中指定关键超参数。
import lightgbm as lgb
lgb_params = {# 平均绝对误差
'metric'： {'mae'}、
# 树中树叶的数量
'num_leaves'： 10,
# 学习日期
learning_rate': 0.02、
# 随机选取 80% 的特征到
'feature_fraction'： 0.8,
# 树的最大深度
'max_depth'：5,
# 忽略训练进度（不显示任何内容）
'verbose': 0、
# 提升迭代次数
num_boost_round'： 15000,
# 如果精度没有提高，就停止训练
early_stopping_rounds': 200、
# 使用计算机上的所有内核
nthread'： -1}
LightGBM 有一个".Dataset()"代码类，用于打包目标变量、回归变量和数据。如下所示，操作非常简单。
# Remove unnecessary fields
cols = [col for col in train.columns if col not in ['date', 'y', "y-1", "y-2", "y-3", "y+1", "y+2", "y+3"]]
# 为训练数据和测试数据创建 y 和 X
Y_train = train['y']
X_train = train[cols]
Y_test = test['y']
X_test = test[cols]
# Use the Dataset class of lightGBM
lgbtrain = lgb.Dataset(data=X_train, label=Y_train, feature_name=cols)
lgbtest = lgb.Dataset(data=X_test, label=Y_test, reference=lgbtrain, feature_name=cols)
model = lgb.train(lgb_params, lgbtrain,
valid_sets=[lgbtrain, lgbtest],
num_boost_round=lgb_params['num_boost_round']
)
预测准确性评估
完成后，我们可以得出训练数据和测试数据的预测值，并评估预测准确度。我们将使用标准指标平均绝对百分比误差 (MAPE) 来评估预测准确度。MAPE 是绝对百分比误差的平均值，10% 的 MAPE 意味着预测值和实际值之间的平均偏差为 10%。
from sklearn.metrics import mean_absolute_percentage_error
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)
mean_absolute_percentage_error(Y_test, y_pred_test)
MAPE 为可接受的 1.90%，这意味着预测值可以高于或低于实际值 1.90%。
实际值与预测值可视化
# 将预测值添加到训练期
train_pred = train.copy()
train_pred['y_pred_train'] = y_pred_train
# 将预测值添加到测试期
test_pred = test.copy()
test_pred['y_pred_test'] = y_pred_test
print([train_pred.shape, test_pred.shape])
# 合并训练期和测试期
actual_pred = pd.concat([train_pred, test_pred], axis=0)
actual_pred.shape
# 用蓝色绘制实际值
# 用橙色绘制训练期的预测值
# 用绿色标出测试期的预测值
plt.figure(figsize=(10,4))
plt.plot(actual_pred['date'], actual_pred[["y",'y_pred_train','y_pred_test']])
plt.xlabel("Date")
plt.ylabel("Actual vs. Predictions")
plt.show()
橙色线是训练期的预测值，绿色线是测试期的预测值。这两条线与实际值非常吻合。
模型可解释性
基于树的模型的优势之一是其可视性。可以通过变量重要性图直观地看到特征对预测的影响。
lgb.plot_importance(model,
max_num_features=20,
figsize=(10, 10),
importance_type="gain")
plt.show()
特征重要性图显示，影响最大的三个变量是 、和。 的出现在意料之中，因为耗电量通常遵循 AR(1) 模式，不会突然变化。12 小时移动平均线 mv_12 出现了，但并不显著。
结论
在本章中，我们探讨了单变量时间序列特征的创建方法，以及如何将其纳入基于树的监督学习框架中。我们利用 lightGBM 模型进行了一步预测，并展示了如何利用变量显著图提高模型可解释性。希望对你有帮助，不妨点个赞❤️
🏴☠️宝藏级🏴☠️ 原创公众号『数据STUDIO』内容超级硬核。公众号以Python为核心语言，垂直于数据科学领域，包括可戳👉 Python｜MySQL｜数据分析｜数据可视化｜机器学习与数据挖掘｜爬虫 等，从入门到进阶！
长按👇关注- 数据STUDIO -设为星标，干货速递