FlyAI赛题预告：时间序列初学者指南

近期，FlyAI服务竞赛平台将上线时间序列相关赛题，看到此文的童鞋们可以多多关注下；推荐大家可以多在FlyAI竞赛服务平台多参加训练和竞赛，以此来提升自己的能力。FlyAI是为AI开发者提供数据竞赛并支持GPU离线训练的一站式服务平台。每周免费提供项目开源算法样例，支持算法能力变现以及快速的迭代算法模型。

翻译来源：A comprehensive beginner’s guide to create a Time Series Forecast (with Codes in Python and R)

简介

时间序列(从现在开始称为TS)被认为是分析领域中鲜为人知的技能之一（甚至几天前我也不知道）。但正如你们所知，我们的迷你黑客马拉松就是基于此，我了解了解决时间序列问题的基本步骤，在这里分享给你们。这些一定会帮助你在今天的黑客马拉松中得到较好的模型。

在细看文章之前，我强烈推荐阅读 A Complete Tutorial on Time Series Modeling in R，这篇文章像是本文的前传。它基于R语言，主要介绍基础概念，而我利用这些概念使用python代码来解决端到端的问题。R语言有很多时间序列的资料，而python却很少。因为在本文中，我会使用python。

本文将会介绍以下几个方面：
1. 时间序列有什么特别之处？
2. 使用pandas加载和处理时间序列
3. 如何检验时间序列的稳定性？
4. 如何使得一个时间序列变得平稳？
5. 实现时间序列的预测

---

1.时间序列有什么特别之处？

顾名思义，时间序列是时间间隔不变的情况下收集的时间点集合。分析这些集合来确定长期趋势，为了预测未来或进行其他形式的分析。但是什么使TS不同于常规回归问题呢？有两个原因：
1. 时间序列是与时间有关的。因此线性模型的基础假设：观察值是独立的是不适应这个场景的。
2. 伴随着增加和减少的趋势，大多数时间序列会存在季节性趋势，比如，特定时间的特定变化。例如，如果你看到羊毛夹克随时间变化的销量，你一定会发现冬季的销量会很高。

由于时间序列固有许多的特性，所有有各种不同的步骤可以对它进行分析。下文将详细讨论。以使用Python加载一个时间序列对象为开始。我们将使用较为流行的飞机乘客数据集，可以在这里下载。

请记住本文的目的是希望使你熟悉关于时间序列的不同使用方法。本文的例子只是用来方便解释时间序列对象，我重点关注题目的广泛性，不会做非常较精确的预测。

2.使用pandas加载和处理时间序列

Pandas有专门处理时间序列对象的库，特别是可以存储时间信息和允许人们执行快速合作的datatime64[ns]类。我们先导入需要的库：

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

现在，我们可以加载数据集和查看一些最初的行以及列的数据类型。

data = pd.read_csv('AirPassengers.csv')
print(data.head())
print('\n Data Types:')
print(data.dtypes)

 

 

数据包含了指定的月份和该月的游客数量。但是时间序列对象的读取和数据类型的“对象”和“整数”的读取是不一样的。为了将读取的数据作为时间序列，我们必须给read_csv函数传入特殊的参数。

dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month',date_parser=dateparse)
print(data.head())

 

 

我们逐个解释这些参数：

1. parse_dates：（指定该列用于日期解析，与date_parser相对应）这是说明这列含有时间数据信息。正如上面所说的，列的名称为“月份”。
2. index_col：（指定数据中哪一列作为Dataframe的行索引）使用pandas处理时间序列数据背后的关键思想是：时间序列的索引必须是描述时间信息的变量。所以该参数告诉pandas使用“月份”那一列数据作为Dataframe行索引。
3. date_parser：指定一个方法，将输入的字符串转换为时间变量。Pandas默认的数据读取格式是‘YYYY-MM-DD HH:MM:SS’。如需要读取的数据不是默认的格式，就要人工定义。类似于这里定义的dataparse函数的东西可以用于此目的。

现在我们看到数据有作为索引的时间对象和作为列的乘客（#Passengers）。我们可以通过以下指令再次检查索引的数据类型。

data.index

 

 

注意： dtype='datetime[ns]'就确认它是一个时间数据对象。个人而言，我会将列转换为序列对象，这样当我每次使用时间序列的时候，就不需要每次都要提及列名称。当然，这因人而异，如果能令你更好工作，可以使用它作为数据框架。

ts = data['#Passengers'] 
ts.head(10)

 

 

在继续进一步的讨论之前，我会探讨一些关于时间序列数据的索引方法。在序列对象中选择一个特定的值。可以通过一下两种方法来实现：

#1. 将字符串常量作为索引:
ts['1949-01-01']

