【译】时间序列预测初学者指南

这篇文章是《 基于R语言的时间序列建模完整教程 》的后续文章，不同的是本文采用 Python 来进行讲解。本文在原文基础上删除和修改了部分内容，如遇到不明白的，请查看 原文 。

在 pandas 加载时间序列数据

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
 
data = pd.read_csv('AirPassengers.csv', parse_dates=['TravelDate'], index_col='TravelDate')
print(data.head())
 
ts = data['Passengers']
ts.head(10)
 
plt.plot(ts)

关于Pandas如何加载数据的，请查看： 如何使用Pandas读取Excel和CSV文件数据

检验时间序列稳定性

前面的文章中已经讲到稳定性的评判标准主要是均值、方差和协商差。具体可以通过下述两种方法进行测试：

1、绘制滚动统计：绘制移动平均数和移动方差，观察它是否随着时间变化。

2、Dickey-Fuller检验：测试结果由测试统计量和一些置信区间的临界值组成。如果“测试统计量”少于“临界值”，并认为序列是稳定的。

以下方法定义了以上两种方式(注意，为了保持单位和平均数相似，这里采用了标准差来代替方差)

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
    
    #Determing rolling statistics
    rolmean = timeseries.rolling(window=12).mean()
    rolstd = timeseries.rolling(window=12).std()
 
    #Plot rolling statistics:
    orig = plt.plot(timeseries, color='blue',label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    #Perform Dickey-Fuller test:
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value (%s)'%key] = value
    print(dfoutput)
test_stationarity(ts)

这里解释一下，DF测试的值怎么看：

• Test Statistic的值如果比Critical Value (5%)小则满足稳定性需求
• p-value越低（理论上需要低于05）证明序列越稳定。

这里这个结果表明这些序列很不稳定，所以接下来考虑如何处理数据，使得序列相对稳定。

使时间序列平稳

有两个主要的因素导致序列不稳定：

• 趋势 Trend
• 季节性 Seasonality

消除趋势

消除趋势的第一个方法是转换。在本例中我们可以清楚地看到有一个显著的趋势。所以我们可以通过变换，惩罚较高值而不是较小值。这可以采用取对数log,平方根,立方跟等。让我们简单在这儿转换一个对数。

ts_log = np.log(ts)
plt.plot(ts_log)

这里我们可以明显看到一个上升的趋势，但是混杂在噪音当中，所以需要去除杂音。这里简单平滑一下数据。平滑的窗口取值12，因为一年有12个月。

在这个简单的例子中,很容易看到一个向前的数据趋势。但是它表现的不是很直观。所以我们可以使用一些技术来估计或对这个趋势建模,然后将它从序列中删除。这里有很多方法,最常用的有:

• 聚合-取一段时间的平均值（月/周平均值）
• 平滑-取滚动平均数
• 多项式回归分析-适合的回归模型

我在这儿讨论将平滑,你也应该尝试其他可以解决的问题的技术。平滑是指采取滚动估计,即考虑过去的几个实例。有各种方法可以解决这些问题，但我将主要讨论以下两个。

移动平均数

在这个方法中，根据时间序列的频率采用“K”连续值的平均数。我们可以采用过去一年的平均数，即过去12个月的平均数。

moving_avg = ts_log.rolling(12).mean()
plt.plot(ts_log)
plt.plot(moving_avg, color='red')

红色表示了滚动平均数。让我们从原始序列中减去这个平均数。注意，从我们采用过去12个月的值开始，滚动平均法还没有对前11个月的值定义：

ts_log_moving_avg_diff = ts_log - moving_avg
ts_log_moving_avg_diff.head(12)

TravelDate
1949-01-01         NaN
1949-02-01         NaN
1949-03-01         NaN
1949-04-01         NaN
1949-05-01         NaN
1949-06-01         NaN
1949-07-01         NaN
1949-08-01         NaN
1949-09-01         NaN
1949-10-01         NaN
1949-11-01         NaN
1949-12-01   -0.065494
Name: Passengers, dtype: float64

注意前11个月是NaN值，现在让我们对这11个月排除后测试稳定性。

ts_log_moving_avg_diff.dropna(inplace=True)
test_stationarity(ts_log_moving_avg_diff)

从上面的测试结果看，已经得到了一个稳定的序列。但是，这个方法有一个缺陷：需要先设定平滑的窗口期，实际在应用的时候像股票这样的数据很难去设定窗口期，所以，我们采取“加权移动平均法”可以对最近的值赋予更高的权重。

加权移动平均法

指数加权移动平均法是很受欢迎的方法，所有的权重被指定给先前的值连同衰减系数。这可以通过pandas实现：

expwighted_avg = ts_log.ewm(halflife=12).mean()
plt.plot(ts_log)
plt.plot(expwighted_avg, color='red')

注意,这里使用了参数“halflife”来定义指数衰减量。这只是一个假设，很大程度上取决于业务领域。其他参数，如跨度和质心也可以用来定义衰减。让我们再来检测下新得到的序列的稳定性：

ts_log_ewma_diff = ts_log - expwighted_avg
test_stationarity(ts_log_ewma_diff)

检测结果比移动平均的效果更好（Test Statistic的值比1%的临界值还小），另外不会出现前面11个月数据遗漏问题。

消除季节性

之前讨论来了简单的趋势减少技术不能在所有情况下使用，特别是在高季节性情况下。让我们谈论一下两种消除趋势和季节性的方法。

• 差分：采用一个特定时间差的差值
• 分解：建立有关趋势和季节性的模型和从模型中删除它们。

差分

处理趋势和季节性的最常见的方法之一就是差分法。在这种方法中,我们采用特定瞬间和它前一个瞬间的不同的观察结果。这主要是在提高平稳性。pandas可以实现一阶差分：

ts_log_diff = ts_log.diff() //ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)

图中可以看出很大程度上减少了趋势。让我们再来检测下：

ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)

