原创：张春阳
用这些模型扩展你的时间序列库
Regularizing, Bagging, Stacking
时间序列数据通常有四个组成部分：
- Autoregression
- Seasonality
- Trend
- Residual
如果能够预测这些成分，你几乎可以预测任何时间序列。听起来很简单，对吧？
但是不完全是这样。关于指定模型的最佳方法有很多不好确定的地方，为了能够更好地解释这些元素，过去几年中已经发布了许多研究来寻找最佳方法，最先进的模型，其中以递归和其他神经网络模型占据了中心地位。除此之外，许多系列还有其他需要考虑的影响，比如假期和规律性的休息。
总而言之，预测任何给定的序列通常不像使用线性回归那么容易。非线性估计和集成方法可以与线性方法相结合，以找到任何序列的最佳方法。在这篇文章中，我从 Scitkit 学习库中拿出了这些模型的示例，以及如何利用它们来最大限度地提高准确性。这些方法应用于每日访客数据集，可在 Kaggle 或 RetressionIt 上找到 2000 多个观测值。
这篇博文的结构相当重复，每个应用模型都有相同的图表和评估指标。如果您已经熟悉了 Scikit 中的机器学习模型和 API，比如 MLR、XGBoost 等，那么可以直接跳到 Bagging 和 Stacking 部分来看它们之间的差异性的对比。你可以在 GitHub 上找到完整的笔记本。
准备模型
所有模型都使用 scalecast 软件包运行，该软件包包含结果，并围绕时间序列数据包装 Scikit learn 和其他模型。默认情况下，其所有预测都使用动态多步预测，与平均一步预测相比，它在整个预测期内返回更合理的准确性/误差指标。
pip install scalecast
我们将使用60天的预测期，也将使用60天集来验证每个模型并调整其超参数。所有模型将在20%的原始数据上进行测试：
f=Forecaster(y=data['First.Time.Visits'],current_dates=data['Date'])
f.generate_future_dates(60)
f.set_test_length(.2)
f.set_validation_length(60)
从 EDA（此处未显示）来看，过去4周内似乎存在自相关，前7个因变量滞后可能是显著的。
f.add_ar_terms(7) # 7 auto-regressive terms
f.add_AR_terms((4,7)) # 4 seasonal terms spaced 7 apart
对于该数据集，必须考虑几个季节性因素，包括每日、每周、每月和季度波动。
f.add_seasonal_regressors(
'month',
'quarter',
'week',
'dayofyear',
raw=False,
sincos=True
) # fourier transformation
f.add_seasonal_regressors(
'dayofweek',
'is_leap_year',
'week',
raw=False,
dummy=True,
drop_first=True
) # dummy vars
最后，我们可以通过添加年份变量来模拟该系列的趋势：
f.add_seasonal_regressors('year')
对于所有这些模型，通常需要向它们提供平稳的时间序列数据。我们可以通过增强的Dickey Fuller测试确认该数据是平稳的：
MLR
我们将从简单开始，尝试应用多元线性回归（MLR）。该模型快速、简单，无需调整超参数。它通常具有极高的准确性，即使使用更先进的方法也很难击败它。
它假设模型中的所有组件可以以线性方式组合，以预测最终输出：
上图公式中，j 是增加的回归数（在我们的例子中，AR、季节和趋势分量），α 是相应的系数。在我们的代码中，调用此函数如下所示：
f.set_estimator('mlr')
f.manual_forecast()
值得注意的是，时间序列的线性方法更常见的模型是 ARIMA，它也使用序列的误差作为回归。MLR假设序列的误差是不相关的，这在时间序列中是虚假的。也就是说，在我们的分析中，使用 MLR 获得了13%的测试集平均绝对百分比误差和76%的 R2。让我们看看是否可以通过增加复杂性来克服这一点。
Lasso
接下来回顾的三个模型，Lasso、Ridge 和 ElasticNet，都使用 MLR 的相同基础函数进行预测，但系数的估计方式不同。有 L1 和 L2 正则化参数，用于减小系数的大小，从而减少过度拟合，并可导致更好的异常值的预测。在我们的案例中，这可能是一种很好的尝试方法，因为 MLR 的样本内 R2 得分为 95%，显著大于样本外 R2 的 76%，表明过度拟合。
有一个参数可以用套 Lasso 估计——L1 惩罚参数的大小，或alpha。我们可以通过对验证数据集进行100个 alpha 值的网格搜索来实现这一点。看起来是这样的：
f.add_sklearn_estimator(Lasso,'lasso')
f.set_estimator('lasso')
lasso_grid = {'alpha':np.linspace(0,2,100)}
f.ingest_grid(lasso_grid)
f.tune()
f.auto_forecast()
Lasso（以及 Ridge 和 ElasticNet）使用缩放输入，以便惩罚参数平衡到所有系数，这一点很重要。默认情况下，Scalecast 使用 MinMax 缩放器。
选择的最佳 α 值为 0.081。该模型既没有提高MLR模型的样本外精度，也没有减少过拟合。我们可以用 Ridge 模型再试一次。
Ridge
Ridge 与 Lasso 类似，只是它使用 L2 惩罚。这里的区别在于，L1惩罚可以将某些系数减少到零，而Ridge 只能将系数减少到接近零。通常，这两种模型都会产生相似的结果。我们可以使用为 Lasso 模型创建的网格来调 Ridge 模型。
f.set_estimator('ridge')
f.ingest_grid(lasso_grid)
f.tune()
f.auto_forecast()
Ridge 的最佳 α 选择为0.384，其结果与 Lasso 模型相似。我们还有一个可以将正则化应用于 MLR 的模型：ElasticNet。
ElasticNet
Scikit learn 提供的 ElasticNet 模型将使用线性模型预测输出，但现在它将混合 L1 和 L2 惩罚。现在需要调整的关键参数包括：
- L1/L2 惩罚系数 (
l1_ratio) - 惩罚值 (
alpha)
Scalecast 为 ElasticNet 提供了一个很好的默认验证网格，所以我们不必创建一个。
f.set_estimator('elasticnet')
f.tune()
f.auto_forecast()
ElasticNet 的最佳参数是 l1_ration 1（相当于套索模型），alpha 为 0.3。它的性能与 Lasso 和 Ridge 模型相当，并且再次没有减少过度拟合。
Random Forest
在用尽了几种线性方法来估计这个序列之后，我们可以转向非线性方法。其中最流行的是随机森林，这是一种基于树的集成方法。它的功能是用原始数据集的自举样本（带替换的采样）聚合多个树估计量。每个样本可以利用原始数据集的不同行和列，最终结果是应用于每个样本的指定数量的底层决策树估计数的平均值。
对于使用随机森林进行时间序列预测，人们不时提出一些批评，例如估计值不能大于最大观测值或小于最小观测值。有时，这个问题可以通过确保平稳数据只提供给估计器来避免。但一般来说，这个模型并不以其时间序列的威力而闻名。让我们看看在这个示例中它是如何工作的。我们可以指定自己的验证网格。
rf_grid= {
'max_depth':[2,3,4,5],
'n_estimators':[100,200,500],
'max_features':['auto','sqrt','log2'],
'max_samples':[.75,.9,1],
}
然后运行预测。
f.set_estimator('rf')
f.ingest_grid(rf_grid)
f.tune()
f.auto_forecast()
从现在开始，我必须粘贴两个表格用于模型基准测试，否则就无法阅读。下表的行与上表中列出的模型相对应。
不幸的是，我们在 Random Forest 上运气不太好。它的表现比之前评估的要差得多。然而，一旦我们用 BaggingRegressionor 模型概述预测，这里介绍的 bootstrapped sampling 的概念将是完整的。
XGBoost
XGBoost 是一个很难简单解释的模型，所以如果读者感兴趣的话，我会按照这里（对于初学者）和这里（对于深入研究）的文章来做。其基本思想是，与随机森林类似，通过一系列决策树进行估计，最终结果是每个基础估计的加权平均数。与随机森林不同的是，树是按顺序建造的，每个后续的树都对之前的树的残差进行建模，希望尽可能减少最终的残差。抽样不是随机的，而是基于前一棵树的弱点。通过这种方式，结果是通过增加样本来获得的，而随机森林使用 bootstrapped aggregation，其中每个树和样本独立于其他树和样本。这是对数据分布下真正发生的事情的过分简化，因此，如果这个解释不能让你满意，请阅读上述 XGBoost 的文章。在这样的模型中有许多超参数需要调整，我们可以构建如下网格：
xgboost_grid = {
'n_estimators':[150,200,250],
'scale_pos_weight':[5,10],
'learning_rate':[0.1,0.2],
'gamma':[0,3,5],
'subsample':[0.8,0.9]
}
评估模型：
f.set_estimator('xgboost')
f.ingest_grid(xgboost_grid)
f.tune()
f.auto_forecast()
XGBoost 在测试集上的 MAPE 结果超过 MLR 12% 。然而，更加过分的是，样本中的 R2 分数接近100%。
到目前为止，希望您了解这些模型是如何构建和评估的。在提供的 notebook 中，您还可以看到LightGBM（微软的增强树模型，类似于XGBoost）、随机梯度下降和K近邻的示例。为了简洁起见，我将跳到 Baggign 和 Stacking 模型来完成这篇文章。
Bagging
来自Scikit learn的 BaggingRegressor 使用了本文随机森林部分中介绍的 Same Sampling 的概念，但我们可以使用任何我们喜欢的模型，而不是每个基础估计量都是一棵决策树。在本例中，我指定了10个多层感知器神经网络模型，每个模型有三层，每层100个单元，以及LBFGS解算器。每个数据子集都可以使用原始数据集大小的90%对观测值进行采样，还可以随机使用原始数据集50%的特征。在代码中，它如下所示：
f.add_sklearn_estimator(BaggingRegressor,'bagging')
f.set_estimator('bagging')
f.manual_forecast(
base_estimator= MLPRegressor(
hidden_layer_sizes=(100,100,100),
solver='lbfgs'
),
max_samples= 0.9,
max_features= 0.5,
)
该模型迄今为止表现最好，测试集 MAPE 为 11%，测试集 R2 得分为 79%。在XGBoost和MLR这两个相同的指标下，这比下一个最好的模型要好得多。
Stacking
最后一个模型是StackingRegressionor。它创建了一个新的估计器，该估计器使用来自其他指定模型的预测来创建最终预测。它通过传递给 final_estimator 参数的估计器进行最终预测。
在我们的例子中，我们将对该数据集上评估最佳样本的四个模型进行叠加：
- K-nearest Neighbors (KNN)
- XGBoost
- LightGBM
- Stochastic Gradient Descent (SGD)
最终估计量是我们在上一节中定义的 bagging 估计量。
使用 BaggingRegressionor 作为最终估计器的优势在于，即使 KNN 和 XGBoost 模型是高度过拟合的，MLP 模型在概念上也应该相应地人为地将其预测加权，因为该模型仅使用输入任何给定MLP模型的一半特征进行训练，它还应该学会相信两个模型的预测，这两个模型并没有过度拟合：LightGBM 和 SGD。因此，对模型的评估如下：
f.add_sklearn_estimator(StackingRegressor,'stacking')
f.set_estimator('stacking')
results = f.export('model_summaries')
estimators = [
('knn',
KNeighborsRegressor(**results.loc[results['ModelNickname'] == 'knn','HyperParams'].values[0])),
('xgboost',
XGBRegressor(**results.loc[results['ModelNickname'] == 'xgboost','HyperParams'].values[0])),
('lightgbm',
LGBMRegressor(**results.loc[results['ModelNickname'] == 'lightgbm','HyperParams'].values[0])),
('sgd',
SGDRegressor(**results.loc[results['ModelNickname'] == 'sgd','HyperParams'].values[0])),
]
final_estimator = BaggingRegressor(
base_estimator = MLPRegressor(
hidden_layer_sizes=(100,100,100),
solver='lbfgs',
),
max_samples = 0.9,
max_features = 0.5,
)
f.manual_forecast(estimators=estimators,final_estimator=final_estimator)
增加模型的复杂性并不总是能改善结果，但这似乎就是这里发生的事情。我们的叠加模型明显优于其他模型，测试集MAPE为10%，测试集R2得分为83%。它还具有与评估的MLR相当的样本内指标。
Backtesting
为了进一步验证我们的模型，我们可以对它们进行反向测试。这是一个迭代评估其在过去n个预测期内的准确性的过程，以查看如果实施该模型，只根据每个预测期之前的观测值进行训练，实际会取得什么结果。默认情况下，scalecast选择10个预测范围，其长度由对象中生成的未来日期的数量决定（在本例中为60）。下面是它在代码中的样子：
f.backtest('stacking')
f.export_backtest_metrics('stacking')
这告诉我们，平均而言，使用60天的预测长度，我们的模型将获得 7% 的 MAPE 分数和 76% 的R2。这是一个很好的方法，可以验证模型不仅仅是幸运地得到了特定的测试集，而且实际上可以推广到看不见的数据。
总结
这篇文章介绍了几种预测的建模概念，包括线性方法（MLR、套索、岭）、树集合（随机森林、XGBoost、LightGBM）、Bagging 和 Stacking。通过增加非线性方法的复杂度，并使用巧妙的抽样方法，我们在一个包含大约 400 个观察值的测试集上获得了 10% 的测试集 MAPE 分数。