#2. 导入datatime库且使用datetime方法:
from datetime import datetime
ts[datetime(1949,1,1)]

两者都会返回“112”，可以从先前的输出中证实正确性。假设我们想要直至1949年5月的数据，可以通过两种方法完成：

#1. 明确前后范围:
ts['1949-01-01':'1949-05-01']

#2. 如果其中的一个索引在末端可以使用':':
ts[:'1949-05-01']

两者均会得到以下的输出：

 

 

此处有两个事情需要注意：
1. 不像数字索引，此处的末尾索引被包含的。比如，如果我们使用a[:5]索引一个列表，，它将返回列表中索引值为[0,1,2,3,4]的值。但是在这里索引‘1949-05-01’是包含在结果输出里面的。
2. 索引必须按照时间排序才能工作。如果你随意打乱索引，它将不会工作。

考虑另一个例子，你需要1949年的所有值，这个可以通过下面的指令实现：

ts['1949']

 

月份部分别省略。相似地，如果你想要某个月所有日期的数据，日期部分也可以被省略。
现在，让我们继续分析这个时间序列。

3.如何检验一个时间序列的平稳性？

如果一个时间序列的统计特征如平均数，方差随着时间保持不变，我们就可以认为它是稳定的。为什么时间序列的稳定性这么重要？大部分时间序列模型是在假设它是稳定的前提下建立的。直观地说，我们可以这样认为，如果一个时间序列随着时间产生特定的行为，就有很高的可能性认为它在未来的行为是一样的。同时，与非稳定序列的相比，根据稳定序列得出的理论是更加成熟的， 也是更容易实现。

稳定性的确定标准是非常严格的。但是，如果时间序列随着时间产生恒定的统计特征，根据实际目的我们可以假设该序列是稳定的。如下：

1. 恒定的平均数
2. 恒定的方差
3. 不随时间变化的自协方差

我会跳过一些细节，因为在这篇文章已经解释的非常清楚。接下来，是关于检测平稳性的方法。首先是绘制出数据，进行可视分析。可以通过下面的指令绘制出数据图。

plt.plot(ts)

 

 

非常清晰的看到，随着季节性的变动，飞机乘客的数量总体上是在不断增长的。但是，不是经常都可以获得这样清晰的视觉推论（下文给出相应例子）。所以，更正式的讲，我们可以通过下面的方法检验平稳性。

1. 绘制滚动统计指标：我们可以绘制移动平均数和移动方差，观察它是否随着时间变化。我认为可以通过任何“t”时刻的移动平均数或移动方差，来获得去年的平均数和方差。如：之前的12个月。但是，这更多的是一种视觉技术。
2. DF检验：这是一种检测数据Critical Values平稳性的统计测试。先做一个待验假设(null hypothesis)：时间序列是不平稳的。测试结果由Test Statistic（测试统计量）和一些Critical Values（置信区间的临界值）组成。如果“测试统计量”小于“临界值”，我们可以拒绝待验假设，并认为序列是稳定的。
   这些概念不是凭直觉得出来的。我推荐大家浏览那篇前传文章。如果你对一些理论数据感兴趣，你可以参考Brockwell和Davis的Introduction to Time Series and Forecasting（时间序列和预测的介绍）。这本书有些晦涩，但是如果你能读懂言外之意，你会明白这些概念和正面接触这些统计方法。

回到检验平稳性这件事上，我们将使用滚动统计指标和DF测试结果，因此我已经定义了一个函数，需要时间序列作为输入，然后生成结果。请注意,我已经绘制标准差来代替方差，为了保持单位和平均数相似。

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):

    #Determing rolling statistics
    #原文中的这两行代码中所使用的API将会弃用
    #rolmean = pd.rolling_mean(timeseries, window=12)
    #rolstd = pd.rolling_std(timeseries, window=12)
    rolmean = pd.Series(timeseries).rolling(window=12).mean()
    rolstd = pd.Series(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)

 

 

虽然标准差的变化很小，但随着时间的推移，均值明显增加，这不是一个平稳的序列。此外，测试统计量远远超过临界值。注意，正负值应该被比较，而不是比较绝对值。

接下来，我们将讨论可用于使这个时间序列趋于平稳的技术。

4.如何使一个时间序列趋于平稳？

虽然许多时间序列模型都采用了平稳性假设，但几乎没有一个实际的时间序列是平稳的。统计学家已经找到了使时间序列平稳的方法，我们现在开始讨论。实际上，几乎不可能使一个序列完全平稳，但我们要尽可能地接近它。

先让我们弄明白是什么导致时间序列不稳定。这里有两个主要原因：

