⽤Python进⾏时间序列预测的7种⽅法
数据准备
数据集（JetRail⾼铁的乘客数量）.
假设要解决⼀个时序问题：根据过往两年的数据（2012 年 8 ⽉⾄ 2014 年 8⽉），需要⽤这些数据预测接下来 7 个⽉的乘客数量。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
df = pd.read_csv('train.csv')
df.head()
df.shape
依照上⾯的代码，我们获得了 2012-2014 年两年每个⼩时的乘客数量。为了解释每种⽅法的不同之处，以每天为单位构造和聚合了⼀个数据集。
从 2012 年 8 ⽉- 2013 年 12 ⽉的数据中构造⼀个数据集。
创建 train and test ⽂件⽤于建模。前 14 个⽉（ 2012 年 8 ⽉- 2013 年 10 ⽉）⽤作训练数据，后两个⽉（2013 年 11 ⽉ – 2013 年 12⽉）⽤作测试数据。
以每天为单位聚合数据集。
import pandas as pd
import matplotlib.pyplot as plt
# 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')  # 4位年⽤Y，2位年⽤y
df.index = df['Timestamp']
df = df.resample('D').mean() #按天采样，计算均值
train['Timestamp'] = pd.to_datetime(train['Datetime'], format='%d-%m-%Y %H:%M')python index函数
train.index = train['Timestamp']
train = sample('D').mean() #
test['Timestamp'] = pd.to_datetime(test['Datetime'], format='%d-%m-%Y %H:%M')
test.index = test['Timestamp']
test = sample('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()
我们将数据可视化（训练数据和测试数据⼀起），从⽽得知在⼀段时间内数据是如何变化的。
⽅法1：朴素法
假设 y 轴表⽰物品的价格，x 轴表⽰时间（天）
如果数据集在⼀段时间内都很稳定，我们想预测第⼆天的价格，可以取前⾯⼀天的价格，预测第⼆天的值。这种假设第⼀个预测点和上⼀个观察点相等的预测⽅法就叫朴素法。即 yt+1^=ytyt+1^=yt
dd = np.asarray(train['Count'])
y_hat = py()
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()
朴素法并不适合变化很⼤的数据集，最适合稳定性很⾼的数据集。我们计算下均⽅根误差，检查模型在测试数据集上的准确率：
ics import mean_squared_error
from math import sqrt
rms = sqrt(mean_squared_error(test['Count'], y_hat['naive']))
print(rms)
最终均⽅误差RMS为：43.91640614391676
⽅法2：简单平均法
物品价格会随机上涨和下跌，平均价格会保持⼀致。我们经常会遇到⼀些数据集，虽然在⼀定时期内出现⼩幅变动，但每个时间段的平均值确实保持不变。这种情况下，我们可以预测出第⼆天的价格⼤致和过去天数的价格平均值⼀致。这种将预期值等同于之前所有观测点的平均值的预测⽅法就叫简单平均法。即y^x+1=1x∑i=1xyiy^x+1=1x∑i=1xyi
y_hat_avg = py()
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()
⽅法3：移动平均法
物品价格在⼀段时间内⼤幅上涨，但后来⼜趋于平稳。我们也经常会遇到这种数据集，⽐如价格或销售额某段时间⼤幅上升或下降。如果我们这时⽤之前的简单平均法，就得使⽤所有先前数据的平均值，但在这⾥使⽤之前的所有数据是说不通的，因为⽤开始阶段的价格值会⼤幅影响接下来⽇期的预测值。因此，我们只取最近⼏个时期的价格平均值。很明显这⾥的逻辑是只有最近的值最要紧。这种⽤某些窗⼝期计算平均值的预测⽅法就叫移动平均法。
计算移动平均值涉及到⼀个有时被称为“滑动窗⼝”的⼤⼩值p。使⽤简单的移动平均模型，我们可以根据之前数值的固定有限数p的平均值预测某个时序中的下⼀个值。这样，对于所有的 i>pi>p：
y^l=1p(yi−1+yi−2+...+yi−p)y^l=1p(yi−1+yi−2+...+yi−p)
在上⽂移动平均法可以看到，我们对“p”中的观察值赋予了同样的权重。但是我们可能遇到⼀些情况，⽐如“p”中每个观察值会以不同的⽅式影响预测结果。将过去观察值赋予不同权重的⽅法就叫做加权移动平均法。加权移动平均法其实还是⼀种移动平均法，只是“滑动窗⼝期”内的值被赋予不同的权重，通常来讲，最近时间点的值发挥的作⽤更⼤了。
即y^l=1m(w1∗yi−1+w2∗yi−2+...+wm∗yi−m)y^l=1m(w1∗yi−1+w2∗yi−2+...+wm∗yi−m)
这种⽅法并⾮选择⼀个窗⼝期的值，⽽是需要⼀列权重值（相加后为1）。例如，如果我们选择[0.40, 0.25, 0.20, 0.15]作为权值，我们会为最近的4个时间点分别赋给40%，25%，20%和15%的权重。
⽅法4：简单指数法
我们注意到简单平均法和加权移动平均法在选取时间点的思路上存在较⼤的差异。我们就需要在这两种⽅法之间取⼀个折中的⽅法，在将所有数据考虑在内的同时也能给数据赋予不同⾮权重。例如，相⽐更早时期内的观测值，它会给近期的观测值赋予更⼤的权重。按照这种原则⼯作的⽅法就叫做简单指数平滑法。它通过加权平均值计算出预测值，其中权重随着观测值从早期到晚期的变化呈指数级下降，最⼩的权重和最早的观测值相关：
y^T+1=αyT+α(1−α)yT−1+α(1−α)2yTt−2+...y^T+1=αyT+α(1−α)yT−1+α(1−α)2yTt−2+...
其中0≤α≤1是平滑参数。对时间点T+1的单步预测值是时序y1,…,yTy1,…,yT的所有观测值的加权平均数。权重下降的速率由参数α控制，预测值y^xy^x是α⋅ytα⋅yt与(1−α)⋅y^x(1−α)⋅y^x的和。
因此，它可以写为：
y^T+1=αyT+(1−α)yT−1y^T+1=αyT+(1−α)yT−1
所以本质上，我们是⽤两个权重α和1−α得到⼀个加权移动平均值，让表达式呈递进形式。
from statsmodels.tsa.api import SimpleExpSmoothing
y_hat_avg = py()
fit = SimpleExpSmoothing(np.asarray(train['Count'])).fit(smoothing_level=0.6, optimized=False)
y_hat_avg['SES'] = fit.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()
模型中使⽤的α值为0.6，我们可以⽤测试集继续调整参数以⽣成⼀个更好的模型。
⽅法5：霍尔特(Holt)线性趋势法
假设y轴表⽰某个物品的价格，x轴表⽰时间（天）。
如果物品的价格是不断上涨的（见上图），我们上⾯的⽅法并没有考虑这种趋势，即我们在⼀段时间内观察到的价格的总体模式。
每个时序数据集可以分解为相应的⼏个部分：趋势（Trend），季节性(Seasonal)和残差(Residual)。任何呈现某种趋势的数据集都可以⽤霍尔特线性趋势法⽤于预测。
import statsmodels.api as sm
sm.tsa.seasonal_decompose(train['Count']).plot()
result = sm.tsa.stattools.adfuller(train['Count'])
plt.show()
我们从图中可以看出，该数据集呈上升趋势。因此我们可以⽤霍尔特线性趋势法预测未来价格。该算法包含三个⽅程：⼀个⽔平⽅程，⼀个趋势⽅程，⼀个⽅程将⼆者相加以得到预测值y^y^：
我们在上⾯算法中预测的值称为⽔平（level）。正如简单指数平滑⼀样，这⾥的⽔平⽅程显⽰它是观测值和样本内单步预测值的加权平均数，趋势⽅程显⽰它是根据ℓ(t)−ℓ(t−1) 和之前的预测趋势 b(t−1) 在时间t处的预测趋势的加权平均值。
我们将这两个⽅程相加，得出⼀个预测函数。我们也可以将两者相乘⽽不是相加得到⼀个乘法预测⽅程。当趋势呈线性增加和下降时，我们⽤相加得到的⽅程；当趋势呈指数级增加或下降时，我们⽤相乘得到的⽅程。实践操作显⽰，⽤相乘得到的⽅程，预测结果会更稳定，但⽤相加得到的⽅程，更容易理解。
from statsmodels.tsa.api import Holt
y_hat_avg = py()
fit = Holt(np.asarray(train['Count'])).fit(smoothing_level=0.3, smoothing_slope=0.1)
y_hat_avg['Holt_linear'] = fit.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()
这种⽅法能够准确地显⽰出趋势，因此⽐前⾯的⼏种模型效果更好。如果调整⼀下参数，结果会更好。
⽅法6：Holt-Winters季节性预测模型
在应⽤这种算法前，我们先介绍⼀个新术语。假如有家酒店坐落在半⼭腰上，夏季的时候⽣意很好，顾
客很多，但每年其余时间顾客很少。因此，每年夏季的收⼊会远⾼于其它季节，⽽且每年都是这样，那么这种重复现象叫做“季节性”（Seasonality）。如果数据集在⼀定时间段内的固定区间内呈现相似的模式，那么该数据集就具有季节性。
我们之前讨论的5种模型在预测时并没有考虑到数据集的季节性，因此我们需要⼀种能考虑这种因素的⽅法。应⽤到这种情况下的算法就叫做Holt-Winters季节性预测模型，它是⼀种三次指数平滑预测，其背后的理念就是除了⽔平和趋势外，还将指数平滑应⽤到季节分量上。
Holt-Winters季节性预测模型由预测函数和三次平滑函数——⼀个是⽔平函数ℓt，⼀个是趋势函数bt，⼀个是季节分量 st，以及平滑参数α,β和γ。
其中 s 为季节循环的长度，0≤α≤ 1, 0 ≤β≤ 1 ， 0≤γ≤ 1。⽔平函数为季节性调整的观测值和时间点t处⾮季节预测之间的加权平均值。趋势函数和霍尔特线性⽅法中的含义相同。季节函数为当前季节指数和去年同⼀季节的季节性指数之间的加权平均值。在本算法，我们同样可以⽤相加和相乘的⽅法。当季节性变化⼤致相同时，优先选择相加⽅法，⽽当季节变化的幅度与各时间段的⽔平成正⽐时，优先选择相乘的⽅法。
from statsmodels.tsa.api import ExponentialSmoothing
y_hat_avg = py()
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()
我们可以看到趋势和季节性的预测准确度都很⾼。我们选择了 seasonal_period = 7作为每周重复的数据。也可以调整其它其它参数，我在搭建这个模型的时候⽤的是默认参数。你可以试着调整参数来优化模型。
⽅法7：⾃回归移动平均模型（ARIMA）
另⼀个场景的时序模型是⾃回归移动平均模型（ARIMA）。指数平滑模型都是基于数据中的趋势和季节性的描述，⽽⾃回归移动平均模型的⽬标是描述数据中彼此之间的关系。ARIMA的⼀个优化版就是季节性ARIMA。它像Holt-Winters季节性预测模型⼀样，也把数据集的季节性考虑在内。
import statsmodels.api as sm
y_hat_avg = py()
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()
我们可以看到使⽤季节性 ARIMA 的效果和Holt-Winters差不多。我们根据 ACF（⾃相关函数）和 PACF（偏⾃相关）图选择参数。如果你为 ARIMA 模型选择参数时遇到了困难，可以⽤ R 语⾔中的 auto.arima。
最后，我们将这⼏种模型的准确度⽐较⼀下：
后话
建议你在解决问题时，可以依次试试这⼏种模型，看看哪个效果最好。我们从上⽂也知道，数据集不同，每种模型的效果都有可能优于其它模型。因此，如果⼀个模型在某个数据集上效果很好，并不代表它在所有数据集上都⽐其它模型好。
参考链接：
1.
2.
3.