有什么方法可以用python中的scikit-learn来预测每月的时间序列？

Question

我想product' sales_index通过使用每月时间序列中的多个特征进行预测。一开始，我开始使用ARMA,ARIMA来做到这一点，但输出对我来说不是很令人满意。在我的尝试中，我只是使用dates和sales列进行预测，输出对我来说并不现实。我想我应该包括更多的特征列来预测sales_index列。但是，我想知道是否有任何方法可以通过使用每月时间序列中的多个特征来进行此预测。我没有使用scikit-learn. 任何人都可以指出我这样做的任何可能方法吗？任何可能的想法？

我尝试使用 ARMA/ARIMA：

这是关于这个要点的可重现的月度时间序列数据，这是我目前的尝试：

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？任何可能的想法？谢谢

Answer 1

概述

通过使用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}']

一旦我们的数据准备好，我们就会得到类似于下图的结果。

数据是静止的吗？

为了控制时间序列的平均分量，我们应该对数据进行一些差分。为了确定是否需要这一步，我们可以进行单位根检验。有几个测试对此做出了不同的假设，可以在此处找到一些单位根测试的列表。为简单起见，我们将考虑 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 的这些幻灯片很有用。

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')

决策树模型

正如我们之前所说，我们正在考虑决策树模型，在使用它们时，应该注意从时间序列中删除趋势。例如，如果您有上升趋势，则 tres 不擅长预测下降趋势。在下面的代码示例中，我选择了AdaBoostRegressor，但您可以自由选择其他树模型。另外，请注意 log_difference_1被认为是由log_difference_2和解释的log_difference_3。

注意。您的数据集有其他协变量aus_avg_rain或slg_adt_ctl，因此要考虑使用它们进行预测，您也可以对它们应用滞后值。

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')

似乎决策树模型正在准确预测真实值。但是，为了评估模型性能，我们应该考虑评估指标，可以在本文中找到有关该主题的良好介绍，您可以随意选择对您的方法更方便的方法。我将使用TimeSeriesSplitfrom 函数scikit-learn通过平均绝对误差来评估模型的误差。

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

Answer 2

这有两个部分：

1. 一个有效的预测代码解决方案。
2. 几个话题的简短讨论

解决方案

我将建议一种稍微不同且更抽象的方法：使用构建在 scikit-learn 之上的Darts 。它包括您期望的“A 列表”库（例如 pandas 和 NumPy），但也包括一些您需要对该领域有相当深入了解才能考虑包含的库（例如 Torch、Prophet 和 Holidays）。此外，它已经包含了许多模型。您可以在这里查看两者的更多信息。一个重要的注意事项：特定的库是 u8darts，而不是 darts。它们很相似——并且具有相似的依赖关系——但它们并不相同。

使用它，您可以轻松开发出具有良好结果的预测模型。例如，整个代码，包括合并重命名的列标题（以及删除未重命名的列和重复的行），只有几行。

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

此外，他们的存储库还有一个带有其他模型和示例的笔记本

结果

讨论

针对我们评论中的讨论，我认为需要指出几点。

去除季节性

如果您正在研究预测模型，那么这是您很少想做的事情。如果您知道每年冬天家庭/办公室供暖的燃料价格都会上涨，但决定将其排除在外，那么您就会被很乐意拿走您的钱的商人活活吃掉。

虚拟数据

这里的技巧是您不知道您的虚拟数据是否看起来与真实数据足够相似。然而，解决这个问题的是另一个问题。

“在 Python 中使用 R 库...”

看看r2py

ARMA 和 ARIMA

...不太可能产生丰硕的结果，部分原因是误差项假设误差为独立同分布，但市场存在偏差。指数平滑效果很好，因为它替代了

“但是当我针对 sales_index 运行它时......”

sales_index 将由几个自变量组成。如您所料，肉类产量（奇怪的是与降雨量成反比）和汇率。但是，您的模型中不包含其他国家的生产数据（100% 纯有机、草饲料、“AAA”阿尔伯塔牛肉）、税收、国内牛肉生产、运输等。这就是为什么您可能不想要与拥有量化分析师和专业领域专业知识的专业交易员进行对抗。最后，我要指出的是，没有任何时间可以表明这是正常还是不正常。

从 2015 年到 2017 年，似乎出现了 1.5 年的上升趋势，随后一直下降到 2019 年。这种情况似乎正在重复，但 2020 年的突然变化让人相信这是一种失常。但数据太小，很难判断其有效性。在这种情况下，趋势线或通道会发挥更好的作用。