我们可以看到平均数和标准差随着时间有小的变化。同时，DF检验统计量小于10% 的临界值，因此该时间序列在90%的置信区间上是稳定的。我们同样可以采取二阶或三阶差分在具体应用中获得更好的结果。这些方法你可以自己尝试。

分解

在这种方法中，趋势和季节性都分别建模，并返回序列的其余部分。

from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log)
 
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
 
plt.subplot(411)
plt.plot(ts_log, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()

趋势和季节性，还有残差值都被分解出来，然后我们计算残差值的稳定性。

ts_log_decompose = residual
ts_log_decompose.dropna(inplace=True)
test_stationarity(ts_log_decompose)

DF测试统计量明显低于1%的临界值，这样时间序列是非常接近稳定。

预测时间序列

我们看到，使用不同的技术都可以是的序列变得稳定，接下里我们以差分处理后的序列搭建模型，因为其相对来说更容易添加噪音及季节性，让其回到预测值。。在执行趋势和季节性预测上，有两种情况：

• 不含依赖值的严格稳定系列。简单的情况下,我们可以建立残差模型作为白噪音（指功率谱密度在整个频域内均匀分布的噪声）。但这是非常罕见的。
• 序列含有明显的依赖值。在这种情况下,我们需要使用一些统计模型像ARIMA（差分自回归移动平均模型）来预测数据。

ARIMA（Auto-Regressive Integrated Moving Averages）模型。平稳时间序列的ARIMA预测的只不过是一个线性方程(如线性回归)。模型有三个主要参数：

• Number of AR (Auto-Regressive) terms (p)： 现在点使用多少个过往数据计算。AR条件仅仅是因变量的滞后。如：如果P等于5，那么预测x（t）将是x（t-1）…x（t-5）。
• Number of MA (Moving Average) terms (q)：使用多少个过往的残余错误值。MA条件是预测方程的滞后预测错误。如：如果q等于5，预测x（t）将是e（t-1）…e(t-5)，e(i)是移动平均叔在第i个瞬间和实际值的差值。
• Number of Differences (d)：为时间序列成为平稳时所做的差分次数。有非季节性的差值，这种情况下我们采用一阶差分。

在这里一个重要的问题是如何确定“p”和“q”的值。我们使用两张图标来确定这些数字。

• 自相关函数（ACF）：这是时间序列和它自身滞后版本之间的相关性的测试。比如在自相关函数可以比较时间的瞬间‘t1’…’t2’以及序列的瞬间‘t1-5’…’t2-5’ (t1-5和t2 是结束点)。
• 部分自相关函数(PACF):这是时间序列和它自身滞后版本之间的相关性测试，但是在预测（已经通过比较干预得到解释）的变量后。如：滞后值为5，它将检查相关性，但是会删除从滞后值1到4得到的结果。

时间序列的自回归函数和部分自回归函数可以在差分后绘制为：

#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
 
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
 
#Plot ACF: 
plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
 
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()

在这个例子中，0刻度线上线的2条虚线为置信区间，用来确定“p”和“q”的值。

• p-部分自相关函数表第一次截断的上层置信区间是滞后值。如果你仔细看，该值是p=0
• q-自相关函数表第一次截断的上层置信区间是滞后值。如果你仔细看，该值是q=1

现在开始搭建ARIMA的三个模型，及计算各自的RSS（这里的RSS是指残差值,而不是实际序列）

自回归（AR）模型

model = ARIMA(ts_log, order=(0, 1, 0))  
results_AR = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))

移动平均数（MA）模型

model = ARIMA(ts_log, order=(0, 1, 1))  
results_MA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2)) 

ARMA组合

在本例中，由于AR的p值为0，所以组合的结果与MA一致：

model = ARIMA(ts_log, order=(0, 1, 1))  
results_ARIMA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))

转化为原数据空间

前面模型采用的数据都是转化后的数据，为了得到最终的结果，需要将数据转化为原数据空间。第一步将预测结果保存为序列。

predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print(predictions_ARIMA_diff.head())

TravelDate
1949-02-01    0.009726
1949-03-01    0.020485
1949-04-01    0.034537
1949-05-01   -0.005924
1949-06-01   -0.006085
dtype: float64

注意这些是从‘1949-02-01’开始而不是第一个月。为什么?这是因为我们将第一个月份取为滞后值，一月前面没有可以减去的元素。一个简单的解决方法是先确定索引的累积和，然后将它添加到基数。累积和的计算方式如下：

predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
print(predictions_ARIMA_diff_cumsum.head())
TravelDate
1949-02-01    0.009726
1949-03-01    0.030211
1949-04-01    0.064748
1949-05-01    0.058824
1949-06-01    0.052739
dtype: float64

将差分转换为对数尺度的方法是这些差值连续地添加到基数。（第一个元素是基数本身）

predictions_ARIMA_log = pd.Series(ts_log.iloc[0], index=ts_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA_log.head()

TravelDate
1949-01-01    4.718499
1949-02-01    4.728225
1949-03-01    4.748710
1949-04-01    4.783247
1949-05-01    4.777323
dtype: float64

最后一步是将指数与原序列比较。

predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))

总结

本周，主要是对一些概念和流程有了更加深入的一些讲解，但是对于如何预测未来一年等方法没有讲到位，后续继续需要学习。