用Python进行时序预测的7种方法

前言

我们很多人现在应该经常听到各种关于市场投资方面的热词，比如虚拟货币，可能不少人还买了这种币。但是，投资这种波动性很大的数字货币安全吗？我们怎样能确保投资这些虚拟货币在未来一定会获得收益？实际上我们没法保证一定会有很好的收益，但可以根据之前的虚拟货币价格估算出近似价值。时序模型就是预测这些值的一种方式。

除了虚拟货币，时序预测还可以应用在很多重要的领域，比如预测销售额、预测呼叫中心的来电数量、预测太阳活动和海洋潮汐、预测股市行为等等。

假设你现在是一家酒店的经理，想预测接下来一年内会有多少顾客，从而相应地调整酒店的进货量，并对酒店的营业额有个合理的预估。

根据之前每年、每月和每天的数据，你就可以用时序预测的方法，对未来一年的顾客量有个接近的预估值。而顾客量的预测值能够帮助酒店合理管理资源和制定计划。

在本文，我们会学习多种时序预测方法，并用同一个数据集对它们进行比较。我们会逐一认识7种方法，以及如何用这些方法优化结果。

我们开始吧！

文章目录

了解问题说明和数据集

• 安装程序库（statsmodels）
• 方法1——先以朴素法开始
• 方法2——简单平均数
• 方法3——移动平均数
• 方法4——指数平滑法
• 方法5——霍尔特线性趋势预测
• 方法6——Holt-Winters季节性预测模型
• 方法7——自回归移动平均模型

了解问题说明和数据集

我们假设要解决一个时序问题：预测搭乘高铁服务商 JetRail 旗下高铁的乘客数量。我们手头有过往两年的数据（2012 年 8 月至 2014 年 8月），需要用这些数据预测接下来 7 个月的乘客数量。

我们首先下载需要的数据集（https://datahack.analyticsvidhya.com/contest/practice-problem-time-series-2/），在本文，我只用到火车数据集。

import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt 
#Importing data
df = pd.read_csv('train.csv')
#Printing head
df.head()

#Printing tail
df.tail()

如上面的输出语句所示，我们获得了 2012-2014 年两年每个小时的乘客数量，然后需要预测未来的乘客数量。

在本文，为了解释每种方法的不同之处，我以每天为单位构造和聚合了一个数据集。

从 2012 年 8 月- 2013 年 12 月的数据中构造一个数据集。

创建 train and test 文件用于建模。前 14 个月（ 2012 年 8 月- 2013 年 10 月）用作训练数据，后两个月（2013 年 11 月 - 2013 年 12 月）用作测试数据。

以每天为单位聚合数据集。

#Subsetting the dataset
#Index 11856 marks the end of year 2013
df = pd.read_csv('train.csv', nrows = 11856)
 
#Creating train and test set 
#Index 10392 marks the end of October 2013 
train=df[0:10392] 
test=df[10392:]
 
#Aggregating the dataset at daily level
df.Timestamp = pd.to_datetime(df.Datetime,format='%d-%m-%Y %H:%M') 
df.index = df.Timestamp 
df = df.resample('D').mean()
train.Timestamp = pd.to_datetime(train.Datetime,format='%d-%m-%Y %H:%M') 
train.index = train.Timestamp 
train = train.resample('D').mean() 
test.Timestamp = pd.to_datetime(test.Datetime,format='%d-%m-%Y %H:%M') 
test.index = test.Timestamp 
test = test.resample('D').mean()

我们将数据可视化（训练数据和测试数据一起），从而得知在一段时间内数据是如何变化的。

#Plotting data
train.Count.plot(figsize=(15,8), title= 'Daily Ridership', fontsize=14)
test.Count.plot(figsize=(15,8), title= 'Daily Ridership', fontsize=14)
plt.show()

安装程序库（statsmodels）

我用来进行时序预测的这个程序库叫 statsmodels。在尝试我们本文介绍的几种方法前，你需要安装这个程序库。可能你的 Python 环境中之前已经安装了 statesmodels，但并不支持我们的预测方法。我们需要将它从程序库中克隆下来，用源代码安装，步骤如下：

• 使用 pip freeze 检查你的 Python 环境中是否已经安装了 statesmodels
• 如果已经存在，用 conda remove statsmodels 将其移除
• 用 git clone git://http://github.com/statsmodels/statsmodels.git 克隆 statesmodels 程序库。在克隆前，用 git init 初始化 Git
• 用cd statsmodels将目录更改为statesmodels
• 用python setup.py build创建设置文件
• 用python setup.py install安装它
• 退出bash/terminal

在你本地环境下重启 bash/terminal，打开 Python，执行 from statsmodels.tsa.api import ExponentialSmoothing 进行验证

