基于树模型的时间序列预测实战

数据STUDIO

发布于 2024-04-18 15:23:29

2260

发布于 2024-04-18 15:23:29

文章被收录于专栏：数据STUDIO

你有没有发现，近期大模型非常火热，但在大多数时间序列预测竞赛中，排名靠前的模型竞赛主要采用基于树的机器学习方法，特别是像LightGBM这样的梯度提升模型表现出色。因此，我将在本文中介绍基于树的时间序列预测。尽管有Prophet和NeuralProphet等方便的工具，但是了解基于树的模型仍然具有很高的价值。

现在，我们将了解一个与经典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以及时间相关和移动平均变量。

y = fn(y_t-1, y_t-2, ..., y_t-25, ... )

该模型可以通过yt-1到yt-25产生下一期的yt，即提前一步预测。在现实应用中，多步预测也很常见，传统方法是建立n个模型来预测接下来的n期。

y+1 = f_{n_1}(y_{t-1}，y_{t-2}，...，y_{t-25}，...)

y+n = fn_n(y_{t-1}，y_{t-2}，...，y_{t-25}，...)

建立 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()

特征重要性图显示，影响最大的三个变量是

y_1、y_24 和 y_23

。

y-1

的出现在意料之中，因为耗电量通常遵循 AR(1) 模式，不会突然变化。12 小时移动平均线 mv_12 出现了，但并不显著。