时间序列分析方法

作者: 井底蛙蛙呱呱呱 | 来源:发表于2020-02-24 00:18 被阅读0次

时间序列是指一组在连续时间上测得的数据，其在数学上的定义是一组向量x(t), t=0,1,2,3,...，其中t表示数据所在的时间点，x(t)是一组按时间顺序（测得）排列的随机变量。包含单个变量的时间序列称为单变量时间序列，而包含多个变量的时间序列则称为多变量。

时间序列在很多方面多有涉及到，如天气预报，每天每个小时的气温，股票走势等等，在商业方面有诸多应用，如：

销售预测，零售产品销量预测;
需求预测，用于定价，库存和劳动力管理；
交通预测，运输和路线优化，道路设施设计；
收入预测，预算和目标设定。

下面我们将通过一个航班数据来说明如何使用已有的工具来进行时间序列数据预测。常用来处理时间序列的包有三个：

最常见的AR、ARMA、ARIMA建模的工具包pmdarima，这个包能非常方便的对平稳序列进行建模，对于但是其无法处理包含 seasonal 成分的序列；
statmodels中的tsa模块，该模块中同样也实现了AR、ARMA、ARIMA方法，除此之外，该模块还实现了SARIMA方法，该方法在ARIMA的基础上又添加了对seasonal 部分建模的方法，因此适用度更高；
fbprophet 包，该包实现的算法与上面的不同，其对序列没有平稳性要求，对于任何序列都可以直接来建模，不过其对训练数据有一定的要求，一般需要几个月的训练数据才能够达到比较好的效果。

对于基于AR、MA的方法一般需要数据预处理，因此本文分为三部分：

1、数据探索分析，包括数据查看数据的一些基本统计量（均值方差等），异常检测（包括缺失值、离群点等），数据分布（是否正态分布）和数据的发展趋势（上升、下降还是基本持平）等；
2、数据预处理。这一步主要是将数据处理成更适合模型建模的数据分布，如将非正态分布转换为正态分布，对于AR、MA相关建模方法，一般会先将非平稳序列转换为平稳序列；
3、建模分析，使用时间序列建模算法对数据进行建模分析。

1、数据探索性分析

1.1 数据基本统计量查看

import warnings
warnings.filterwarnings('ignore')

import itertools
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

import pandas as pd
pd.set_option('display.expand_frame_repr', False)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)

import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.ar_model import AR
from statsmodels.tsa.arima_model import ARMA, ARIMA
from pmdarima.arima import auto_arima
from statsmodels.tsa.statespace.sarimax import SARIMAX
from fbprophet import Prophet

from math import sqrt

import matplotlib
matplotlib.rcParams['axes.labelsize'] = 14
matplotlib.rcParams['xtick.labelsize'] = 12
matplotlib.rcParams['ytick.labelsize'] = 12
matplotlib.rcParams['text.color'] = 'k'

import seaborn as sns

from random import random

from sklearn.metrics import mean_squared_error, r2_score, \
    mean_absolute_error, median_absolute_error, mean_squared_log_error

1.2 数据可视化分析

通过简单的初步处理以及可视化可以帮助我们有效快速的了解数据的分布（以及时间序列的趋势）。

观察数据的频率直方图以及密度分布图以洞察数据结构，从下图可以看出：

1. 该数据不是一个完美的正态分布；
1. 该数据分布左移，因此在后续分析前进行转换可能会有用。
  
  观察数据的箱线图：
每一年的中间值都呈上升趋势证明了这个数据随着时间的推移呈现稳定的上升趋势；
随着时间的推移中间数据量（25%-50%区间）稳定上升；
不同月份的数值呈现一个先上升后下降的趋势，结合上面图的每年数值总体上升的整体趋势，可以猜测一年之中不同月份的数据可能存在一个规律，因此考虑季节性（或者说部分周期性）模型会比较好。

1.3 序列分解

使用statsmodels对该时间序列进行分解，以了解该时间序列数据的各个部分，每个部分都代表着一种模式类别。借用statsmodels序列分解我们可以看到数据的主要趋势成分、季节成分和残差成分，这与我们上面的推测相符合。

2、数据预处理

2.1 时间序列平稳性检验

如果一个时间序列的均值和方差随着时间变化保持稳定，则可以说这个时间序列是稳定的。

大多数时间序列模型都是在平稳序列的前提下进行建模的。造成这种情况的主要原因是序列可以有许多种（复杂的）非平稳的方式，而平稳性只有一种，更加的易于分析，易于建模。

在直觉上，如果一段时间序列在某一段时间序列内具有特定的行为，那么将来很可能具有相同的行为。譬如已连续观察一个星期都是六点出太阳，那么可以推测明天也是六点出太阳，误差非常小。

