教程1：线性回归
target = weight_1 * feature_1 + weight_2 * feature_2 + bias时间步(Time-step)特征
Hardcover = weight * time + bias把这里的时间称为：时间虚拟变量(time dummy)，因为这是假的时间。
滞后(lag)特征
Hardcover = weight * lag_1 + bias这里是滞后一步，这里的线性回归就是根据前一个值预测现在的值。
时间步(time-step)特征和滞后(lag)特征的可以组合在一起。
教程2：趋势(Trend)
时间序列的趋势部分表示该序列平均值的持续、长期变化。 趋势是一系列中移动最慢的部分，代表了最大时间尺度的重要性。
一个序列中任何持续的和缓慢移动的变化都可能构成一个趋势——例如，时间序列通常具有变化的趋势。
移动平均图(moving average plot)
要查看时间序列可能具有什么样的趋势，我们可以使用移动平均图。 为了计算时间序列的移动平均值，我们计算某个定义的宽度的滑动窗口内的值的平均值。 图表上的每个点代表位于任一侧窗口内的系列中所有值的平均值。 这个想法是为了消除序列中的任何短期波动，以便只保留长期的变化。
线性趋势的移动平均图。 曲线上的每个点（蓝色）是大小为 12 的窗口内的点（红色）的平均值。
上面的 Mauna Loa 系列如何年复一年地重复上下运动——这是一种短期的季节性变化。 要使变化成为趋势的一部分，它应该比任何季节性的变化发生的时间更长。 因此，为了可视化趋势，我们在比该序列中任何的季节性的周期更长的时间段内取平均值。 对于 Mauna Loa 系列，我们选择了大小为 12 的窗口来平滑每年的季节。
线性回归
多项式回归
预测未来30天的趋势:
趋势模型之所以有用，有很多原因。 除了作为更复杂模型的基线或起点之外，我们还可以将它们用作有些无法学习趋势的算法（如 XGBoost 和随机森林）的“混合模型”中的一个组件。
sklearn：LinearRegression
教程 3-季节性(Seasonality)
什么是季节性？
只要序列的平均值有规律的、周期性的变化，我们就说时间序列表现出季节性。 季节性变化通常遵循时钟和日历——一般一天、一周或一年的重复。 季节性通常是由自然界在几天和几年内的循环或围绕的日期和时间的社会行为惯例的循环驱动的。
四个时间序列的季节模式
两种关于季节性的特征。 第一种，指示器(indicators)，最适合一个季节性周期中有少量的观察值，例如在每天的观察值中找到以周为周期的季节性。 第二种，傅里叶特征(Fourier features)，最适合一个季节性周期中有许多的观察值，例如在每天的观察值中找到以年为周期的季节性。
季节性图(Seasonal Plots)和季节性指示器(Seasonal Indicators)
就像我们使用移动平均图来发现系列中的趋势一样，我们可以使用季节性图来发现季节性。
季节性图显示了针对某个常见时期绘制的时间序列片段，该时期是您要观察的“季节”。 该图显示了维基百科关于 三角学(Trigonometry) 的文章的每日浏览量的季节性图：文章的每日浏览量是在一个共同的 每周 期间绘制的。
该序列有明显的以周为周期季节性，工作日较高，周末下降。
季节性指示器(Seasonal indicators)
季节性指示器是表示时间序列水平的季节性差异的二元特征(Seasonal indicators are binary features that represent seasonal differences in the level of a time series)。 如果您将季节性周期视为分类特征并进行 one-hot 编码，则可以得到季节性指示器。
通过对一周中的每一天进行 one-hot 编码，我们得到每周的季节性指示器。 为 三角学(Trigonometry) 系列创建每周每周的季节性指示器将为我们提供六个新的“虚拟”特征。
（如果删除其中一个指示器，线性回归效果会更好；所以我们在下表中选择删除了星期一。）
| Date | Tuesday | Wednesday | Thursday | Friday | Saturday | Sunday |
|---|---|---|---|---|---|---|
| 2016-01-04 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2016-01-05 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2016-01-06 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2016-01-07 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 |
| 2016-01-08 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 |
| 2016-01-09 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 |
| 2016-01-10 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 |
| 2016-01-11 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... | ... |
向训练数据中添加季节性指示器有助于模型识别季节性周期内的平均值：
普通线性回归学习季节中每个时间的平均值。
指示器就像一个开关一样。 在任何时候，这些指示器中最多有一个的值为“1”（开）。 线性回归给星期一学习到一个基准值为2379，然后根据当天哪个指示器是 “开 ” 对值进行调整；其余的指示器由于值是0, 所以值不会被计算。
傅里叶特征(Fourier Features )和周期图
特征更适合有许多观察值的长季节周期，这种情况使用指示器就很不明智（回忆我们之前的以周为周期的指示器，一周七天就会多出来6个特征，如果观察值过多，就会多出很多特征！）。 傅立叶特征不是为每个日期创建一个特征，而是尝试用几个特征来捕捉季节性曲线的整体形状。
让我们看一下 三角学(Trigonometry) 中的年度季节图。 注意各种频率的重复：每年3次长的上下运动，每年52次的短周运动，也许还有其他。
三角学(Trigonometry)序列的年度季节性变化
我们试图用傅里叶特征捕捉一个季节内的这些频率。 这个想法是在我们的训练数据中包含与我们试图建模的季节具有相同频率的周期性曲线。 我们使用的曲线是三角函数正弦和余弦的曲线。
傅里叶特征是成对的正弦和余弦曲线，从最长的季节开始，每个潜在频率都对应一对。 对年度季节性进行建模的傅立叶对将具有频率：每年一次、每年两次、每年三次，依此类推。
年度季节性的前两个傅立叶对。 上: 频率每年一次。 下: 频率每年两次。
如果我们将一组这些正弦/余弦曲线添加到我们的训练数据中，线性回归算法将计算出适合目标序列中季节性分量的权重。 该图说明了线性回归如何使用四个 傅立叶对 来模拟 三角学(Trigonometry)序列中的年度季节性。
上: 四个傅立叶对的曲线，正弦和余弦的总和以及回归系数。 每条曲线模拟不同的频率。 下: 这些曲线之和近似于季节性模式。
请注意，我们只需要八个特征（四个正弦/余弦对）就可以很好地估计年度季节性。 与需要数百个特征（一年中的每一天一个）的季节性指示器方法相比较。 通过仅使用傅立叶特征对季节性的“主效应”进行建模，向训练数据中添加特征更少，这意味着减少了计算时间并降低了过度拟合的风险。
使用周期图选择傅里叶特征
我们实际上应该在我们的特征集中包含多少傅里叶对呢？ 我们可以用周期图来回答这个问题。 周期图告诉您时间序列中频率的强度。 具体来说，周期图的 y 轴上的值为 (a ** 2 + b ** 2) / 2，其中 a 和 b 是该频率下正弦和余弦的系数（如 在上面的 Fourier Components 图中）。
从左到右，周期图在 Quarterly 之后下降，一年四次。 这就是我们选择四对傅立叶对来模拟年度季节的原因。 我们忽略了Weekly频率，因为它使用季节性指示器来建模更好。
计算傅里叶特征（可选）
了解傅里叶特征的计算方式对于使用它们并不是必不可少的，但如果看到细节可以更好的理解它，下面的单元格说明了如何从时间序列的索引中导出一组傅里叶特征。 （不过，我们将在应用程序中使用来自 statsmodels 的库函数。）
import numpy as np
def fourier_features(index, freq, order):
time = np.arange(len(index), dtype=np.float32)
k = 2 * np.pi * (1 / freq) * time
features = {}
for i in range(1, order + 1):
features.update({
f"sin_{freq}_{i}": np.sin(i * k),
f"cos_{freq}_{i}": np.cos(i * k),
})
return pd.DataFrame(features, index=index)
# Compute Fourier features to the 4th order (8 new features) for a
# series y with daily observations and annual seasonality:
#
# fourier_features(y, freq=365.25, order=4)示例 - Tunnel Traffic
来看看一周和一年的季节性图。
看一下周期图：
周期图与上面的季节图一致：每周季节性较强和每年季节性较弱。 我们将用指示器来对每周季节性进行建模，用傅里叶特征对每年的每年季节性。 从右到左，周期图在双月 (6) 和每月 (12) 之间递减，所以让我们使用 10 个傅立叶对。
我们将使用 DeterministicProcess 创建我们的季节性特征，我们在第 2 课中用于创建趋势(Trend)特征的相同方法。 要使用两个季节性时段（每周和每年），我们需要将其中一个实例化为“附加项”：
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
fourier = CalendarFourier(freq="A", order=10) # 10 sin/cos pairs for "A"nnual seasonality
dp = DeterministicProcess(
index=tunnel.index,
constant=True, # dummy feature for bias (y-intercept)
order=1, # trend (order 1 means linear)
seasonal=True, # weekly seasonality (indicators)
additional_terms=[fourier], # annual seasonality (fourier)
drop=True, # drop terms to avoid collinearity
)
X = dp.in_sample() # create features for dates in tunnel.index
创建特征集后，我们就可以拟合模型并进行预测了。 我们将添加一个 90 天的预测，以了解我们的模型如何在训练数据之外进行推断。 这里的代码与前面课程中的代码相同。
y = tunnel["NumVehicles"]
model = LinearRegression(fit_intercept=False)
_ = model.fit(X, y)
y_pred = pd.Series(model.predict(X), index=y.index)
X_fore = dp.out_of_sample(steps=90)
y_fore = pd.Series(model.predict(X_fore), index=X_fore.index)
ax = y.plot(color='0.25', style='.', title="Tunnel Traffic - Seasonal Forecast")
ax = y_pred.plot(ax=ax, label="Seasonal")
ax = y_fore.plot(ax=ax, label="Seasonal Forecast", color='C3')
_ = ax.legend()
在下一课中，我们将学习如何将时间序列本身用作特征。 使用时间序列作为预测的输入可以让我们对序列中经常出现的另一个情况进行建模：周期。
教程 4-时间序列作为特征(Time Series As Features)
什么是序列依赖（Serial Dependence）？
在前面的课程中，我们研究了像*时间依赖(time dependent)属性这种容易建模的时间序列性质 ，也就是说，我们可以直接从时间索引中获得特征。 然而，一些时间序列只能使用序列依赖(serially dependent)*的属性来建模，即 使用目标序列的过去值作为特征。 根据时间来绘图，这些时间序列的结构可能并不明显； 然而，根据过去的值来绘制，结构就会变得清晰——如下图所示。
这两个系列具有序列依赖，但不具有时间依赖。 右边的点的坐标为 (时间t-1的值, 时间t的值).
借助趋势和季节性，我们训练模型以将曲线拟合到上图左侧的图上 – 模型学习了时间依赖性。 本节课的目标是训练模型以将曲线拟合到右侧的图上 – 我们希望它们学习序列依赖。
周期(Cycles)
周期是指在一个时间序列中，某个时间点的值的增长或衰减取决于它之前时间的值，和时间步长本身不一定有关。它是一种特别常见的序列依赖表现方式。 周期行为是可以影响自身或随时间存在持续影响的系统特点。 经济、流行病、动物种群、火山爆发和类似的自然现象经常表现出周期行为。
四个具有周期行为的时间序列。
周期性行为与季节性的区别在于，周期不一定像季节那样依赖于时间。 一个周期中发生的事情与特定的发生日期无关，而更多地与最近发生的事情有关。 与时间的（至少是相对的）独立性意味着周期行为可能比季节性更不规则。
滞后序列(Lagged Series)和滞后图(Lag Plots)
为了调查时间序列中可能存在的序列依赖性（如周期），我们需要创建序列的“滞后”副本。 Lagging 一个时间序列意味着将其值向前移动一个或多个时间步长，或者等效地将其索引中的时间向后移动一个或多个步骤。 在任何一种情况下，结果都是滞后序列中的观察值似乎发生在较晚的时间。
这显示了美国的月度失业率 (y) 及其第一个和第二个滞后序列（分别为 y_lag_1 和 y_lag_2）。 注意滞后序列的值是如何及时向前移动的。
| Date | y | y_lag_1 | y_lag_2 |
|---|---|---|---|
| 1954-07 | 5.8 | NaN | NaN |
| 1954-08 | 6.0 | 5.8 | NaN |
| 1954-09 | 6.1 | 6.0 | 5.8 |
| 1954-10 | 5.7 | 6.1 | 6.0 |
| 1954-11 | 5.3 | 5.7 | 6.1 |
通过滞后时间序列，我们可以使其过去的值与我们试图预测的值同时出现（换句话说，在同一行中）。 这使得滞后序列可用作建模序列依赖的特征。 为了预测美国失业率序列，我们可以使用y_lag_1和y_lag_2作为特征来预测目标y。 这使用前两个月失业率来预测未来失业率。
滞后图(Lag plots)
时间序列的滞后图显示了其值与滞后值的关系。 通过查看滞后图，时间序列中的序列依赖性通常会变得明显。 我们可以从 US Unemployment 的滞后图中看到，当前失业率与过去的失业率之间存在强烈且明显的线性关系。
自相关的美国失业滞后图。
最常用的序列依赖性度量称为自相关，它只是时间序列与其滞后之一的相关性。 US Unemployment 在滞后1时的自相关为0.99，在滞后2时为0.98，依此类推。\
选择滞后
在选择滞后作为特征时，包含每个具有较大自相关的滞后 通常没有用处。 例如，在 US Unemployment 中，滞后2处的自相关可能完全来自滞后1的“衰减”信息 – 只是从上一步结转的相关性。 如果滞后 2 不包含任何新内容，那么如果我们已经有了滞后 1，就没有理由包含它。
偏自相关告诉您滞后与所有先前滞后的相关性 – 也可以说，是滞后贡献的“新”相关性的数量。 绘制偏自相关可以帮助您选择要使用的滞后特征。 在下图中，滞后1到滞后6不在“无相关性”区间（蓝色），因此我们可以选择滞后 1到滞后6作为 US Unemployment 的特征。 （滞后 11 可能是误报。）
美国失业率的部分自相关通过滞后12与95% 的置信区间没有相关性。
像上面这样的图被称为相关图。 相关图适用于滞后特征，本质上就像周期图适用于傅里叶特征。
最后，我们需要注意自相关和偏自相关是线性依赖性的度量。 因为现实世界的时间序列通常有很大的非线性依赖，在选择滞后特征时最好看滞后图（或者使用一些更一般的依赖度量，比如互信息）。 太阳黑子系列具有非线性相关的滞后，我们可能会忽略自相关。
太阳黑子 序列滞后图.
像这样的非线性关系可以转换为线性关系，也可以通过适当的算法学习。
示例 - 流感趋势
Flu Trends 数据集包含 2009 年至 2016 年间数周因流感而就诊医生的记录。我们的目标是预测未来几周的流感病例数。
我们将采取两种方法。 第一种方法是，我们使用滞后特征预测就诊次数。 第二种方法是使用另一组(谷歌趋势捕获的与流感相关的搜索词)时间序列的滞后来预测医生的就诊次数。
我们的流感趋势数据显示了不规则的周期而不是常规的季节性：高峰往往出现在新年前后，但有时更早或更晚，有时更大或更小。 使用滞后特征对这些周期进行建模将使我们的预测能够对不断变化的条件做出动态反应，而不是像季节性特征那样受限于确切的日期和时间。
让我们先看一下滞后和自相关图：
滞后图表明 FluVisits与其滞后的关系主要是线性的，而偏自相关图表明可以使用滞后 1、2、3 和 4 来捕获依赖关系。我们可以在 Pandas 中使用 shift 方法来滞后时间序列。 同时，我们将用 0.0 填充滞后创建的缺失值。
在之前的课程中，我们能够为训练数据之外的任意多个步长创建预测。 然而，当使用滞后特征时，我们仅限于预测滞后值可用的时间步。 在星期一使用lag1的特征，我们无法对星期三进行预测，因为它所需的lag值是星期二,此时还没发生。
我们将在第6课中看到处理这个问题的策略。对于这个例子，我们将只使用来自测试集的值。
只看预测值，我们可以看到我们的模型需要一个时间步来对目标序列的突然变化做出反应。 这是仅使用目标值的滞后作为特征的模型的一个常见限制。
为了改进预测，我们可以尝试找到领先指标，即可以为流感病例变化提供“预警”的时间序列。 对于我们的第二种方法，我们将在我们的训练数据中添加一些由谷歌趋势测量的与流感相关的搜索词的流行度。
将搜索词组“FluCough”与目标“FluVisits”绘制成图表表明，此类搜索词可用作领先指标：与流感相关的搜索往往在就诊前几周变得更多。
该数据集包含 129 个这样的术语，但我们只会使用其中的几个。
我们的预测有点粗略，但我们的模型似乎能够更好地预测流感访问量的突然增加，这表明搜索流行度的几个时间序列作为领先指标确实是有效的。
本课中说明的时间序列可以称为“纯周期”：它们没有明显的趋势或季节性。 尽管时间序列同时拥有趋势、季节性和周期——这三个部分的并不少见。 您可以通过为每个组件添加适当的特征来使用线性回归对此类序列进行建模。 您甚至可以将经过训练以分别学习组件的模型结合起来，我们将在下一课中学习如何使用预测混合模型进行操作。
教程 5-混合模型(Hybrid Models)
引言
线性回归擅长推断趋势，但无法学习相互作用。 XGBoost 擅长学习相互互，但无法推断趋势。 在本课中，我们将学习如何创建“混合”预测器，将互补的学习算法结合起来，让一个的优势弥补另一个的弱点。
组件(Components)和残差(Residuals)
为了设计出有效的混合模型，我们需要更好地理解时间序列是如何构建的。 到目前为止，我们已经研究了三种依赖模式：趋势、季节和周期。 许多时间序列可以通过仅由这三个组件加上一些本质上不可预测的、完全随机的残差的相加模型来描述：
series = trend + seasons + cycles + error
这个模型中的每一项我们都可以称为时间序列的组件。
模型的残差是模型训练的目标与模型做出的预测之间的差异 – 换句话说，真实曲线和拟合曲线之间的差异。 针对一个特征绘制残差，你会得到目标值的“剩余”部分，或者叫模型未能从该特征中学习到的内容。
目标序列和预测（蓝色）之间的差异给出了序列的残差。
上图左侧是Tunnel Traffic系列的一部分和第3课中的趋势-季节性曲线。减去拟合曲线后，残差在右侧。 残差包含趋势季节性模型没有学习到的Tunnel Traffic 中的所有内容。
我们可以将学习时间序列的组件想象为一个迭代过程：首先学习趋势并将其从序列中减去，然后从去趋势的残差中学习季节性并减去季节，然后学习周期并减去周期 ，最后只剩下不可预知的错误。
逐步学习 Mauna Loa CO2的组件。 从系列中减去拟合曲线（蓝色）以获得下一步中的系列。
将我们学习到的所有组件加在一起，就得到了完整的模型。 如果您在一组完整的趋势、季节和周期特征上对其进行训练建模，这基本上就是线性回归要做的事情。
添加学习的组件以获得完整的模型。
使用残差进行混合预测
在之前的课程中，我们使用单一算法（线性回归）一次学习所有组件。 但也可以对某些组件使用一种算法，而对其余组件使用另一种算法。 这样，我们总是可以为每个组件选择最佳算法。 为此，我们使用一种算法来拟合原始序列，然后使用第二种算法来拟合残差序列。
详细来说，流程是这样的：
# 1. Train and predict with first model
model_1.fit(X_train_1, y_train)
y_pred_1 = model_1.predict(X_train)
# 2. Train and predict with second model on residuals
model_2.fit(X_train_2, y_train - y_pred_1)
y_pred_2 = model_2.predict(X_train_2)
# 3. Add to get overall predictions
y_pred = y_pred_1 + y_pred_2 我们通常会根据我们希望每个模型学习的内容，来使用不同的特征集（上面的X_train_1和X_train_2）。 例如，如果我们使用第一个模型来学习趋势，那我们的第二个模型通常不需要趋势的特征。
虽然可以使用两个以上的模型，但实际上它似乎并不是特别有用。 事实上，构建混合模型的最常见策略就是我们刚刚描述的：一个简单（通常是线性）学习算法，然后是一个复杂的非线性学习器，如 GBDT 或深度神经网络，简单模型通常设计为 后面强大算法的“助手”。
设计混合模型
除了我们在本课程中概述的方式之外，还有很多方法可以组合机器学习模型。然而，成功地组合模型需要我们更深入地研究这些算法是如何运作的。
回归算法通常可以通过两种方式进行预测：要么转换特征，要么转换目标。 特征转换学习算法一些将特征作为输入的数学函数，然后将它们组合和转换以产生与训练集中的目标值匹配的输出。 线性回归和神经网络就是这种类型。
目标转换算法使用特征对训练集中的目标值进行分组，并通过对组中的值进行平均来进行预测； 一组特征只是表明要平均哪个组。 决策树和最近邻就是这种类型。
重要的是：在给定适当的特征作为输入的情况下，特征变换通常可以推断超出训练集的目标值，但目标变换的预测将始终限制在训练集的范围内。 如果时间虚拟变量(time dummy)继续计算时间步长，则线性回归继续绘制趋势线。 给定相同的时间虚拟变量(time dummy)，决策树将按照训练数据的最后一步所表明的趋势来预测未来的趋势。 决策树无法推断趋势。 随机森林和梯度提升决策树（如 XGBoost）是决策树的集合，因此它们也无法推断趋势。
决策树无法推断训练集之外的趋势。
这种差异是本课使用混合设计(hybrid design)的动机：使用线性回归来推断趋势，转换目标以消除趋势，并将 XGBoost 应用于去掉趋势的残差。 要混合神经网络（特征转换）的话，您可以将另一个模型的预测作为特征包含在内，然后神经网络会将其作为其自身预测的一部分。 对残差的拟合方法其实和梯度提升算法使用的方法是一样的，所以我们称这些为boosted hybrids； 使用预测作为特征的方法称为“堆叠”，因此我们将这些称为stacked hybrids。
在Kaggle 比赛中获胜的混合模型
为了获得灵感，这里有一些过去比赛中得分最高的解决方案：
- STL boosted with exponential smoothing - Walmart Recruiting - Store Sales Forecasting
- ARIMA and exponential smoothing boosted with GBDT - Rossmann Store Sales
- An ensemble of stacked and boosted hybrids - Web Traffic Time Series Forecasting
- Exponential smoothing stacked with LSTM neural net - M4 (non-Kaggle)
示例 - 美国零售销售额(US Retail Sales)
US Retail Sales 数据集包含美国人口普查局收集的 1992 年至 2019 年各个零售行业的月度销售数据。 我们的目标是根据早些年的销售额预测 2016-2019 年的销售额。 除了创建线性回归 + XGBoost 混合，我们还将了解如何设置用于 XGBoost 的时间序列数据集。
from pathlib import Path
from warnings import simplefilter
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
from xgboost import XGBRegressor
simplefilter("ignore")
# Set Matplotlib defaults
plt.style.use("seaborn-whitegrid")
plt.rc(
"figure",
autolayout=True,
figsize=(11, 4),
titlesize=18,
titleweight='bold',
)
plt.rc(
"axes",
labelweight="bold",
labelsize="large",
titleweight="bold",
titlesize=16,
titlepad=10,
)
plot_params = dict(
color="0.75",
style=".-",
markeredgecolor="0.25",
markerfacecolor="0.25",
)
data_dir = Path("./ts-course-data/")
industries = ["BuildingMaterials", "FoodAndBeverage"]
retail = pd.read_csv(
data_dir / "us-retail-sales.csv",
usecols=['Month'] + industries,
parse_dates=['Month'],
index_col='Month',
).to_period('D').reindex(columns=industries)
retail = pd.concat({'Sales': retail}, names=[None, 'Industries'], axis=1)
retail.head()
| Sales | ||
|---|---|---|
| Industries | BuildingMaterials | FoodAndBeverage |
| Month | ||
| 1992-01-01 | 8964 | 29589 |
| 1992-02-01 | 9023 | 28570 |
| 1992-03-01 | 10608 | 29682 |
| 1992-04-01 | 11630 | 30228 |
| 1992-05-01 | 12327 | 31677 |
首先让我们使用线性回归模型来学习每个序列的趋势。 为了演示，我们将使用二次（2 阶）趋势。 （这里的代码和上一课的代码基本相同。）虽然fit并不完美，但足以满足我们的需要。
y = retail.copy()
# Create trend features
dp = DeterministicProcess(
index=y.index, # dates from the training data
constant=True, # the intercept
order=2, # quadratic trend
drop=True, # drop terms to avoid collinearity
)
X = dp.in_sample() # features for the training data
# Test on the years 2016-2019. It will be easier for us later if we
# split the date index instead of the dataframe directly.
idx_train, idx_test = train_test_split(
y.index, test_size=12 * 4, shuffle=False,
)
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]
# Fit trend model
model = LinearRegression(fit_intercept=False)
model.fit(X_train, y_train)
# Make predictions
y_fit = pd.DataFrame(
model.predict(X_train),
index=y_train.index,
columns=y_train.columns,
)
y_pred = pd.DataFrame(
model.predict(X_test),
index=y_test.index,
columns=y_test.columns,
)
# Plot
axs = y_train.plot(color='0.25', subplots=True, sharex=True)
axs = y_test.plot(color='0.25', subplots=True, sharex=True, ax=axs)
axs = y_fit.plot(color='C0', subplots=True, sharex=True, ax=axs)
axs = y_pred.plot(color='C3', subplots=True, sharex=True, ax=axs)
for ax in axs: ax.legend([])
_ = plt.suptitle("Trends")
虽然线性回归算法能够进行多输出回归，但 XGBoost 算法却不能。为了使用 XGBoost 一次预测多个系列，我们将把这些系列从 wide 格式（每列一个时间系列）转换为 long 格式，序列按行的类别索引。
# The `stack` method converts column labels to row labels, pivoting from wide format to long
X = retail.stack() # pivot dataset wide to long
display(X.head())
y = X.pop('Sales') # grab target series| Sales | ||
|---|---|---|
| Month | Industries | |
| 1992-01-01 | BuildingMaterials | 8964 |
| FoodAndBeverage | 29589 | |
| 1992-02-01 | BuildingMaterials | 9023 |
| FoodAndBeverage | 28570 | |
| 1992-03-01 | BuildingMaterials | 10608 |
为了让 XGBoost 能够学会区分我们的两个时间序列，我们将“行业”的行标签转换为带有标签编码的分类特征。 我们还将通过从时间索引中提取月份数字来创建年度季节性特征。
# Turn row labels into categorical feature columns with a label encoding
X = X.reset_index('Industries')
# Label encoding for 'Industries' feature
for colname in X.select_dtypes(["object", "category"]):
X[colname], _ = X[colname].factorize()
# Label encoding for annual seasonality
X["Month"] = X.index.month # values are 1, 2, ..., 12
# Create splits
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]现在我们将之前所做的趋势预测转换为长格式，然后从原始系列中减去它们。 这将为我们提供 XGBoost 可以学习的去趋势（残差）系列。
# Pivot wide to long (stack) and convert DataFrame to Series (squeeze)
y_fit = y_fit.stack().squeeze() # trend from training set
y_pred = y_pred.stack().squeeze() # trend from test set
# Create residuals (the collection of detrended series) from the training set
y_resid = y_train - y_fit
# Train XGBoost on the residuals
xgb = XGBRegressor()
xgb.fit(X_train, y_resid)
# Add the predicted residuals onto the predicted trends
y_fit_boosted = xgb.predict(X_train) + y_fit
y_pred_boosted = xgb.predict(X_test) + y_pred 拟合看起来相当不错，尽管我们可以看到 XGBoost 学习的趋势与线性回归学习的趋势一样好 – 特别注意的是 XGBoost 无法弥补'BuildingMaterials' 序列中拟合不佳的趋势。
axs = y_train.unstack(['Industries']).plot(
color='0.25', figsize=(11, 5), subplots=True, sharex=True,
title=['BuildingMaterials', 'FoodAndBeverage'],
)
axs = y_test.unstack(['Industries']).plot(
color='0.25', subplots=True, sharex=True, ax=axs,
)
axs = y_fit_boosted.unstack(['Industries']).plot(
color='C0', subplots=True, sharex=True, ax=axs,
)
axs = y_pred_boosted.unstack(['Industries']).plot(
color='C3', subplots=True, sharex=True, ax=axs,
)
for ax in axs: ax.legend([])
教程 6-使用机器学习预测(forecasting-with-machine-learning)
引言
在第2课和第3课中，我们将预测视为一个简单的回归问题，我们的所有特征都来自单个输入，即时间索引。 只需生成我们想要的趋势和季节性特征，我们就可以轻松地创建未来任何时间的预测。
然而，当我们在第4课中添加滞后特征时，问题的性质发生了变化。 滞后特征要求在预测时已知滞后目标值。 一个lag_1将时间序列向前移动1步，这意味着您可以预测未来1步，但不能预测2步。
在第4课中，我们只是假设我们总是可以产生延迟到我们想要预测的时期（换句话说，每个预测都只是向前迈出一步）。 现实世界的预测通常需要更多的东西，因此在本课中，我们将学习如何针对各种情况进行预测。
定义预测任务
在设计预测模型之前，需要确定两件事：
- 预测时可获得哪些信息（特征），以及，
- 您需要预测的时间段。
**预测起点(forecast origin)**是您开始预测的时间点。 实际上，您可能会认为预测起点是您在预测时拥有的最后一条训练数据。 直到起点的一切都可以用来创建特征。
**预测范围(forecast horizon)**是您要预测的时间范围。 我们经常通过时间步长的长度来描述预测范围：例如，“1 步”预测或“5 步”预测。 预测范围描述的是目标值。
使用四个滞后特征。具有两步的前置时间(lead time)和三步预测范围 该图表示单行训练数据 -- 换句话说，用于单个预测的数据。
在起点和范围之间的时间叫预测的前置时间(lead time)(有时也叫 延迟(latency))。 预测的前置时间(lead time)指的是从起点到范围的步数：例如“向前1步”或“向前3步”预测。 在实践中，由于数据采集或处理的延迟，预测可能需要向前多个步骤再开始。
为预测准备数据
为了使用机器学习算法预测时间序列，我们需要将序列转换为可以与这些算法一起使用的dataframe。 （除非您只使用趋势和季节性等确定性特征。）
我们在第4课中看到了这个过程的前半部分，当时我们创建了一个滞后的特征集。 后半部分是准备目标值。 我们如何做到这一点取决于预测任务。
dataframe中的每一行代表一个预测。 该行的时间索引是预测范围内的第一个时间，但我们将整个范围的值安排在同一行中。 对于多步预测，这意味着我们需要一个模型来产生多个输出，每步输出一个。
import numpy as np
import pandas as pd
N = 20
ts = pd.Series(
np.arange(N),
index=pd.period_range(start='2010', freq='A', periods=N, name='Year'),
dtype=pd.Int8Dtype,
)
# Lag features
X = pd.DataFrame({
'y_lag_2': ts.shift(2),
'y_lag_3': ts.shift(3),
'y_lag_4': ts.shift(4),
'y_lag_5': ts.shift(5),
'y_lag_6': ts.shift(6),
})
# Multistep targets
y = pd.DataFrame({
'y_step_3': ts.shift(-2),
'y_step_2': ts.shift(-1),
'y_step_1': ts,
})
data = pd.concat({'Targets': y, 'Features': X}, axis=1)
data.head(10).style.set_properties(['Targets'], **{'background-color': 'LavenderBlush'}) \
.set_properties(['Features'], **{'background-color': 'Lavender'})| Targets | Features | |||||||
|---|---|---|---|---|---|---|---|---|
| y_step_3 | y_step_2 | y_step_1 | y_lag_2 | y_lag_3 | y_lag_4 | y_lag_5 | y_lag_6 | |
| Year | ||||||||
| 2010 | 2 | 1 | 0 | nan | nan | nan | nan | nan |
| 2011 | 3 | 2 | 1 | nan | nan | nan | nan | nan |
| 2012 | 4 | 3 | 2 | 0 | nan | nan | nan | nan |
| 2013 | 5 | 4 | 3 | 1 | 0 | nan | nan | nan |
| 2014 | 6 | 5 | 4 | 2 | 1 | 0 | nan | nan |
| 2015 | 7 | 6 | 5 | 3 | 2 | 1 | 0 | nan |
| 2016 | 8 | 7 | 6 | 4 | 3 | 2 | 1 | 0 |
| 2017 | 9 | 8 | 7 | 5 | 4 | 3 | 2 | 1 |
| 2018 | 10 | 9 | 8 | 6 | 5 | 4 | 3 | 2 |
| 2019 | 11 | 10 | 9 | 7 | 6 | 5 | 4 | 3 |
上面说明了如何通过类似于定义预测图来准备数据集：使用五个滞后特征的一个具有2步长度的前置时间的3步预测长度的预测任务。起始时间序列是y_step_1。 我们可以填写或删除缺失值。
多步预测策略
有许多策略可以生成预测所需的多个目标步骤。 我们将概述四种常见的策略，每种策略都有优点和缺点。
多输出模型
使用能够产生多个输出的模型。 线性回归和神经网络都可以产生多个输出。 这种策略简单而有效，但不可能适用于您可能想要使用的每种算法。 例如，XGBoost 无法做到这一点。
直接策略
为预测范围的每一步训练一个单独的模型：一个模型预测向前1步，另一个预测向前2步，依此类推。 向前1步预测与向前2步（等等）是不同的问题，因此它可以让不同的模型对每一步进行预测。 缺点是训练大量模型的计算成本可能很高。
递归策略
训练一个单步模型并使用其预测值来更新下一步的滞后特征。 使用递归方法，我们将模型的1步预测反馈回同一模型，以用作下一个预测步骤的滞后特征。 这样我们只需要训练一个模型，但由于错误会一步一步地传播，因此长期预测可能不准确。（长期预测 误差）
DirRec策略（直接递归）
直接策略和递归策略的组合：为每个步骤训练一个模型，并使用来自先前步骤的预测作为新的滞后特征。 逐步地，每个模型都会获得一个额外的滞后输入。 由于每个模型总是有一组最新的滞后特征，DirRec 策略可以比直接策略更好地捕获序列依赖，但它也可能遭受像递归策略这样的错误传播。
示例 - 流感趋势
在此示例中，我们将对第4课中的流感趋势数据应用多输出和直接策略，对训练期之后的多个星期进行真实预测。
我们将我们的预测任务定义为8周的时间范围和1周的前置时间。 换句话说，我们将从下周开始预测接下来八周的流感病例数。
下面设置示例并定义辅助函数 plot_multistep。
from pathlib import Path
from warnings import simplefilter
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from xgboost import XGBRegressor
simplefilter("ignore")
# Set Matplotlib defaults
plt.style.use("seaborn-whitegrid")
plt.rc("figure", autolayout=True, figsize=(11, 4))
plt.rc(
"axes",
labelweight="bold",
labelsize="large",
titleweight="bold",
titlesize=16,
titlepad=10,
)
plot_params = dict(
color="0.75",
style=".-",
markeredgecolor="0.25",
markerfacecolor="0.25",
)
%config InlineBackend.figure_format = 'retina'
def plot_multistep(y, every=1, ax=None, palette_kwargs=None):
palette_kwargs_ = dict(palette='husl', n_colors=16, desat=None)
if palette_kwargs is not None:
palette_kwargs_.update(palette_kwargs)
palette = sns.color_palette(**palette_kwargs_)
if ax is None:
fig, ax = plt.subplots()
ax.set_prop_cycle(plt.cycler('color', palette))
for date, preds in y[::every].iterrows():
preds.index = pd.period_range(start=date, periods=len(preds))
preds.plot(ax=ax)
return ax
data_dir = Path("../input/ts-course-data")
flu_trends = pd.read_csv(data_dir / "flu-trends.csv")
flu_trends.set_index(
pd.PeriodIndex(flu_trends.Week, freq="W"),
inplace=True,
)
flu_trends.drop("Week", axis=1, inplace=True)首先，我们将准备我们的目标系列（每周流感的访问）以进行多步预测。 一旦完成，训练和预测将非常简单。
def make_lags(ts, lags, lead_time=1):
return pd.concat(
{
f'y_lag_{i}': ts.shift(i)
for i in range(lead_time, lags + lead_time)
},
axis=1)
# Four weeks of lag features
y = flu_trends.FluVisits.copy()
X = make_lags(y, lags=4).fillna(0.0)
def make_multistep_target(ts, steps):
return pd.concat(
{f'y_step_{i + 1}': ts.shift(-i)
for i in range(steps)},
axis=1)
# Eight-week forecast
y = make_multistep_target(y, steps=8).dropna()
# Shifting has created indexes that don't match. Only keep times for
# which we have both targets and features.
y, X = y.align(X, join='inner', axis=0)
多输出模型
我们将使用线性回归作为多输出策略。 一旦我们为多个输出准备好数据，训练和预测就和往常一样了。
# Create splits
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle=False)
model = LinearRegression()
model.fit(X_train, y_train)
y_fit = pd.DataFrame(model.predict(X_train), index=X_train.index, columns=y.columns)
y_pred = pd.DataFrame(model.predict(X_test), index=X_test.index, columns=y.columns)
请记住，多步模型将为用作输入的每个实例生成完整的预测。 训练集中有 269 周，测试集中有 90 周，我们现在对这些周中的每一周都有一个 8 周的预测。
train_rmse = mean_squared_error(y_train, y_fit, squared=False)
test_rmse = mean_squared_error(y_test, y_pred, squared=False)
print((f"Train RMSE: {train_rmse:.2f}\n" f"Test RMSE: {test_rmse:.2f}"))
palette = dict(palette='husl', n_colors=64)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(11, 6))
ax1 = flu_trends.FluVisits[y_fit.index].plot(**plot_params, ax=ax1)
ax1 = plot_multistep(y_fit, ax=ax1, palette_kwargs=palette)
_ = ax1.legend(['FluVisits (train)', 'Forecast'])
ax2 = flu_trends.FluVisits[y_pred.index].plot(**plot_params, ax=ax2)
ax2 = plot_multistep(y_pred, ax=ax2, palette_kwargs=palette)
_ = ax2.legend(['FluVisits (test)', 'Forecast'])Train RMSE: 389.12 Test RMSE: 582.33
直接策略
XGBoost 不能为回归任务生成多个输出。 但是通过简单的方式，我们仍然可以使用它来生成多步预测。 这就是使用 scikit-learn 的 MultiOutputRegressor 来包装它。
from sklearn.multioutput import MultiOutputRegressor
model = MultiOutputRegressor(XGBRegressor())
model.fit(X_train, y_train)
y_fit = pd.DataFrame(model.predict(X_train), index=X_train.index, columns=y.columns)
y_pred = pd.DataFrame(model.predict(X_test), index=X_test.index, columns=y.columns)
这里的 XGBoost 显然在训练集上过度拟合。 但在测试集上，它似乎能够比线性回归模型更好地捕捉到流感季节的一些动态。 通过一些超参数调整它可能会做得更好。
train_rmse = mean_squared_error(y_train, y_fit, squared=False)
test_rmse = mean_squared_error(y_test, y_pred, squared=False)
print((f"Train RMSE: {train_rmse:.2f}\n" f"Test RMSE: {test_rmse:.2f}"))
palette = dict(palette='husl', n_colors=64)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(11, 6))
ax1 = flu_trends.FluVisits[y_fit.index].plot(**plot_params, ax=ax1)
ax1 = plot_multistep(y_fit, ax=ax1, palette_kwargs=palette)
_ = ax1.legend(['FluVisits (train)', 'Forecast'])
ax2 = flu_trends.FluVisits[y_pred.index].plot(**plot_params, ax=ax2)
ax2 = plot_multistep(y_pred, ax=ax2, palette_kwargs=palette)
_ = ax2.legend(['FluVisits (test)', 'Forecast'])
要使用 DirRec 策略，您只需将 MultiOutputRegressor 替换为另一个 scikit-learn 包装器 RegressorChain。 递归策略就需要我们需要自己编写代码来实现。
轮到你了
为商店销售额创建预测数据集 并应用 DirRec 策略。
2002
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
