有什么方法可以在 python 中使用 scikit-learn 预测月度时间序列？
any way to predict monthly time series with scikit-learn in python?
我想通过使用每月时间序列中的多个特征来预测 product' sales_index。一开始，我开始使用 ARMA、ARIMA 来执行此操作，但输出对我来说不是很令人满意。在我的尝试中，我只是使用 dates 和 sales 列进行预测，输出对我来说并不现实。我认为我应该包含更多功能列来预测 sales_index 列。但是，我想知道是否有任何方法可以通过使用每月时间序列中的多个特征来进行此预测。我没有使用 scikit-learn 完成很多时间序列。谁能指出我这样做的任何可能方法？有什么想法吗？
我尝试使用 ARMA/ARIMA:
这是reproducible monthly time series data on this gist，这是我目前的尝试：
import pandas as pd
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import matplotlib.pyplot as plt
df = pd.read_csv("tsdf.csv", sep=",")
dates = pd.date_range(start='2015-01', freq='MS', periods=len(df))
df.set_index(dates,inplace=True)
train = df[df.index < '2019-01']
test = df[df.index >= '2019-01']
model = ARMA(train['sales_index'],order=(2,0))
model_fit = model.fit()
predictions = model_fit.predict(start=len(train), end=len(train)+len(test)-1, dynamic=False)
# plot results
plt.figure(figsize=(12,6))
plt.plot(test['sales_index'])
plt.plot(predictions, color='red')
plt.show()
这是我当前尝试的输出：
在我的尝试中，我只是简单地使用 df['sales_index] 和 df['dates'] 作为 ARMA 模型。显然这样做，预测输出不是很真实和信息丰富。我在想是否有任何方法可以将 df['sales_index'] 以外的所有特征列提供给模型以预测 df['sales_index']。我想不出更好的方法来使用 ARMA 模型。
也许scikit-learn 可以为这个预测起到更好的作用。我不确定如何使用 sklearn 来完成此时间序列分析。谁能指出这个时间序列可能的 sklearn 解决方案？在 sklearn 中是否有可能这样做？任何可能的想法？谢谢
您可以在初始化模型和进行预测时使用可选的 exog 参数向 ARMA 模型添加额外的特征。
例如，要添加一些您的功能：
# initialize the model
model = ARMA(train['sales_index'],
exog=train[['slg_adt_ctl', 'slg_bbs', 'retail_sales_index']],
order=(2,0))
model_fit = model.fit() # fit the model
# make predictions
predictions = model_fit.predict(start=len(train),
end=len(train)+len(test)-1,
exog=test[['slg_adt_ctl', 'slg_bbs', 'retail_sales_index']],
dynamic=False)
当我们生成预测图时，我们现在获得了一些额外的预测能力。
概览
通过使用Scikit-Learn 库，可以考虑使用不同的决策树来预测数据。在此示例中，我们将使用 AdaBoostRegressor，但也可以切换到 RandomForestRegressor 或任何其他可用的树。因此，通过选择树，我们应该意识到去除数据的趋势，通过这种方式，我们举例说明通过分别对数据进行差分和应用对数变换来控制时间序列的均值和方差。
正在准备数据
时间序列有两个基本组成部分，均值和方差。理想情况下，我们希望控制这些组件，对于可变性，我们可以简单地对数据应用对数变换，对于趋势我们可以区分它，我们稍后会看到。
此外，对于这种方法，我们认为实际值 y_t 可以用两个先验值 y_t-1 和 y_t-2 来解释。您可以通过修改 range 函数的输入来使用这些滞后值。
# Load data
tsdf = pd.read_csv('tsdf.csv', sep="\t")
# For simplicity I just take the target variable and the date
tsdf = tsdf[['sales_index', 'dates']]
# Log value
tsdf['log_sales_index'] = np.log(tsdf.sales_index)
# Add previous n-values
for i in range(3):
tsdf[f'sales_index_lag_{i+1}'] = tsdf.sales_index.shift(i+1)
# For simplicity we drop the null values
tsdf.dropna(inplace=True)
tsdf[f'log_sales_index_lag_{i+1}'] = np.log(tsdf[f'sales_index_lag_{i+1}'])
tsdf[f'log_difference_{i+1}'] = tsdf.log_sales_index - tsdf[f'log_sales_index_lag_{i+1}']
数据准备就绪后，我们会得到类似于下图的结果。
数据是否平稳？
要控制时间序列的平均分量，我们应该对数据进行一些差分。为了确定是否需要这一步，我们可以进行单位根检验。有几个对此做出不同假设的测试，可以找到一些单位根测试的列表 here。
为简单起见，我们将考虑 KPSS 测试，我们假设数据是平稳的零假设，特别是，它假设围绕均值或线性趋势保持平稳。
from arch.unitroot import KPSS
# Test for stationary
kpss_test = KPSS(tsdf.sales_index)
# Test summary
print(kpss_test.summary().as_text())
我们看到 P-value = .280 比通常约定的 0.05 大。因此，我们需要对数据应用一阶差分。作为旁注，可以迭代地执行此测试以了解应将多少差异应用于数据。
在下图中，我们看到了原始数据与对数一阶差分的比较，注意时间序列的最后这些值突然发生变化，这似乎是结构性变化，但我们不会深入研究它。如果您想深入了解此主题，布鲁斯·汉森 (Bruce Hansen) 的这些 slides 很有用。
plt.figure(figsize=(12, 6))
plt.subplot(1,2,1)
plt.plot(tsdf.sales_index)
plt.title('Original Time Series')
plt.subplot(1,2,2)
plt.plot(tsdf.log_difference_1)
plt.title('Log first difference Time Series')
决策树模型
正如我们之前所说，我们正在考虑决策树模型，在使用它们时应该注意从时间序列中删除趋势。例如，如果你有上升趋势，tress 不擅长预测下降趋势。在下面的代码示例中，我选择了 AdaBoostRegressor，但您可以自由选择其他树模型。此外，请注意 log_difference_1 被认为是由 log_difference_2 和 log_difference_3 解释的。
Note. Your dataset has other covariates as aus_avg_rain or slg_adt_ctl, so to consider them for predicting you could apply as well lag values on them.
from sklearn.model_selection import train_test_split
from sklearn.ensemble import AdaBoostRegressor
# Forecast difference of log values
X, Y = tsdf[['log_difference_2', 'log_difference_3']], tsdf['log_difference_1']
# Split in train-test
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, shuffle=False, random_state=0)
# Initialize the estimator
mdl_adaboost = AdaBoostRegressor(n_estimators=500, learning_rate=0.05)
# Fit the data
mdl_adaboost.fit(X_train, Y_train)
# Make predictions
pred = mdl_adaboost.predict(X_test)
test_size = X_test.shape[0]
评估预测
test_size = X_test.shape[0]
plt.plot(list(range(test_size)), np.exp(tsdf.tail(test_size).log_sales_index_lag_1 + pred), label='predicted', color='red')
plt.plot(list(range(test_size)), tsdf.tail(test_size).sales_index, label='real', color='blue')
plt.legend(loc='best')
plt.title('Predicted vs Real with log difference values')
决策树模型似乎准确地预测了真实值。然而，为了评估模型性能，我们应该考虑一个评估指标，可以在这个 article 上找到关于这个主题的很好的介绍，随意选择一个对你的方法更方便的。我将使用 scikit-learn 中的 TimeSeriesSplit 函数，通过平均绝对误差来评估模型的误差。
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error
X, Y = np.array(tsdf[['log_difference_2', 'log_difference_3']]), np.array(tsdf['log_difference_1'])
# Initialize a TimeSeriesSplitter
tscv = TimeSeriesSplit(n_splits=5)
# Retrieve log_sales_index and sales_index to unstransform data
tsdf_log_sales_index = np.array(tsdf.copy().reset_index().log_sales_index_lag_1)
tsdf_sales_index = np.array(tsdf.copy().reset_index().sales_index_lag_1)
# Dict to store metric value at every iteration
metric_iter = {}
for idx, val in enumerate(tscv.split(X)):
train_i, test_i = val
X_train, X_test = X[train_i], X[test_i]
Y_train, Y_test = Y[train_i], Y[test_i]
# Initialize the estimator
mdl_adaboost = AdaBoostRegressor(n_estimators=500, learning_rate=0.05)
# Fit the data
mdl_adaboost.fit(X_train, Y_train)
# Make predictions
pred = mdl_adaboost.predict(X_test)
# Unstransform predictions
pred_untransform = [np.exp(val_test + val_pred) for val_test, val_pred in zip(tsdf_log_sales_index[test_i], pred)]
# Real value
real = tsdf_sales_index[test_i]
# Store metric
metric_iter[f'iter_{idx + 1}'] = mean_absolute_error(real, pred_untransform)
现在我们看到平均 MAE 误差非常低。
print(f'Average MAE error: {np.mean(list(metric_iter.values()))}')
>>> Average MAE error: 17.631090959806535
这有两部分：
- 一个有效的预测代码解决方案。
- 几个话题的简单讨论
解决方案
我将建议一种略有不同且更抽象的方法：使用 Darts which is built on top of scikit-learn. It includes 'A list' libraries (e.g., pandas and NumPy) you'd expect but also a few that you would need to have fairly deep knowledge of the field to think to include (e.g., Torch, Prophet, and Holidays). Furthermore, it has a number of models already included. You can see a bit more for both here。重要说明：特定库是 u8darts，而不是飞镖。它们很相似——并且具有相似的依赖关系——但它们并不相同。
使用它，你可以很容易地开发出一个具有很好结果的预测模型。例如，整个代码，包括合并重命名的列标题（以及删除未重命名的列和重复的行）只有几行。
import pandas as pd
import matplotlib.pyplot as plt
from darts import TimeSeries
from darts.models import ExponentialSmoothing
df = pd.read_csv("../data/tsdf.csv", sep="\t")
df.drop(['Unnamed: 0', 'aus_slg_fmCatl'], axis=1, inplace = True)
df.columns = ['Date', 'adult_cattle(head)','Bulls Bullocks & Steers(head)', 'Cows&Heifers(head)',
'Calves(head)', 'Beef(tonnes)', 'Veal(tonnes)','Adult cattle(kg/head)', 'Calves(kg/head)',
'aust-avg_rain','US/AUS_ExchangeRate', 'sales_index', 'retail_sales_index']
df.drop_duplicates(subset=['Date'], inplace=True)
series = TimeSeries.from_dataframe(df, 'Date', 'Beef(tonnes)').values
train, val = series.split_before(pd.Timestamp('2019-01-01'))
model = ExponentialSmoothing()
model.fit(train)
prediction = model.predict(len(val))
series.plot(label='actual')
prediction.plot(label='forecast', lw=2)
plt.legend()
此外，他们的回购有一个 notebook with additional models and examples
结果
讨论
针对我们收到的评论中的讨论，我认为有几点需要说明。
去除季节性
如果您正在为预测模型拍摄，这是您很少想做的事情。如果您知道每个冬天 home/office 供暖燃料价格都会上涨，但决定将其排除在外 - 您会被很乐意拿走您的钱的商人活活吃掉。
虚拟数据
这里的诀窍是您不知道您的虚拟数据是否与真实数据看起来足够相似。然而，解决方案是一个不同的问题。
“在 Python 中使用 R 库...”
看看r2py
ARMA 和 ARIMA
... 不太可能产生丰硕的结果，部分原因是误差项假定误差为 IID，但市场存在偏差。指数平滑效果很好，因为它替代了
“但是当我 运行 它反对 sales_index...”
sales_index 将由几个自变量组成。正如您所期望的那样，肉类生产（奇怪的是与降雨量呈负相关）和汇率。但是，您的模型中不包括其他国家/地区的生产数据 (100% pure organic, grass feed, 'AAA' Alberta beef)、税收、国内牛肉生产、运输等。这就是为什么您可能不想与拥有这两者的专业交易员对抗宽客和认真的领域专业知识。最后，我要指出的是，在时间上没有任何迹象表明这是正常的还是不正常的。
从 2015 年到 2017 年似乎有 1.5 年的上升期，随后一直下降到 2019 年。似乎 正在重复，但 2020 年的突然变化导致人们相信这是一种失常。但数据太小，很难判断任何有效性。在这种情况下，趋势线或通道会更好。
我想通过使用每月时间序列中的多个特征来预测 product' sales_index。一开始，我开始使用 ARMA、ARIMA 来执行此操作，但输出对我来说不是很令人满意。在我的尝试中，我只是使用 dates 和 sales 列进行预测，输出对我来说并不现实。我认为我应该包含更多功能列来预测 sales_index 列。但是，我想知道是否有任何方法可以通过使用每月时间序列中的多个特征来进行此预测。我没有使用 scikit-learn 完成很多时间序列。谁能指出我这样做的任何可能方法？有什么想法吗？
我尝试使用 ARMA/ARIMA:
这是reproducible monthly time series data on this gist，这是我目前的尝试：
import pandas as pd
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import matplotlib.pyplot as plt
df = pd.read_csv("tsdf.csv", sep=",")
dates = pd.date_range(start='2015-01', freq='MS', periods=len(df))
df.set_index(dates,inplace=True)
train = df[df.index < '2019-01']
test = df[df.index >= '2019-01']
model = ARMA(train['sales_index'],order=(2,0))
model_fit = model.fit()
predictions = model_fit.predict(start=len(train), end=len(train)+len(test)-1, dynamic=False)
# plot results
plt.figure(figsize=(12,6))
plt.plot(test['sales_index'])
plt.plot(predictions, color='red')
plt.show()
这是我当前尝试的输出：
在我的尝试中，我只是简单地使用 df['sales_index] 和 df['dates'] 作为 ARMA 模型。显然这样做，预测输出不是很真实和信息丰富。我在想是否有任何方法可以将 df['sales_index'] 以外的所有特征列提供给模型以预测 df['sales_index']。我想不出更好的方法来使用 ARMA 模型。
也许scikit-learn 可以为这个预测起到更好的作用。我不确定如何使用 sklearn 来完成此时间序列分析。谁能指出这个时间序列可能的 sklearn 解决方案？在 sklearn 中是否有可能这样做？任何可能的想法？谢谢
您可以在初始化模型和进行预测时使用可选的 exog 参数向 ARMA 模型添加额外的特征。
例如，要添加一些您的功能：
# initialize the model
model = ARMA(train['sales_index'],
exog=train[['slg_adt_ctl', 'slg_bbs', 'retail_sales_index']],
order=(2,0))
model_fit = model.fit() # fit the model
# make predictions
predictions = model_fit.predict(start=len(train),
end=len(train)+len(test)-1,
exog=test[['slg_adt_ctl', 'slg_bbs', 'retail_sales_index']],
dynamic=False)
当我们生成预测图时，我们现在获得了一些额外的预测能力。
概览
通过使用Scikit-Learn 库，可以考虑使用不同的决策树来预测数据。在此示例中，我们将使用 AdaBoostRegressor，但也可以切换到 RandomForestRegressor 或任何其他可用的树。因此，通过选择树，我们应该意识到去除数据的趋势，通过这种方式，我们举例说明通过分别对数据进行差分和应用对数变换来控制时间序列的均值和方差。
正在准备数据
时间序列有两个基本组成部分，均值和方差。理想情况下，我们希望控制这些组件，对于可变性，我们可以简单地对数据应用对数变换，对于趋势我们可以区分它，我们稍后会看到。
此外，对于这种方法，我们认为实际值 y_t 可以用两个先验值 y_t-1 和 y_t-2 来解释。您可以通过修改 range 函数的输入来使用这些滞后值。
# Load data
tsdf = pd.read_csv('tsdf.csv', sep="\t")
# For simplicity I just take the target variable and the date
tsdf = tsdf[['sales_index', 'dates']]
# Log value
tsdf['log_sales_index'] = np.log(tsdf.sales_index)
# Add previous n-values
for i in range(3):
tsdf[f'sales_index_lag_{i+1}'] = tsdf.sales_index.shift(i+1)
# For simplicity we drop the null values
tsdf.dropna(inplace=True)
tsdf[f'log_sales_index_lag_{i+1}'] = np.log(tsdf[f'sales_index_lag_{i+1}'])
tsdf[f'log_difference_{i+1}'] = tsdf.log_sales_index - tsdf[f'log_sales_index_lag_{i+1}']
数据准备就绪后，我们会得到类似于下图的结果。
数据是否平稳？
要控制时间序列的平均分量，我们应该对数据进行一些差分。为了确定是否需要这一步，我们可以进行单位根检验。有几个对此做出不同假设的测试，可以找到一些单位根测试的列表 here。 为简单起见，我们将考虑 KPSS 测试，我们假设数据是平稳的零假设，特别是，它假设围绕均值或线性趋势保持平稳。
from arch.unitroot import KPSS
# Test for stationary
kpss_test = KPSS(tsdf.sales_index)
# Test summary
print(kpss_test.summary().as_text())
我们看到 P-value = .280 比通常约定的 0.05 大。因此，我们需要对数据应用一阶差分。作为旁注，可以迭代地执行此测试以了解应将多少差异应用于数据。
在下图中，我们看到了原始数据与对数一阶差分的比较，注意时间序列的最后这些值突然发生变化，这似乎是结构性变化，但我们不会深入研究它。如果您想深入了解此主题，布鲁斯·汉森 (Bruce Hansen) 的这些 slides 很有用。
plt.figure(figsize=(12, 6))
plt.subplot(1,2,1)
plt.plot(tsdf.sales_index)
plt.title('Original Time Series')
plt.subplot(1,2,2)
plt.plot(tsdf.log_difference_1)
plt.title('Log first difference Time Series')
决策树模型
正如我们之前所说，我们正在考虑决策树模型，在使用它们时应该注意从时间序列中删除趋势。例如，如果你有上升趋势，tress 不擅长预测下降趋势。在下面的代码示例中，我选择了 AdaBoostRegressor，但您可以自由选择其他树模型。此外，请注意 log_difference_1 被认为是由 log_difference_2 和 log_difference_3 解释的。
Note. Your dataset has other covariates as
aus_avg_rainorslg_adt_ctl, so to consider them for predicting you could apply as well lag values on them.
from sklearn.model_selection import train_test_split
from sklearn.ensemble import AdaBoostRegressor
# Forecast difference of log values
X, Y = tsdf[['log_difference_2', 'log_difference_3']], tsdf['log_difference_1']
# Split in train-test
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, shuffle=False, random_state=0)
# Initialize the estimator
mdl_adaboost = AdaBoostRegressor(n_estimators=500, learning_rate=0.05)
# Fit the data
mdl_adaboost.fit(X_train, Y_train)
# Make predictions
pred = mdl_adaboost.predict(X_test)
test_size = X_test.shape[0]
评估预测
test_size = X_test.shape[0]
plt.plot(list(range(test_size)), np.exp(tsdf.tail(test_size).log_sales_index_lag_1 + pred), label='predicted', color='red')
plt.plot(list(range(test_size)), tsdf.tail(test_size).sales_index, label='real', color='blue')
plt.legend(loc='best')
plt.title('Predicted vs Real with log difference values')
决策树模型似乎准确地预测了真实值。然而，为了评估模型性能，我们应该考虑一个评估指标，可以在这个 article 上找到关于这个主题的很好的介绍，随意选择一个对你的方法更方便的。我将使用 scikit-learn 中的 TimeSeriesSplit 函数，通过平均绝对误差来评估模型的误差。
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error
X, Y = np.array(tsdf[['log_difference_2', 'log_difference_3']]), np.array(tsdf['log_difference_1'])
# Initialize a TimeSeriesSplitter
tscv = TimeSeriesSplit(n_splits=5)
# Retrieve log_sales_index and sales_index to unstransform data
tsdf_log_sales_index = np.array(tsdf.copy().reset_index().log_sales_index_lag_1)
tsdf_sales_index = np.array(tsdf.copy().reset_index().sales_index_lag_1)
# Dict to store metric value at every iteration
metric_iter = {}
for idx, val in enumerate(tscv.split(X)):
train_i, test_i = val
X_train, X_test = X[train_i], X[test_i]
Y_train, Y_test = Y[train_i], Y[test_i]
# Initialize the estimator
mdl_adaboost = AdaBoostRegressor(n_estimators=500, learning_rate=0.05)
# Fit the data
mdl_adaboost.fit(X_train, Y_train)
# Make predictions
pred = mdl_adaboost.predict(X_test)
# Unstransform predictions
pred_untransform = [np.exp(val_test + val_pred) for val_test, val_pred in zip(tsdf_log_sales_index[test_i], pred)]
# Real value
real = tsdf_sales_index[test_i]
# Store metric
metric_iter[f'iter_{idx + 1}'] = mean_absolute_error(real, pred_untransform)
现在我们看到平均 MAE 误差非常低。
print(f'Average MAE error: {np.mean(list(metric_iter.values()))}')
>>> Average MAE error: 17.631090959806535
这有两部分：
- 一个有效的预测代码解决方案。
- 几个话题的简单讨论
解决方案
我将建议一种略有不同且更抽象的方法：使用 Darts which is built on top of scikit-learn. It includes 'A list' libraries (e.g., pandas and NumPy) you'd expect but also a few that you would need to have fairly deep knowledge of the field to think to include (e.g., Torch, Prophet, and Holidays). Furthermore, it has a number of models already included. You can see a bit more for both here。重要说明：特定库是 u8darts，而不是飞镖。它们很相似——并且具有相似的依赖关系——但它们并不相同。
使用它，你可以很容易地开发出一个具有很好结果的预测模型。例如，整个代码，包括合并重命名的列标题（以及删除未重命名的列和重复的行）只有几行。
import pandas as pd
import matplotlib.pyplot as plt
from darts import TimeSeries
from darts.models import ExponentialSmoothing
df = pd.read_csv("../data/tsdf.csv", sep="\t")
df.drop(['Unnamed: 0', 'aus_slg_fmCatl'], axis=1, inplace = True)
df.columns = ['Date', 'adult_cattle(head)','Bulls Bullocks & Steers(head)', 'Cows&Heifers(head)',
'Calves(head)', 'Beef(tonnes)', 'Veal(tonnes)','Adult cattle(kg/head)', 'Calves(kg/head)',
'aust-avg_rain','US/AUS_ExchangeRate', 'sales_index', 'retail_sales_index']
df.drop_duplicates(subset=['Date'], inplace=True)
series = TimeSeries.from_dataframe(df, 'Date', 'Beef(tonnes)').values
train, val = series.split_before(pd.Timestamp('2019-01-01'))
model = ExponentialSmoothing()
model.fit(train)
prediction = model.predict(len(val))
series.plot(label='actual')
prediction.plot(label='forecast', lw=2)
plt.legend()
此外，他们的回购有一个 notebook with additional models and examples
结果
讨论
针对我们收到的评论中的讨论，我认为有几点需要说明。
去除季节性
如果您正在为预测模型拍摄，这是您很少想做的事情。如果您知道每个冬天 home/office 供暖燃料价格都会上涨，但决定将其排除在外 - 您会被很乐意拿走您的钱的商人活活吃掉。
虚拟数据
这里的诀窍是您不知道您的虚拟数据是否与真实数据看起来足够相似。然而，解决方案是一个不同的问题。
“在 Python 中使用 R 库...”
看看r2py
ARMA 和 ARIMA
... 不太可能产生丰硕的结果，部分原因是误差项假定误差为 IID，但市场存在偏差。指数平滑效果很好，因为它替代了
“但是当我 运行 它反对 sales_index...”
sales_index 将由几个自变量组成。正如您所期望的那样，肉类生产（奇怪的是与降雨量呈负相关）和汇率。但是，您的模型中不包括其他国家/地区的生产数据 (100% pure organic, grass feed, 'AAA' Alberta beef)、税收、国内牛肉生产、运输等。这就是为什么您可能不想与拥有这两者的专业交易员对抗宽客和认真的领域专业知识。最后，我要指出的是，在时间上没有任何迹象表明这是正常的还是不正常的。
从 2015 年到 2017 年似乎有 1.5 年的上升期，随后一直下降到 2019 年。似乎 正在重复，但 2020 年的突然变化导致人们相信这是一种失常。但数据太小，很难判断任何有效性。在这种情况下，趋势线或通道会更好。