而且，与非平稳序列相比，平稳序列相关的理论更加成熟且易于实现。

一般可以通过以下几种方式来检验序列的平稳性：

1. ACF与PACF曲线。如果时间序列是稳定的，则ACF/PACF曲线中某一观测数据与其滞后数据点呈现很小的相关性并且这个相关性会迅速下降消失；
1. 滑动均值和方差的变化趋势。通过绘制滑动均值和方差的的曲线我们可以看出时间序列的均值和方差是否随着时间序列变化；
1. adf检验 (Augmented Dickey-Fuller Test)。这是用于检验数据平稳性的统计检验方法之一。在adf检验中，零假设为时间序列是非平稳的。测试结果包括一些统计量、置信度和临界值。如果检验统计量小于关键值，则我们可以拒绝原假设并认为该序列是平稳的。

对于ACF和PCF的一些解释：
ACF(Autocorrelation Function):
我们可以假设每个变量的分布都符合高斯分布（钟形曲线）。如果是这种情况，我们可以使用Pearson的相关系数来总结变量之间的相关性。皮尔逊相关系数是介于-1和1之间的数字，分别描述了负相关或正相关，零值表示没有相关性。

我们可以计算时间序列观测值与先前时间步长的观测值的相关性，称为滞后。因为时间序列观测值的相关性是使用同一序列某一时刻值域先前时间的值计算的，所以这称为序列相关性或自相关。时间序列按时滞的自相关图称为自相关函数，或简称ACF。该图有时称为相关图或自相关图。

PACF(Partial Autocorrelation Function):
PACF相当于在计算X(t)和X(t-h)的相关性的时候，挖空在(t-h,t)上所有数据点对X(t)的影响，去反映X(t-h)和X(t)之间真正的相关性（直接相关性）。

某一观测值与先前时间观测值之间的相关性包括直接相关和间接相关。这些间接性是观测值相关性的线性函数。偏自相关函数试图消除的就是这些间接相关。

个人理解，这里的直接相关可以看作是真实数据值，而间接相关则可以看作是随机误差值。因此从这点来说AR模型看PACF(直接相关性，直接对数据进行建模)，而MA模型看ACF（对误差进行建模）。

参考:
A Gentle Introduction to Autocorrelation and Partial Autocorrelation
计量经济学中，ACF和PACF函数有什么区别？

2.1.1 通过ACF和PACF曲线来查看时间序列的平稳性

如果时间序列是平稳性的，那么在ACF/PACF中观测点数据与之前数据点的相关性会急剧下降。

下图中的圆锥形阴影是置信区间，区间外的数据点说明其与观测数据本身具有强烈的相关性，这种相关性并非来自于统计波动。

PACF在计算X(t)和X(t-h)的相关性的时候，挖空在（t-h，t）上所有数据点对X(t)的影响，反应的是X(t)和X(t-h)之间真实的相关性（直接相关性）。

从下图可以看出，数据点的相关性并没有急剧下降，因此该序列是非平稳的。

2.1.2 通过滑动均值/方差来检验序列的平稳性

如果序列是平稳的，那么其滑动均值/方差会随着时间的变化保持稳定。

但是从下图我们可以看到，随着时间的推移，均值呈现明显的上升趋势，而方差也呈现出波动式上升的趋势，因此该序列是非平稳的。

2.1.3 adf检验序列的稳定性

一般来讲p值小于0.05我们便认为其是显著性的，可以拒绝零假设。但是这里的p值为0.99明显是非显著性的，因此接受零假设，该序列是非平稳的。

下面的统计量是adf的结果，上面的图是滑动均值/方差

2.2 数据转换——使非平稳序列转换为平稳序列

从上面的平稳性检验我们可以知道该时间序列为非平稳序列。此外，通过上面1.3部分的序列分解我们也可以看到，该序列可分解为3部分：

长期趋势（Trend）因素。一般来说该部分的均值会随着时间变化；
季节性（seasonal）因素。一般来说该部分变化是有规律的，呈周期性的变化。如大家每年都会在冬天买棉袄夏天买短袖T恤而不会反过来。

我们可以使用数据转换来对那些较大的数据施加更大的惩罚，如取对数、开平方根、立方根、差分等，以达到序列平稳的目的。

2.2.1 对数转换（log变换）

看来取对数在这里貌似没什么用- -

2.2.2 滑动平均

我们对上面取log后的数据进行滑动平均。

原理

滑动平均后数据失去了其原来的特点（波动式上升），这样损失的信息过多，肯定是无法作为后续模型的输入的。

2.2.3 差分