方法1：先以朴素法开始

考虑一下下图。我们假设 y 轴表示一种币的价格，x 轴表示时间（天）。

我们可以从图表中得知币值从一开始就很稳定。我们经常遇到数据集在一段时间内始终很稳定。如果我们想预测第二天的价格，可以简单的取前面一天的价格，预测第二天的值。这种假设第一个预测点和上一个观察点相等的预测方法就叫朴素法。

现在，我们用朴素法预测测试集的价格。

dd= np.asarray(train.Count)
y_hat = test.copy()
y_hat['naive'] = dd[len(dd)-1]
plt.figure(figsize=(12,8))
plt.plot(train.index, train['Count'], label='Train')
plt.plot(test.index,test['Count'], label='Test')
plt.plot(y_hat.index,y_hat['naive'], label='Naive Forecast')
plt.legend(loc='best')
plt.title("Naive Forecast")
plt.show()

我们下面计算均方根误差，检查模型在测试数据集上的准确率。

from sklearn.metrics import mean_squared_error
from math import sqrt
rms = sqrt(mean_squared_error(test.Count, y_hat.naive))
print(rms)
 
RMSE = 43.9164061439

从均方根误差值以及上面的图表可以看出，朴素法并不适合变化很大的数据集，最适合稳定性很高的数据集。我们还可以用不同的方法优化结果。下面我们试试其它方法。

方法2：简单平均法

考虑一下下图。我们假设y轴表示某个货币的价格，x轴表示时间（天）。

从上图我们可以看到货币价格会随机上涨和下跌，价格线出现小小的切口，这样以来平均价格会保持一致。我们经常会遇到一些数据集，虽然在一定时期内出现小幅变动，但每个时间段的平均值确实保持不变。这种情况下，我们可以预测出第二天的价格大致和过去天数的价格平均值一致。

这种将预期值等同于之前所有观测点的平均值的预测方法就叫简单平均法。

我们用之前全部已知的值计算出它们的平均值，将它作为要预测的下一个值。当然这不会很准确，但已经很接近了。这种预测方法在某些情况下效果是最好的。