1. 趋势-均值随着时间变化而变化。举例：在飞机乘客这个案例中，我们看到总体上，飞机乘客的数量是在不断增长的。
2. 季节性-特定时间帧的变化。举例：在特定的月份购买汽车的人数会有增加的趋势，因为工资上涨或者节假日到来。

模型的根本原理使建模或估计序列中的的趋势和季节性，从序列中删除这些因素，得到一个平稳的序列。然后统计预测技术可以在这个序列上完成。最后一步是通过运用趋势和季节性约束将预测值转换成原来的区间。

注意：我将探讨一系列方法。可能有些对文中情况有用，有些不能。但是我的目的是得到一系列可用方法，而不是仅仅关注目前的问题。

让我们从趋势部分开始。

4.1估计和消除趋势

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

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

 

 

在这个简单的例子中，很容易看到一个向前的数据趋势。但在有噪音的情况下，这不是很直观。因此我们可以使用一些技术来估计或对这个趋势建模,然后将它从序列中删除。这里有很多方法,最常用的有:
1. 聚合——取一段时间的平均值（月/周平均值）
2. 平滑——取滚动平均数
3. 多项式拟合——拟合一个回归模型

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

4.1.1移动平均

在这个方法中，根据时间序列的频率，我们采用K个连续值的平均数。我们可以采用过去一年的平均数，即过去12个月的平均数。关于确定滚动数据，pandas有特定函数。

# 原文代码中使用的API将要被废弃
# moving_avg = pd.rolling_mean(ts_log,12)
moving_avg = pd.Series(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)

 

 

注意前11个月的值是NaN非数字的，现在我们删除这些NaN数据并且绘制出图像来测试平稳性。

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

 

 

这看起来像是个更好的序列。滚动平均值出现轻微的变化，但是没有明显的趋势。同时，检验统计量比5%的临界值小，所以我们在95%的置信区间认为它是平稳序列。

4.1.2指数加权移动平均

然而，上个方法有一个缺陷：要严格定义时段。在这种情况下，我们可以采用年平均数，但是对于复杂的情况的像预测股票价格，是很难得到一个数字的。所以，我们采取“加权移动平均法”可以对最近的值赋予更高的权重。关于指定权重有很多技巧。指数加权移动平均法是很受欢迎的方法，权值被分配给所有之前的值，其中有一个衰减系数。可以在这里找到详细细节。这可以通过pandas实现：

# 原文代码中使用的API将要被废弃
# expwighted_avg = pd.ewma(ts_log, halflife=12)
expweighted_avg = pd.Series(ts_log).ewm(halflife=12).mean()
plt.plot(ts_log)
plt.plot(expwighted_avg, color='red')

 

 

注意,这里使用了参数“halflife（半衰期）”来定义指数衰减量。这只是一个假设，将很大程度上取决于业务领域。其他参数，如span（跨度）和center of mass（质心）也可以用来定义衰减，正如上面分享的链接的探讨。现在让我们从序列中移除这个移动平均，继续检验稳定性。

ts_log_ewma_diff = ts_log - expwighted_avg
test_stationarity(ts_log_ewma_diff)

 

 

这个时间序列的平均值和标准差变化更小。同时，test statistic（检验统计量）小于1%的critical value（临界值），这比之前的情况好。请注意在这种情况下就不会有遗漏值因为所有的值在一开始就被赋予了权重。所以即使没有之前的值它也能工作。

4.2消除趋势和季节性

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

1. 差分——取一个特定时间间隔的差值
2. 分解——建立有关趋势和季节性的模型然后从模型中删除它们。

4.2.1差分

处理趋势和季节性的最常见的方法之一就是差分法。在这种方法中,我们取特定时刻和它前一个时刻的差值。这个方法在大多数情况下都可以很好地提高平稳性。pandas可以实现一阶差分：

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%的置信区间内认为时间序列是平稳的。我们还可以取二阶或三阶差分，这在某些应用中可能会得到更好的结果。你们可以继续尝试。

4.2.2分解

在这种方法中，趋势和季节性被分开建模，这个序列剩下的部分被返回。我将会跳过统计指标，直接给出结果。

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%的临界值，这样时间序列是非常接近平稳的。你也可以尝试更高级的分解技术从而产生更好的结果。同时，你应该注意到，在这种情况下将剩余值（残差）转换为原始值对未来数据不是很直观。

5.时间序列预测

我们看到不同的技术，它们可以使时间序列稳定。让我们建立差分后的时间序列模型，因为差分是很受欢迎的技术，也相对更容易添加噪音和对预测的残差进行季节性回退。在完成了趋势和季节性估计和消除技术后,有两种情况:

1. 不含依赖值的严格平稳序列。简单的情况下,我们可以将残差作为白噪音。但这是非常罕见的。
2. 序列含有明显的依赖值。在这种情况下,我们需要使用一些统计模型像ARIMA来预测数据。

让我给你简要介绍一下ARIMA，我不会介绍技术细节,但如果你希望更有效地应用它们，你应该理解这些概念的细节。ARIMA代表差分自回归移动平均。平稳时间序列的ARIMA预测的只不过是一个线性方程(如线性回归)。预测依赖于ARIMA模型参数(p,d,q)：

1. 自回归项的数目(p)：AR项仅仅是因变量的时滞。例如，如果p等于5，那么预测x(t)将是x(t-1)…x(t-5)。
2. 移动平均的数目(q)：MA项是预测误差的时滞。例如，如果q等于5，那么预测x(t)将是e(t-1)…e(t-5)，e(i)是第i个时刻的移动平均和真实值的差值。
3. 差分次数(d)：这里指的是非季节性的差分次数，即这种情况下我们采用一阶差分。因为无论我们传递差分后的变量且d=0，还是传递原始变量且d=1，得到的结果是一样的。

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

1. 自相关函数(ACF)：这是时间序列和它自身滞后版本之间的相关性的测量。比如时滞为5，自相关函数会比较‘t1’到‘t2’时刻的序列和‘t1-5’到‘t2-5’时刻的序列。
2. 偏自相关函数(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旁边的两条虚线之间代表置信区间。这些可以用来确定“p”和“q”的值:
1. p——偏自相关函数图中，第一次截断上层置信区间的所在位置是滞后值。如果你仔细看，该值是p=2。
2. q——自相关函数图中，第一次截断上层置信区间的所在位置是滞后值。如果你仔细看，该值是q=2。

现在，考虑个体以及组合的效果，建立3个不同的ARIMA模型。我也会打印出各自的RSS。请注意,这里的RSS是指残差值,而不是实际序列。

首先，我们需要加载ARIMA模型。

from statsmodels.tsa.arima_model import ARIMA

p，d，q值可以使用ARIMA的order参数来指定，即一个元组(p,d,q)。建立三种情况下的模型：

5.1 AR模型

model = ARIMA(ts_log, order=(2, 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))

 

 

5.2 MA模型

model = ARIMA(ts_log, order=(0, 1, 2))  
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))

 

 

5.3 组合模型

model = ARIMA(ts_log, order=(2, 1, 2))  
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))

 

 

在这里我们可以看到,AR模型和MA模型几乎有相同的RSS，但组合模型效果明显更好。现在,我们只剩下最后一步,即把这些值回退到原始区间。

5.4 回退至原始区间

既然组合模型获得更好的结果，让我们将它倒回原始值，看看它的预测表现如何。第一步是作为一个独立的序列存储预测结果，然后观察它。

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

 

 

注意，这些是从‘1949-02-01’开始，而不是第一个月。为什么?这是因为我们将第一个月份取为滞后值，一月前面没有可以减去的元素。将差分后的区间转换为取对数后的区间的方法是这些差值连续地添加到基本值。（这句话就是想要回退差分）一个简单的方法就是首先确定现在区间的累计总和,然后将其添加到基本值。累计总和可以由下面得到:

predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
print(predictions_ARIMA_diff_cumsum.head())

 

你可以对输出结果进行回算，检查这些是否正确的。接下来我们将它们添加到基本值。为此我们将使用所有的值创建一个序列作为基本值，并添加差值。我们这样做：

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

 

第一个元素是基本值本身，从基本值开始值累计添加。最后一步是做指数操作，然后与原序列比较。

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

 

 

最后我们获得一个原始区间的预测结果。虽然不是一个很好的预测。但是你获得了思路对吗?现在，我把它留给你去进一步改进，实现一个更好的方案。

最后注意

在本文中，我试图给你们提供一个标准方法去解决时间序列问题。这个来得正好，因为今天是我们的小型编程马拉松，就是挑战你们是否可以解决类似的问题。我们广泛的讨论了平稳性的概念，怎样使时间序列趋近平稳和最终的预测残差。这是一个漫长的过程，我跳过了一些统计细节，我鼓励大家使用文中给出的链接作为参考材料。如果您不想复制粘贴，可以从我的GitHub下载带有所有代码的iPython notebook。

（本翻译中对应的代码相较于原文有部分改动，比如某些废弃api的修改，以及python2写法改为python3写法）

 

————————————————

更多精彩内容请访问FlyAI-AI竞赛服务平台；为AI开发者提供数据竞赛并支持GPU离线训练的一站式服务平台；每周免费提供项目开源算法样例，支持算法能力变现以及快速的迭代算法模型。

挑战者，都在FlyAI！！！