差分是常用的将非平稳序列转换平稳序列的方法。ARIMA中的 'I' 便是指的差分，因此ARIMA是可以对非平稳序列进行处理的，其相当于先将非平稳序列通过差分转换为平稳序列再来使用ARMA进行建模。

一般差分是用某时刻数值减去上一时刻数值来得到新序列。但这里有一点区别，我们是使用当前时刻数值来减去其对应时刻的滑动均值。

由于滑动均值其前11个位置没有值，因此这里差分后序列的前11个也没有值

我们来看看刚刚差分的结果怎么样。

蓝线即为取log+差分的结果
让我们在这里稍微停一下，我们刚刚干了什么？p值变为0.02，这意味着我们的时间序列竟然转换成平稳序列了（尽管其滑动均值看起来略有波动）。

让我们稍微总结下我们刚刚的步骤：

1. 对原始数据取对数；
1. 求取各个位置的滑动平均值；
1. 各个位置取对数后的值减去该位置上对应的滑动均值，得到转换后的平稳序列。

通过上面的3步我们成功的将一个非平稳序列转换成了一个平稳序列。上面使用的是最简单的滑动均值，下面我们试试指数滑动平均怎么样。

2.2.4 指数滑动平均

上面是最常用的指数滑动平均的定义，但是pandas实现的指数滑动平均好像与这个有一点区别，详细区别还得去查pandas文档。

指数滑动均值的效果看起来也很差。我们使用差分+指数滑动平均再来试试吧。

2.3 分解处理后的序列

在上面我们通过取log+(指数)滑动平均+差分已经成功将非平稳序列转换为了平稳序列。

下面我们看看，转换后的平稳序列的各个成分是什么样的。不过这里我们使用的是最简单的差分，当前时刻的值等于原始序列当前时刻的值减去原始序列中上一时刻的值，即： x'(t) = x(t) - x(t-1)。

看起来挺不错，是个平稳序列的样子。不过，还是检验一下吧。

经过上面的转换后序列已经基本可以认为是平稳序列了。然我们来看看其各个组分现在是什么样子。

可以看到，趋势(Trend)部分已基本被去除，但是季节性(seasonal)部分还是很明显，而ARIMA是无法对含有seasonal的序列进行建模分析的。

3、时间序列建模预测

在一开始我们提到了3个包均可以对时间序列进行建模。

为了简便，这里pmdarima和statsmodels.tsa直接使用最好的建模方法即SARIMA，该方法在ARIMA的基础上添加了额外功能，可以拟合seasonal部分以及额外添加的数据。

3.1 使用`statsmodels.tsa`中的方法来进行建模

在使用ARIMA（Autoregressive Integrated Moving Average）模型前，我们先简单了解下这个模型。这个模型其实可以包括三部分，分别对应着三个参数（p, d, q）：

AR（Autoregressive，p）,自回归模型就是利用滞后的时间来预测当前时间点的数据，如使用x(t-1), x(t-2)和x(t-3)来拟合预测x(t)。这里使用的滞后变量的数量就是p； AR模型
MA（Moving Averages，q），滑动平均模型使用滞后序列的白噪声来拟合当前数据，同样的滞后数据数量即是q； MA模型
Difference（即差分，d），上面讲了通过差分可以将非平稳序列转变为平稳的，这个差分的作用与其相同，因为AR和MA无法对非平稳序列进行拟合。差分

因此ARIMA模型就是将AR和MA模型结合起来然后加上差分，克服了不能处理非平稳序列的问题。但是，需要注意的是，其仍然无法对seasonal进行拟合。

ARMA模型
更详细的介绍可以参考数据分析技术：时间序列分析的AR/MA/ARMA/ARIMA模型体系。

下面开始使用ARIMA来拟合数据。

（1）先分训练集和验证集。需要注意的是这里使用的原始数据来进行建模而非转换后的数据。

（2）ARIMA一阶差分建模并预测

查看一下预测的未还原差分的结果：

（3）对差分结果进行还原

可视化一下看看结果怎么样。

结果很差呀 - -
（4）使用模型评估指标对模型进行评估试试

（5）使用SARIMA来建模试试。上面我们使用ARIMA虽然使用了差分能将非平稳序列转换为平稳序列，但是ARIMA是无法处理seasonal部分的，而SARIMA则可以，下面我们使用SARIMA试试。

Seasonal Autoregressive Integrated Moving-Average (SARIMA)

SARIMA是ARIMA的扩展，他明确支持具有季节性成分的单变量时间序列数据。该实现称为SARIMAX而不是SARIMA是因为该实现还支持可选的外部变量，通过exog参数指定，外部变量的实例：人口，假期，航空公司数量，重大事件。