y_hat_avg = test.copy()
y_hat_avg['avg_forecast'] = train['Count'].mean()
plt.figure(figsize=(12,8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['avg_forecast'], label='Average Forecast')
plt.legend(loc='best')
plt.show()

我们现在计算均方根误差值，检查模型的准确率。

rms = sqrt(mean_squared_error(test.Count, y_hat_avg.avg_forecast))
print(rms)
 
RMSE = 109.545990803

我们可以看到这种模型并没有改善准确率。因此我们可以从中推断出当每个时间段的平均值保持不变时，这种方法的效果才能达到最好。虽然朴素法的准确率高于简单平均法，但这并不意味着朴素法在所有的数据集上都比简单平均法好。

方法3——移动平均法

看一下下图，我们还是假设y轴为货币价格，x轴为时间（天）。

我们从图中可以知道货币价格在一段时间内大幅上涨，但后来又趋于平稳。我们也经常会遇到这种数据集，比如价格或销售额某段时间大幅上升或下降。如果我们这时用之前的简单平均法，就得使用所有先前数据的平均值，但在这里使用之前的所有数据是说不通的，因为用开始阶段的价格值会大幅影响接下来日期的预测值。

因此，我们只取最近几个时期的价格平均值。很明显这里的逻辑是只有最近的值最要紧。这种用某些窗口期计算平均值的预测方法就叫移动平均法。计算移动平均值涉及到一个有时被称为“滑动窗口”的大小值n。使用简单的移动平均模型，我们可以根据之前数值的固定有限数p的平均值预测某个时序中的下一个值。这样，对于所有的 i > p

移动平均法实际上很有效，特别是当你为时序选择了正确的p值时。

y_hat_avg = test.copy()
y_hat_avg['moving_avg_forecast'] = train['Count'].rolling(60).mean().iloc[-1]
plt.figure(figsize=(16,8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['moving_avg_forecast'], label='Moving Average Forecast')
plt.legend(loc='best')
plt.show()

我们只选择过去两个月的数据。现在计算均方根误差值，检查模型的准确度。

rms = sqrt(mean_squared_error(test.Count, y_hat_avg.moving_avg_forecast))
print(rms)
RMSE = 46.7284072511

我们可以看到，对于这个数据集，朴素法比简单平均法和移动平均法的表现要好。下面我们试一试简单指数平滑法，它比移动平均法的一个进步之处就是相当于对移动平均法进行了加权。在上文移动平均法可以看到，我们对“n”中的观察值赋予了同样的权重。但是我们可能遇到一些情况，比如“n”中每个观察值会以不同的方式影响预测结果。将过去观察值赋予不同权重的方法就叫做加权移动平均法。

加权移动平均法其实还是一种移动平均法，只是“滑动窗口期”内的值被赋予不同的权重，通常来讲，最近时间点的值发挥的作用更大了。

这种方法并非选择一个窗口期的值，而是需要一列权重值（相加后为1）。例如，如果我们选择[0.40, 0.25, 0.20, 0.15]作为权值，我们会为最近的4个时间点分别赋给40%，25%，20%和15%的权重。

方法4——简单指数平滑法

在认识上面这几种方法后，我们可以注意到简单平均法和加权移动平均法依靠的时间点完全相反。那么我们就需要在这两种方法之间取一个折中的方法，在将所有数据考虑在内的同时也能给数据赋予不同非权重。例如，相比更早时期内的观测值，它会给近期的观测值赋予更大的权重。按照这种原则工作的方法就叫做简单指数平滑法。它通过加权平均值计算出预测值，其中权重随着观测值从早期到晚期的变化呈指数级下降，最小的权重和最早的观测值相关：

其中0≤ α ≤1是平滑参数。对时间点T+1的单步预测值是时序y1,···，yT的所有观测值的加权平均数。权重下降的速率由参数α控制，预测值 ŷx 是α⋅y 与(1−α)·ŷ t-1 的和。

因此，它可以写为：

所以本质上，我们是用两个权重 α 和 1−α 得到一个加权移动平均值。

我们可以看到 ŷ x−1 和1−α 相乘，让表达式呈递进形式，这也是该方法被称为“指数”的原因。时间 t+1 处的预测值为最近观测值 yt 和最近预测值 ŷ t|t−1 之间的加权平均值。

from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt
y_hat_avg = test.copy()
fit2 = SimpleExpSmoothing(np.asarray(train['Count'])).fit(smoothing_level=0.6,optimized=False)
y_hat_avg['SES'] = fit2.forecast(len(test))
plt.figure(figsize=(16,8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['SES'], label='SES')
plt.legend(loc='best')
plt.show()

我们下面计算一下均方根误差，验证一下模型的准确度。

rms = sqrt(mean_squared_error(test.Count, y_hat_avg.SES))
print(rms)
 
RMSE = 43.3576252252

我们可以看到 α 为 0.6 时，简单指数平滑模型生成了迄今为止最好的结果。我们可以用测试集继续调整参数以生成一个更好的模型。

方法5——霍尔特线性趋势法

我们现在已经认识了几种预测方法，但是我们也可以看到上面这些模型在波动性较大的数据集上表现并不是很好。

如果货币的价格是不断上涨的（见上图），我们上面的方法并没有考虑这种趋势，即我们在一段时间内观察到的价格的总体模式。在上图例子中，我们可以看到货币的价格呈上涨趋势。

虽然上面这些方法都可以应用于这种趋势，但我们仍需要一种方法可以在无需假设的情况下，准确预测出价格趋势。这种考虑到数据集变化趋势的方法就叫做霍尔特线性趋势法。

每个时序数据集可以分解为相应的几个部分：趋势，季节性和残差。任何呈现某种趋势的数据集都可以用霍尔特线性趋势法用于预测。

import statsmodels.api as sm
sm.tsa.seasonal_decompose(train.Count).plot()
result = sm.tsa.stattools.adfuller(train.Count)
plt.show()

我们从图中可以看出，该数据集呈上升趋势。因此我们可以用霍尔特线性趋势法预测未来价格。

该算法包含三个方程：一个水平方程，一个趋势方程，一个方程将二者相加以得到预测值 ŷ ：

我们在上面算法中预测的值称为水平（level）。在上面三个方程中，你可以看到我们将水平和趋势相加以生成预测方程。

正如简单指数平滑一样，这里的水平方程显示它是观测值和样本内单步预测值的加权平均数，趋势方程显示它是根据 ℓ(t)−ℓ(t−1) 和之前的预测趋势 b(t−1) 在时间t处的预测趋势的加权平均值。

我们将这两个方程相加，得出一个预测函数。我们也可以将两者相乘而不是相加得到一个乘法预测方程。当趋势呈线性增加和下降时，我们用相加得到的方程；当趋势呈指数级增加或下降时，我们用相乘得到的方程。实践操作显示，用相乘得到的方程，预测结果会更稳定，但用相加得到的方程，更容易理解。

y_hat_avg = test.copy()
 
fit1 = Holt(np.asarray(train['Count'])).fit(smoothing_level = 0.3,smoothing_slope = 0.1)
y_hat_avg['Holt_linear'] = fit1.forecast(len(test))
 
plt.figure(figsize=(16,8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['Holt_linear'], label='Holt_linear')
plt.legend(loc='best')
plt.show()

我们计算一下均方根误差，检查模型的准确率：

rms = sqrt(mean_squared_error(test.Count, y_hat_avg.Holt_linear))
print(rms)
RMSE = 43.0562596115

我们可以看到这种方法能够准确地显示出趋势，因此比前面的几种模型效果更好。如果调整一下参数，结果会更好。

方法6——Holt-Winters季节性预测模型

在应用这种算法前，我们先介绍一个新术语。假如有家酒店坐落在半山腰上，夏季的时候生意很好，顾客很多，但每年其余时间顾客很少。因此，每年夏季的收入会远高于其它季节，而且每年都是这样，那么这种重复现象叫做“季节性”（Seasonality）。如果数据集在一定时间段内的固定区间内呈现相似的模式，那么该数据集就具有季节性。

我们之前讨论的5种模型在预测时并没有考虑到数据集的季节性，因此我们需要一种能考虑这种因素的方法。应用到这种情况下的算法就叫做Holt-Winters季节性预测模型，它是一种三次指数平滑预测，其背后的理念就是除了水平和趋势外，还将指数平滑应用到季节分量上。

Holt-Winters季节性预测模型由预测函数和三次平滑函数——一个是水平函数ℓt，一个是趋势函数bt，一个是季节分量 st，以及平滑参数 α, β 和 γ。

其中 s 为季节循环的长度，0 ≤ α ≤ 1, 0 ≤ β ≤ 1 ， 0 ≤ γ ≤ 1.

水平函数为季节性调整的观测值和时间点t处非季节预测之间的加权平均值。趋势函数和霍尔特线性方法中的含义相同。季节函数为当前季节指数和去年同一季节的季节性指数之间的加权平均值。

在本算法，我们同样可以用相加和相乘的方法。当季节性变化大致相同时，优先选择相加方法，而当季节变化的幅度与各时间段的水平成正比时，优先选择相乘的方法。

y_hat_avg = test.copy()
fit1 = ExponentialSmoothing(np.asarray(train['Count']) ,seasonal_periods=7 ,trend='add', seasonal='add',).fit()
y_hat_avg['Holt_Winter'] = fit1.forecast(len(test))
plt.figure(figsize=(16,8))
plt.plot( train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['Holt_Winter'], label='Holt_Winter')
plt.legend(loc='best')
plt.show()

我们现在计算一下均方根误差，看看模型的准确度。

rms = sqrt(mean_squared_error(test.Count, y_hat_avg.Holt_Winter))
print(rms)
 
RMSE = 23.9614925662

我们可以看到趋势和季节性的预测准确度都很高。我们选择了 seasonal_period = 7作为每周重复的数据。也可以调整其它其它参数，我在搭建这个模型的时候用的是默认参数。你可以试着调整参数来优化模型。

方法7——自回归移动平均模型

另一个场景的时序模型是自回归移动平均模型（ARIMA）。指数平滑模型都是基于数据中的趋势和季节性的描述，而自回归移动平均模型的目标是描述数据中彼此之间的关系。ARIMA的一个优化版就是季节性ARIMA。它像Holt-Winters季节性预测模型一样，也把数据集的季节性考虑在内。

y_hat_avg = test.copy()
fit1 = sm.tsa.statespace.SARIMAX(train.Count, order=(2, 1, 4),seasonal_order=(0,1,1,7)).fit()
y_hat_avg['SARIMA'] = fit1.predict(start="2013-11-1", end="2013-12-31", dynamic=True)
plt.figure(figsize=(16,8))
plt.plot( train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['SARIMA'], label='SARIMA')
plt.legend(loc='best')
plt.show()

下面我们计算一下均方根误差，看看模型的准确度。

rms = sqrt(mean_squared_error(test.Count, y_hat_avg.SARIMA))
print(rms)
 
RMSE = 26.035582877

我们可以看到使用季节性 ARIMA 的效果个 Holt-Winters 一样好。我们根据 ACF 和 PACF 图选择参数。如果你为 ARIMA 模型选择参数时遇到了困难，可以用 R 语言中的 auto.arima。

最后，我们将这几种模型的准确度比较一下：

后话

希望本文对你有所帮助，在解决时许问题的时候能从容以对。我建议你在解决问题时，可以依次试试这几种模型，看看哪个效果最好。

我们从上文也知道，数据集不同，每种模型的效果都有可能优于其它模型。因此，如果一个模型在某个数据集上效果很好，并不代表它在所有数据集上都比其它模型好。

参考资料：
https://www.analyticsvidhya.com/blog/2018/02/time-series-forecasting-methods/