Time series Analysis美赛数学建模数据挖掘

时间序列分析方法

2020-02-24  本文已影响0人  井底蛙蛙呱呱呱

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

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

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

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

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.3 序列分解

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

2、数据预处理

2.1 时间序列平稳性检验

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

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

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

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

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

对于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部分：

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

2.2.1 对数转换（log变换）

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

2.2.2 滑动平均

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

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

2.2.3 差分

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

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


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

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


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

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

通过上面的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）：

因此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）不满足这些属性，则表明它可以进一步改善。模型诊断根据下面的几个方面来判断残差是否符合正态分布：

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

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


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

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

3.3 使用Prophet进行建模

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



上式中，

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

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

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

Prophet对输入数据有要求：

关于 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

上一篇下一篇

猜你喜欢

热点阅读