他添加了3个新的超参数，以指定序列的季节性分量的自回归（AR），差分（I）和移动平均值（MA），以及季节性周期的附加参数m。

Trend因素有三个参数控制，这几个参数与ARIMA中是一样的：

p，trend AR阶数（order）；

d，trend差分阶数；

q，trend MA阶数。

Seasonal因素有四个参数控制：

P，季节性AR阶数；

D，季节性差分阶数；

Q，季节性MA阶数；

m，单个季节周期的时间步数。例如，月度数据的S为12表示每年的季节周期。

SARIMA表示法：SARIMA（p，d，q）（P，D，Q，m）。

auto_arima documentation for selecting best model

先手动选择几组参数，然后参数搜索找到最佳值。需要注意的是，为了避免过拟合，这里的阶数一般不太建议取太大。

可视化看看结果怎么样吧。

效果相当不错呢。再来看看评估指标怎么样。

各个指标都还不错呢，至少比上面的ARIMA是进不了不少，可喜可贺。看来模型对seasonal部分的拟合还是非常重要的啊。

（6）最后，我们还能对拟合好的模型进行诊断看看结果怎么样。

我们主要关心的是确保模型的残差（residual）部分互不相关，并且呈零均值正态分布。若季节性ARIMA（SARIMA）不满足这些属性，则表明它可以进一步改善。模型诊断根据下面的几个方面来判断残差是否符合正态分布：

下图的右上角图中，红色的KDE曲线与N(0, 1)曲线形状大致相同，说明残差是基本符合正态分布的；
左下角的QQ图表明，蓝色（残差）序列的分布遵循正态分布样本的线性趋势。这也说明残差符合正态分布；
左上图随着时间的推移，没有显示出任何的季节性变化，说明其是白噪声；
右下角的自相关图表明观测序列与其之前的相关性非常小。

3.2 使用 `pmdarima` 中的包来进行建模

同样的，为了方便，我们这里使用pmdarima 中一个可以自动搜索最佳参数的方法auto_arima来进行建模。

可以看到，这里是对seasonal部分进行建模了的。看看结果怎么样吧。

结果很不错啊，虽然比上面的最好的要差些，但是比起ARIMA还是要好上很多的。

3.3 使用`Prophet`进行建模

一般来说，在实际生活和生产环节中，除了季节项，趋势项，剩余项之外，通常还有节假日的效应。所以，在prophet算法里面，作者同时考虑了以上四项，即：

上式中，

g(t) 表示趋势项，它表示时间序列在非周期上面的变化趋势；
s(t) 表示周期项，或者称为季节项，一般来说是以周或者年为单位；
h(t) 表示节假日项，表示在当天是否存在节假日；
ξ_t 表示误差项或者称为剩余项。

更多详细Prophet算法内容可以参考 Facebook 时间序列预测算法 Prophet 的研究。

Prophet算法就是通过拟合这几项，然后把它们累加起来得到时间序列的预测值。

Prophet提供了直观且易于调整的参数：

趋势参数（Trend Params）：
- growth：指定线性或逻辑趋势；
- changepoints：包括潜在变更点的日期列表（如未指定，则为自动）；
- n_chagepoints：如果未提供更改点（changepoint），可以提供自动寻找changgepoint的数量；
- chagepoint_prior_scale：用于控制自动选择changepoint的灵活性；
季节性和假期参数（Seasonality and Holiday Params）：
- yearly_seasonality: 用来拟合年度季节性；
- weekly_seasonality: 用来拟合每周的季节性；
- daily_seasonality: 用来拟合每日的季节性；
- holidays: 形式为dataframe，包含假日名和日期；
- seasonality_prior_scale: 用来调整模型季节性强度的参数；
- holiday_prior_scale: 用于调整模型假期强度的参数。

Prophet对输入数据有要求：

y, target, numeric;
ds, datetime

关于 Prophet 的使用例子可以参考 Prophet example notebooks

下面使用 Prophet 来进行处理数据。

蓝色表示预测结果
查看各个部分的变化情况。由于我们上面仅选择了yearly_seasonality=True，因此下面的seasonal中只有year。

最后评估一下模型的预测结果怎么样。

这结果可以说是相当好了，而且是没有附加节假日信息，没有调参的结果。

参考：
Facebook 时间序列预测算法 Prophet 的研究
 Prophet example notebooks
auto_arima documentation for selecting best model
数据分析技术：时间序列分析的AR/MA/ARMA/ARIMA模型体系
 https://github.com/advaitsave/Introduction-to-Time-Series-forecasting-Python
时间序列分析
 My First Time Series Comp (Added Prophet)
Prophet官方文档：https://facebookincubator.github.io

时间序列分析方法