利用信号处理技术提高时间序列预测性能
时间序列数据建模和预测是复杂的课题。有很多技术可以用来提高预测任务的模型性能。在这里，我们将讨论一种技术，它可以改进预测模型对时间特征的吸收和利用，从而使模型更具普适性。主要关注的是在训练中为时间序列预测模型创建季节性特征的过程 – 如果在特征创建过程中包含傅里叶变换，可以轻松提高预测模型的性能。
这篇文章假设你熟悉时间序列预测的基本知识 – 我们不会讨论这个主题的总体，而只会讨论其某个方面的改进。如果想了解时间序列预测的介绍，请参考Kaggle课程 – 这里讨论的技术是基于他们关于季节性的课程的。
数据集
考虑澳大利亚季度水泥生产数据集，显示了澳大利亚的波特兰水泥的总季度生产量（以百万吨为单位）从1956年第一季度到2014年第一季度。
df = pd.read_csv('Quarterly_Australian_Portland_Cement_production_776_10.csv', usecols=['time', 'value'])# 将年份浮点数转换为适当的日期时间格式df['time'] = df['time'].apply(lambda x: str(int(x)) + '-' + str(int(1 + 12 * (x % 1))).rjust(2, '0'))df['time'] = pd.to_datetime(df['time'])df = df.set_index('time').to_period()df.rename(columns={'value': 'production'}, inplace=True)
productiontime 1956Q1 0.4651956Q2 0.5321956Q3 0.5611956Q4 0.5701957Q1 0.529... ...2013Q1 2.0492013Q2 2.5282013Q3 2.6372013Q4 2.5652014Q1 2.229[233 rows x 1 columns]
这是带有趋势、季节性成分和其他属性的时间序列数据。观测结果是按季度进行的，跨越了几十年。让我们首先查看季节性成分，使用周期图（所有代码都包含在链接的伴随笔记本中）。
周期图显示了频谱分量（季节性分量）的功率密度。最强的分量是周期为4个季度或1年的分量。这确认了视觉印象，即图中最强的上升和下降变化大约每十年发生10次。还有一个周期为2个季度的分量 – 这是相同的季节性现象，只是意味着4个季度的周期性不是简单的正弦波，而是具有更复杂的形状。为简单起见，我们将忽略周期大于4的分量。
如果您尝试对这些数据拟合模型，也许是为了预测超出数据集末尾的时间的水泥生产情况，那么捕捉这种年度季节性并在训练中为模型提供这些特征或一组特征是一个好主意。让我们看一个例子。
季节性特征
季节性分量可以以多种不同的方式建模。通常，它们被表示为时间虚拟变量或正弦余弦对。这些是具有与要建模的季节性相等的周期或其子倍数的合成特征。
年度时间虚拟变量：
seasonal_year = DeterministicProcess(index=df.index, constant=False, seasonal=True).in_sample()print(seasonal_year)
s(1,4) s(2,4) s(3,4) s(4,4)time 1956Q1 1.0 0.0 0.0 0.01956Q2 0.0 1.0 0.0 0.01956Q3 0.0 0.0 1.0 0.01956Q4 0.0 0.0 0.0 1.01957Q1 1.0 0.0 0.0 0.0... ... ... ... ...2013Q1 1.0 0.0 0.0 0.02013Q2 0.0 1.0 0.0 0.02013Q3 0.0 0.0 1.0 0.02013Q4 0.0 0.0 0.0 1.02014Q1 1.0 0.0 0.0 0.0[233 rows x 4 columns]每年的正弦-余弦对：
cfr = CalendarFourier(freq='Y', order=2)seasonal_year_trig = DeterministicProcess(index=df.index, seasonal=False, additional_terms=[cfr]).in_sample()with pd.option_context('display.max_columns', None, 'display.expand_frame_repr', False): print(seasonal_year_trig)
sin(1,freq=A-DEC) cos(1,freq=A-DEC) sin(2,freq=A-DEC) cos(2,freq=A-DEC)time 1956Q1 0.000000 1.000000 0.000000 1.0000001956Q2 0.999963 0.008583 0.017166 -0.9998531956Q3 0.017166 -0.999853 -0.034328 0.9994111956Q4 -0.999963 -0.008583 0.017166 -0.9998531957Q1 0.000000 1.000000 0.000000 1.000000... ... ... ... ...2013Q1 0.000000 1.000000 0.000000 1.0000002013Q2 0.999769 0.021516 0.043022 -0.9990742013Q3 0.025818 -0.999667 -0.051620 0.9986672013Q4 -0.999917 -0.012910 0.025818 -0.9996672014Q1 0.000000 1.000000 0.000000 1.000000[233 rows x 4 columns]时间虚拟变量可以捕捉季节性现象的非常复杂的波形。另一方面，如果周期是N，则你需要N个时间虚拟变量，所以如果周期很长，则需要大量的虚拟列，这可能是不可取的。对于非惩罚模型，通常只使用N-1个虚拟变量 —— 为了避免多重共线性问题而将其中一个虚拟变量移除（我们在此忽略了这个问题）。
正弦/余弦对可以模拟任意长的时间周期。每对将模拟一个纯正弦波形与一些初始相位。随着季节特征波形变得更加复杂，你将需要添加更多的对（增加CalendarFourier的输出阶数）。
（副注：对于我们想要建模的每个周期，我们使用一个正弦/余弦对。实际上，我们只希望得到一列A*sin(2*pi*t/T + phi)其中A是模型分配给该列的权重，t是时间序列的时间索引，T是组成周期。但在拟合数据过程中，模型无法调整初始相位phi —— 正弦值是预先计算的。幸运的是，上面的公式等同于：A1*sin(2*pi*t/T) + A2*cos(2*pi*t/T)，模型只需要找到权重A1和A2）。
TLDR：这两种技术的共同之处在于它们通过一组重复特征来表示季节性，其值的周期性与正在建模的时间周期一样频繁（时间虚拟变量和基本正弦/余弦对），或在每个周期内循环多次（高阶正弦/余弦）。其中每一个特征的值在一个固定区间内变化：从0到1，或从-1到1。这里所有要建模季节性的特征都在这些特征集合中。
让我们看看在这些季节性特征上拟合线性模型会发生什么。但首先，我们需要向用于训练模型的特征集合中添加趋势特征。
拟合线性模型
让我们创建趋势特征，然后将其与上面生成的seasonal_year时间虚拟变量组合在一起：
trend_order = 2trend_year = DeterministicProcess(index=df.index, constant=True, order=trend_order).in_sample()X = trend_year.copy()X = X.join(seasonal_year)
const trend trend_squared s(1,4) s(2,4) s(3,4) s(4,4)time 1956Q1 1.0 1.0 1.0 1.0 0.0 0.0 0.01956Q2 1.0 2.0 4.0 0.0 1.0 0.0 0.01956Q3 1.0 3.0 9.0 0.0 0.0 1.0 0.01956Q4 1.0 4.0 16.0 0.0 0.0 0.0 1.01957Q1 1.0 5.0 25.0 1.0 0.0 0.0 0.0... ... ... ... ... ... ... ...2013Q1 1.0 229.0 52441.0 1.0 0.0 0.0 0.02013Q2 1.0 230.0 52900.0 0.0 1.0 0.0 0.02013Q3 1.0 231.0 53361.0 0.0 0.0 1.0 0.02013Q4 1.0 232.0 53824.0 0.0 0.0 0.0 1.02014Q1 1.0 233.0 54289.0 1.0 0.0 0.0 0.0[233 rows x 7 columns]这是将用于训练/验证模型的 X 数据帧（特征）。我们使用二次趋势特征对数据进行建模，再加上年度季节性所需的 4 个时间虚变量。y 数据帧（目标）将仅包含水泥产量数字。
让我们从数据中划分出一个包含 2010 年观察值的验证集。我们将在上面显示的 X 特征和代表水泥产量的 y 目标（测试部分）上拟合一个线性模型，然后在验证中评估模型性能。我们还将使用仅趋势的模型执行上述所有操作，该模型仅适配趋势特征，而忽略季节性。
def do_forecast(X, index_train, index_test, trend_order): X_train = X.loc[index_train] X_test = X.loc[index_test] y_train = df['production'].loc[index_train] y_test = df['production'].loc[index_test] model = LinearRegression(fit_intercept=False) _ = model.fit(X_train, y_train) y_fore = pd.Series(model.predict(X_test), index=index_test) y_past = pd.Series(model.predict(X_train), index=index_train) trend_columns = X_train.columns.to_list()[0 : trend_order + 1] model_trend = LinearRegression(fit_intercept=False) _ = model_trend.fit(X_train[trend_columns], y_train) y_trend_fore = pd.Series(model_trend.predict(X_test[trend_columns]), index=index_test) y_trend_past = pd.Series(model_trend.predict(X_train[trend_columns]), index=index_train) RMSLE = mean_squared_log_error(y_test, y_fore, squared=False) print(f'RMSLE: {RMSLE}') ax = df.plot(**plot_params, title='AUS 水泥产量 - 预测') ax = y_past.plot(color='C0', label='向后拟合') ax = y_fore.plot(color='C3', label='预测') ax = y_trend_past.plot(ax=ax, color='C0', linewidth=3, alpha=0.333, label='趋势过去') ax = y_trend_fore.plot(ax=ax, color='C3', linewidth=3, alpha=0.333, label='趋势未来') _ = ax.legend()do_forecast(X, index_train, index_test, trend_order)
RMSLE: 0.03846449744356434
蓝色是训练，红色是验证。完整模型是清晰而细线，仅趋势模型是宽而模糊的线。
这不算太糟糕，但有一个明显的问题：模型已经学会了一个恒定振幅的年度季节变化。尽管年度波动实际上是递增的，但模型只能坚持固定振幅变化。显然，这是因为我们只给了模型固定振幅的季节性特征，而特征中没有其他内容（目标滞后等）来帮助克服这个问题。
让我们深入挖掘信号（数据），看看是否有任何内容可以帮助解决固定振幅问题。
傅里叶光谱图
周期图将突出显示信号中的所有谱成分（数据中的所有季节性成分），并提供其整体强度的概览，但它是一个总数，是在整个时间间隔上进行的。它无法说明每个季节性成分的幅度在整个数据集中是否随时间变化。
为了捕捉这种变化，您必须使用傅里叶光谱图。它类似于周期图，但是在整个数据集上重复执行，涵盖许多时间窗口。周期图和光谱图都可以在 scipy 库中使用。
让我们绘制具有 2 和 4 个季度周期的季节性成分的光谱图，如上所述。和往常一样，完整的代码在结尾链接的伴随笔记本中。
compute_spectrum(df['production'], 4, 0.1)plot_spectrogram(spectrum, figsize_x=10)
这个图表显示的是时间上2个季度和4个季度的成分的振幅，即振幅的“强度”。它们最初相当弱，但在2010年左右变得非常强，与本文顶部的数据集绘图中所看到的变化相匹配。
如果我们不是将固定振幅的季节特征提供给模型，而是根据光谱图所告诉我们的来调整这些特征的振幅，那么最终的模型在验证过程中会表现得更好吗？
调整季节成分
让我们选择4个季度的季节成分。我们可以对这个成分的order=2趋势在训练数据上拟合一个线性模型（称为包络模型），并将该趋势延伸到验证/测试数据（model.predict()）：
envelope_features = DeterministicProcess(index=X.index, constant=True, order=2).in_sample()spec4_train = compute_spectrum(df['production'].loc[index_train], max_period=4)spec4_trainspec4_model = LinearRegression()spec4_model.fit(envelope_features.loc[spec4_train.index], spec4_train['4.0'])spec4_regress = pd.Series(spec4_model.predict(envelope_features), index=X.index)ax = spec4_train['4.0'].plot(label='分量包络', color='灰色')spec4_regress.loc[spec4_train.index].plot(ax=ax, color='C0', label='包络回归：过去')spec4_regress.loc[index_test].plot(ax=ax, color='C3', label='包络回归：未来')_ = ax.legend()
蓝线是训练数据中4个季度成分变化的强度，拟合为二次趋势（order=2）。红线是相同的东西，在验证数据上延伸（预测）。
我们已经建模了4个季度季节成分的时间变化。我们可以将包络模型的输出乘以对应于4个季度季节成分的时间虚拟变量：
spec4_regress = spec4_regress / spec4_regress.mean()season_columns = ['s(1,4)', 's(2,4)', 's(3,4)', 's(4,4)']for c in season_columns: X[c] = X[c] * spec4_regressprint(X)
const trend trend_squared s(1,4) s(2,4) s(3,4) s(4,4)time 1956Q1 1.0 1.0 1.0 0.179989 0.000000 0.000000 0.0000001956Q2 1.0 2.0 4.0 0.000000 0.181109 0.000000 0.0000001956Q3 1.0 3.0 9.0 0.000000 0.000000 0.182306 0.0000001956Q4 1.0 4.0 16.0 0.000000 0.000000 0.000000 0.1835811957Q1 1.0 5.0 25.0 0.184932 0.000000 0.000000 0.000000... ... ... ... ... ... ... ...2013Q1 1.0 229.0 52441.0 2.434701 0.000000 0.000000 0.0000002013Q2 1.0 230.0 52900.0 0.000000 2.453436 0.000000 0.0000002013Q3 1.0 231.0 53361.0 0.000000 0.000000 2.472249 0.0000002013Q4 1.0 232.0 53824.0 0.000000 0.000000 0.000000 2.4911392014Q1 1.0 233.0 54289.0 2.510106 0.000000 0.000000 0.000000[233 rows x 7 columns]这4个时间哑变量现在不再是0或1了。它们已经乘以了组件包络，现在它们的振幅随时间变化，就像包络一样。
重新训练模型
让我们再次训练主模型，现在使用修改后的时间哑变量。
do_forecast(X, index_train, index_test, trend_order)
RMSLE: 0.02546321729737165
我们不追求完美的匹配，因为这只是一个玩具模型，但很明显模型表现更好，更贴近目标变量的四季度变化。性能指标也提高了51%，这一点也不错。
最后的评论
提高模型性能是一个广泛的课题。任何模型的行为不仅取决于单个特征或单一技术。然而，如果你想尽可能地提取给定模型的全部潜力，你可能应该提供有意义的特征给它。时间哑变量或正弦/余弦傅里叶对在反映调节它们的季节性在时间上变化时更有意义。
通过傅里叶变换来调整季节性分量的包络对于线性模型尤其有效。增强树的受益并不多，但无论使用何种模型，你仍然可以看到改进。如果你正在使用混合模型，其中一个线性模型处理确定性特征（日历等），而增强树模型处理更复杂的特征，那么线性模型“做对了”就尤为重要，这样其他模型就会减少工作量。
你使用来调整季节特征的包络模型只能在训练数据上进行训练，并且在训练过程中不能看到任何测试数据，就像任何其他模型一样。一旦你调整了包络，时间哑变量中就包含了来自目标变量的信息——它们不再纯粹是确定性特征，不能提前在任意预测时间范围上计算。因此，在这一点上，如果不小心，信息可能从验证/测试数据中泄漏回训练数据。
链接
本文章使用的数据集在此处以公共领域（CC0）许可证提供：
KEY2STATS© 2023 KEY2STATS – RStudio®是RStudio, Inc.的注册商标。AP®是The College的注册商标…
www.key2stats.com
本文章中使用的代码可以在这里找到：
misc/seasonal_features_fourier_article/cement.ipynb at master · FlorinAndrei/misc
random stuff. Contribute to FlorinAndrei/misc development by creating an account on GitHub.
github.com
一个提交到Kaggle Store Sales – Time Series Forecasting竞赛的笔记本，使用了本文中描述的想法：
Fourier spectrogram, ensemble models
Explore and run machine learning code with Kaggle Notebooks | Using data from Store Sales – Time Series Forecasting
www.kaggle.com
一个包含向Kaggle提交的笔记本的开发版本的GitHub存储库在这里：
GitHub – FlorinAndrei/timeseries-forecasting-fourier-time-dummies: Time series forecasting with…
时间序列预测与傅立叶调整时间虚拟变量 – GitHub…
github.com
本文中所使用的所有图片和代码均由作者创建。