机器学习（11）——时间序列分析_时间序列 分阶段 机器学习

目录

1 时间序列数据的相关检验

1.1 白噪声检验

1.2 平稳性检验

1.3 自相关分析和偏自相关分析

2 移动平均算法

2.1 简单移动平均法

2.2 简单指数平滑法

2.3 霍尔特线性趋势法

2.4 Holt-Winters 季节性预测模型

3 ARIMA 模型

3.1 AR 模型

3.2 ARMA 模型

3.3 ARIMA 模型

3.3.1 训练ARIMA模型

3.3.2自动搜索ARIMA模型的参数

4 SARIMA 模型

5 Prophet 模型预测时间序列

5.1 数据准备

5.2 模型建立与数据预测

6 多元时间序列ARIMAX模型

6.1 数据准备与可视化

6.2 ARIMAX模型建立与预测

7 时序数据的异常值检测

7.1 数据准备与可视化

7.2 时序数据异常值检测

---

时间序列数据是常见的数据类型之一,时间序列分析基于随机过程理论和数理统计学方法,研究时间序列数据所遵从的统计规律,常用于系统描述、系统分析、预测未来等。
时间序列数据主要是根据时间先后,对同样的对象按照等时间间隔收集的数据,比如每日的平均气温、每天的销售额、每月的降水量等。虽然有些序列所描述的内容取值是连续的,比如气温的变化可能是连续的,但是由于观察的时间段并不是连续的,所以可以认为是离散的时间序列数据。
一般地,对任何变量做定期记录就能构成一个时间序列。根据所研究序列数量的不同,可以将时间序列数据分为一元时间序列数据和多元时间序列数据。
时间序列的变化可能受一个或多个因素的影响,导致它在不同时间的取值有差异,这些影响因素分别是长期趋势、季节变动、循环波动(周期波动)和不规则波动(随机波动)。时间序列分析主要有确定性变化分析和随机性变化分析。确定性变化分析包括趋势变化分析、周期变化分析、循环变化分析。随机性变化分析主要有AR、MA、ARMA、ARIMA模型等。

首先导入本章会使用到的库和模块,程序如下:

## 图像显示中文的问题
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False
 
import seaborn as sns 
sns.set(font= "Kaiti",style="ticks",font_scale=1.4)
 
## 导入会使用到的相关库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import *
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.tsa.api import SimpleExpSmoothing,Holt,ExponentialSmoothing,AR,ARIMA,ARMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
 
import pmdarima as pm
from sklearn.metrics import mean_absolute_error
 
import pyflux as pf
from fbprophet import Prophet
 
 
## 忽略提醒
import warnings
warnings.filterwarnings("ignore")

时间序列模型的预测主要可以通过statsmodels 库的tsa模块来完成。针对时间序列数据,常用的分析流程如下:
(1)根据时间序列的散点图、自相关函数和偏自相关函数图等识别序列是否是非随机序列,如果是非随机序列,则观察其平稳性。
(2)对非平稳的时间序列数据采用差分进行平稳化处理,直到处理后序列是平稳的非随机序列。
(3)根据所识别出来的特征建立相应的时间序列模型。
(4)参数估计,检验是否具有统计意义。
(5)假设检验,判断模型的残差序列是否为白噪声序列。
(6)利用已通过检验的模型进行预测。

1 时间序列数据的相关检验

对于时间序列数据,最重要的检验就是时间序列数据是否为白噪声数据、时间序列数据是否平稳,以及对时间序列数据的自相关系数和偏自相关系数进行分析。如果时间序列数据是白噪声数据, 说明其没有任何有用的信息。针对时间序列数据的很多分析方法,都要求所研究的时间序列数据是平稳的,所以判断时间序列数据是否平稳,以及如何将非平稳的时间序列数据转化为平稳序列数据,对时间序列数据的建模研究是非常重要的。

1.1 白噪声检验

本节将会利用两个时间序列数据进行相关的检验分析,首先读取数据并使用折线图将两组时间序列进行可视化,运行下面的程序后,结果如图6-1所示。

文件提取：

链接：https://pan.baidu.com/s/1EJEOq_FuJpf2y2Rw-lcC5g 
提取码：whj6

## 读取时间序列数据,该数据包含:X1为飞机乘客数据,X2为一组随机数据
df = pd.read_csv("E:/PYTHON/timeserise.csv")
## 查看数据的变化趋势
df.plot(kind = "line",figsize = (10,6))
plt.grid()
plt.title("时序数据")
plt.show()

运行结果如下：

图6-1 序列的波动情况

如果一个序列是白噪声(即独立同分布的随机数据),那么就无须再对其建立时间序列模型来预测,因为预测随机数是无意义的。因此在建立时间序列分析之前,需要先对其进行白噪声检验。
常用的白噪声检验方法是Ljung-Box检验(简称LB检验),其原假设和备择假设分别为HO:延迟期数小于或等于m期的序列之间相互独立(序列是白噪声); H1:延迟期数小于或等于m期的序列之间有相关性(序列不是白噪声)。 Ljung-Box 检验可以使用 sm.stats.diagnostic.acorr_ljungbox()函数,对两个序列进行白噪声检验,程序如下: 

## 白噪声检验Ljung-Box检验
## 该检验用来检查序列是否为随机序列，如果是随机序列，那它们的值之间没有任何关系
## 使用LB检验来检验序列是否为白噪声，原假设为在延迟期数内序列之间相互独立。
lags = [4,8,16,32]
LB = sm.stats.diagnostic.acorr_ljungbox(df["X1"],lags = lags,return_df = True)
print("序列X1的检验结果:\n",LB)
LB = sm.stats.diagnostic.acorr_ljungbox(df["X2"],lags = lags,return_df = True)
print("序列X2的检验结果:\n",LB)
 
## 如果P值小于0.05，说明序列之间不独立，不是白噪声

运行结果如下：

序列X1的检验结果:
         lb_stat      lb_pvalue
4    427.738684   2.817731e-91
8    709.484498  6.496271e-148
16  1289.037076  1.137910e-264
32  1792.523003   0.000000e+00
序列X2的检验结果:
       lb_stat  lb_pvalue
4    1.822771   0.768314
8    8.452830   0.390531
16  15.508599   0.487750
32  28.717743   0.633459

从上面的结果中可以看出,在延迟阶数为[4,8,16,32]的情况下,序列X1的LB检验P值均小于0.05,说明可以拒绝序列为白噪声的原假设,认为该数据不是随机数据，即该数据不是随机的，是有规律可循的,有分析价值。而序列X2的LB检验P值均大于0.05,说明该序列为白噪声,没有分析价值。

1.2 平稳性检验

时间序列是否是平稳的,对选择预测的数学模型非常关键。如果一组时间序列数据是平稳的, 就可以直接使用自回归移动平均模型(ARMA)进行预测,如果数据是不平稳的,就需要尝试建立差分移动自回归平均模型(ARIMA)等进行预测。
判断序列是否平稳有两种检验方法:一种是根据时序图和自相关图显示的特征做出判断；另一种是构造检验统计量进行假设检验,如单位根检验。第一种判断方法比较主观,第二种方法则是客观的判断方法。
常用的单位根检验方法是ADF检验,它能够检验时间序列中单位根的存在性,其检验的原假设和备择假设分别为HO:序列是非平稳的(序列有单位根); H1:序列是平稳的(序列没有单位根)。
Python中sm.tsa模块的adfuller()函数可以进行单位根检验,针对序列X1和X2可以使用下面的程序进行单位根检验。

## 序列的单位根检验，即检验序列的平稳性
dftest = adfuller(df["X2"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X2单位根检验结果:\n",dfoutput)
 
dftest = adfuller(df["X1"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X1单位根检验结果:\n",dfoutput)
 
## 对X1进行一阶差分后的序列进行检验
X1diff = df["X1"].diff().dropna()
dftest = adfuller(X1diff,autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X1一阶差分单位根检验结果:\n",dfoutput)
 
## 一阶差分后 P值大于0.05, 小于0.1，可以认为其是平稳的

运行结果如下：

X2单位根检验结果:
 adf                           -1.124298e+01
p-value                        1.788000e-20
usedlag                        0.000000e+00
Number of Observations Used    1.430000e+02
dtype: float64
X1单位根检验结果:
 adf                              0.815369
p-value                          0.991880
usedlag                         13.000000
Number of Observations Used    130.000000
dtype: float64
X1一阶差分单位根检验结果:
 adf                             -2.829267
p-value                          0.054213
usedlag                         12.000000
Number of Observations Used    130.000000
dtype: float64

从上面的单位根检验的输出结果中可以发现,序列X2的检验P值小于0.05,说明X2是一个平稳时间序列(注意该序列属于白噪声,白噪声序列是平稳序列)。针对序列X1的单位根检验, 可发现其P值远大于0.05,说明其实不平稳,而针对其一阶差分后的结果可以发现,一阶差分后P值大于0.05,但是小于0.1,可以认为其是平稳序列。
针对数据的平稳性检验,还可以使用KPSS检验,其原假设为检测的序列是平稳的。该检验可以使用kpss()函数来完成,使用该函数对序列进行检验的程序如下：

## KPSS检验的原假设为：序列x是平稳的。
 
## 对序列X2使用KPSS检验平稳性
dfkpss = kpss(df["X2"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X2 KPSS检验结果:\n",dfoutput)
## 接受序列平稳的原假设
 
 
## 对序列X1使用KPSS检验平稳性
dfkpss = kpss(df["X1"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1 KPSS检验结果:\n",dfoutput)
## 拒绝序列平稳的原假设
 
## 对序列X1使用KPSS检验平稳性
dfkpss = kpss(X1diff)
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1一阶差分KPSS检验结果:\n",dfoutput)
## 接受序列平稳的原假设

运行结果如下：

X2 KPSS检验结果:
 kpss_stat     0.087559
 p-value      0.100000
 usedlag     14.000000
dtype: float64
X1 KPSS检验结果:
 kpss_stat     1.052175
 p-value      0.010000
 usedlag     14.000000
dtype: float64
X1一阶差分KPSS检验结果:
 kpss_stat     0.05301
 p-value      0.10000
 usedlag     14.00000
dtype: float64

从输出的检验结果中可以知道,序列X2是平稳序列,序列X1是不平稳序列,X1一阶差分后的序列是平稳序列。
针对时间序列ARIMA(p,d,g)模型,参数d可以通过差分次数来确定,也可以利用pm.arima模块的ndiffs()函数进行相应的检验来确定。如果对序列建立ARIMA模型可以使用下面的程序确定参数d的取值:

## 检验ARIMA模型的参数d
X1d = pm.arima.ndiffs(df["X1"], alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法对序列X1的参数d取值进行预测,d = ",X1d)
 
X1diffd = pm.arima.ndiffs(X1diff, alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法对序列X1一阶差分后的参数d取值进行预测,d = ",X1diffd)
 
X2d = pm.arima.ndiffs(df["X2"], alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法对序列X2的参数d取值进行预测,d = ",X2d)

运行结果如下：

使用KPSS方法对序列X1的参数d取值进行预测,d =  1
使用KPSS方法对序列X1一阶差分后的参数d取值进行预测,d =  0
使用KPSS方法对序列X1的参数d取值进行预测,d =  0

从输出的结果中可以发现,针对平稳序列获得的参数d取值为0,而针对不平稳的时间序列X1其参数d的预测结果为1。

针对时间序列SARIMA模型,还有一个季节周期平稳性参数D需要确定,同时也可以利用pm.arima模块中的nsdiffs()函数进行相应的检验来确定,使用该函数的程序示例如下:

## 检验SARIMA模型的参数季节阶数D
X1d = pm.arima.nsdiffs(df["X1"], 12, max_D=2)
print("对序列X1的季节阶数D取值进行预测,D = ",X1d)
 
X1diffd = pm.arima.nsdiffs(X1diff, 12, max_D=2)
print("序列X1一阶差分后的季节阶数D取值进行预测,D = ",X1diffd)

运行结果如下：

对序列X1的季节阶数D取值进行预测,D =  1
序列X1一阶差分后的季节阶数D取值进行预测,D =  1

从程序的输出结果中可以发现,序列X1和序列X1一阶差分后的序列,检验结果都为D=1。

1.3 自相关分析和偏自相关分析

自相关分析和偏自相关分析,是用来确定ARMA(p,q)模型中两个参数p和g的一种方法,在 确定序列为平稳的非白噪声序列后,可以通过序列的自相关系数和偏自相关系数取值的大小来分析 序列的截尾情况。
对于一个时间序列,如果样本的自相关系数ACF不等于0,直到滞后期s=q,而滞后期 s>q时ACF几乎为0,那么可以认为真实的数据生成过程是MA(q)。如果样本的偏自相关系数 PACF 不等于0,直到滞后期s=p,而滞后期s>p时PACF几乎为0,那么可以认为真实的数据生 成过程是AR(p)。更一般的情况是,根据样本的ACF 和PACF 的表现,可拟合出一个较合适的ARMA(p,g)模型。表6-1展示了如何确定模型中的参数p和q。

针对时间序列的自相关系数和偏自相关系数的情况,可以使用 plot_acf()函数和 plot_pacf()函 数进行可视化,运行下面的程序可获得时间序列X2的自相关系数和偏自相关系数的情况,得到的 结果如图6-2所示。

## 对随机序列X2进行自相关和偏相关分析可视化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X2"],"r-")
plt.grid()
plt.title("X2序列波动")
ax = fig.add_subplot(1,3,2)
plot_acf(df["X2"], lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(df["X2"], lags=60,ax = ax)
plt.grid()
plt.tight_layout()
plt.show()
 
#在图像中滞后0表示自己和自己的相关性，恒等于1。不用于确定p和q。

 运行结果如下：

图6-2 X2的自相关系数和偏自相关系数

从图6-2中可以发现,针对白噪声的平稳序列,参数p和g的取值均可以为0。

使用下面的程序可以将序列X1进行自相关分析可视化,结果如图6-3所示。

## 对非随机序列X1进行自相关和偏相关分析可视化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X1"],"r-")
plt.grid()
plt.title("X1序列波动")
ax = fig.add_subplot(1,3,2)
plot_acf(df["X1"], lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(df["X1"], lags=60,ax = ax)
plt.ylim([-1,1])
plt.grid()
plt.tight_layout()
plt.show()

 运行结果如下：

图6-3 X1的自相关系数和偏自相关系数 

从图6-3中可以发现,序列X1具有一定的周期性。

针对序列X1一阶差分后的序列,其自相关和偏自相关分析可视化可以使用下面的程序，运行后的结果如图6-4所示。

## 对非随机序列X1一阶差分后的序列进行自相关和偏相关分析可视化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(X1diff,"r-")
plt.grid()
plt.title("X1序列一阶差分后波动")
ax = fig.add_subplot(1,3,2)
plot_acf(X1diff, lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(X1diff, lags=60,ax = ax)
plt.grid()
plt.tight_layout()
plt.show()
 
#ARMA(p,q)中，自相关系数的滞后，对应着参数q；偏相关系数的滞后对应着参数p。

  运行结果如下：

图6-4 X1一阶差分后的自相关系数和偏自相关系数 

从图6-4中可以发现,序列X1一阶差分后同样具有一定的周期性。
pm.arima模块的decompose()函数可以对时间序列数据进行分解,使用参数multiplicative可以获得乘法模型的分解结果,使用参数additive 可以获得加法模型的分解结果。运行下面的程序 可获得对序列X1乘法模型分解的结果,可视化结果如图6-5所示。

## 时间序列的分解
## 通过观察序列X1，可以发现其既有上升的趋势，也有周期性的趋势，所以可以将该序列进行分解
## 使用乘法模型分解结果(通常适用于有增长趋势的序列)
X1decomp = pm.arima.decompose(df["X1"].values,"multiplicative", m=12)
## 可视化出分解的结果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
                              show=False)
ax[0].set_title("乘法模型分解结果")
plt.show()

 运行结果如下：

图6-5 序列X1的分解结果

通过观察序列X1的分解结果,可以发现其既有上升趋势,也有周期性的变化趋势
使用下面的程序可以对序列X1一阶差分后的序列使用加法模型进行分解,程序运行后的结果如图6-6所示。

## 使用加法模型分解结果(通常适用于平稳趋势的序列)
X1decomp = pm.arima.decompose(X1diff.values,"additive", m=12)
## 可视化出分解的结果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
                              show=False)
ax[0].set_title("加法模型分解结果")
plt.show()

  运行结果如下：

 图6-6 序列X1一阶差分后的分解结果

2 移动平均算法

移动平均算法是一种简单有效的时间序列的预测方法,它的基本思想是:根据时间序列逐项推 移,依次计算包含一定项数的序时平均值,以反映长期趋势。该预测方法中最简单的是简单移动平 均法和简单指数平滑法,较复杂的有霍尔特线性趋势法和Holt-Winters季节性预测模型方法。
本节将使用前面导入的时间序列X1,结合多种移动平均算法对其进行建模与预测,建模时会将 数据后面的24个样本作为测试集,将前面的样本作为训练集,数据切分程序如下,程序运行后的结 果如图6-7所示。训练集中包含120个样本,测试集中包含24个样本。

## 数据准备
## 对序列X1进行切分，后面的24个数据用于测试集
train = pd.DataFrame(df["X1"][0:120])
test = pd.DataFrame(df["X1"][120:])
## 可视化切分后的数据
train["X1"].plot(figsize=(14,7), title= "乘客数量数据",label = "X1 train")
test["X1"].plot(label = "X1 test")
plt.legend()
plt.grid()
plt.show()
 
print(train.shape)
print(test.shape)
df["X1"].shape

运行结果如下：

(120, 1)
(24, 1)
(144,)

 图6-7 训练集和测试集的划分

2.1 简单移动平均法

简单移动平均法中各元素的权重都相等。Python中可以使用时间序列的rolling()和mean()方法进行计算和预测,对切分后的序列进行预测的程序如下,程序中同时将训练集、测试集和预测数 据进行了可视化对比分析, rolling(24).mean()表示计算最近的24个数据的均值,作为待预测数据 的结果,程序运行后的结果如图6-8所示。

## 简单移动平均进行预测
y_hat_avg = test.copy(deep = False)
y_hat_avg["moving_avg_forecast"] =  train["X1"].rolling(24).mean().iloc[-1]
## 可视化出预测结果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["moving_avg_forecast"].plot(style="g--o", lw=2,
                                      label="移动平均预测")
plt.legend()
plt.grid()
plt.title("简单移动平均预测")
plt.show()

运行结果如下：

图6-8 简单移动平均法预测的结果

从图6-8中可以发现,使用简单移动平均法对数据进行预测的效果并不好,使用下面的程序可以计算在测试集上的平均绝对值误差,可知平均绝对值误差为82.55。

## 计算预测结果和真实值的误差
print("预测绝对值误差:",mean_absolute_error(test["X1"],y_hat_avg["moving_avg_forecast"]))

 运行结果如下：

预测绝对值误差: 82.55208333333336

2.2 简单指数平滑法

简单指数平滑又称指数移动平均值,是以指数式递减加权的移动平均。各数据的权重随时间呈 指数式递减,越近期的数据权重越大,但较旧的数据也给予一定的权重。在Python中可以使用 SimpleExpSmoothing()函数对时间序列数据进行简单指数平滑法的建模和预测,对切分后的序列 进行预测的程序如下。在下面的程序中,通过训练获得了两个指数平滑模型,分别对应着参数 smoothing_level=0.15 和 smoothing_level=0.5。同时将训练集、测试集和预测数据进行了对比 可视化,程序运行后的结果如图6-9所示。

## 数据准备
y_hat_avg = test.copy(deep = False)
## 模型构建
model1 = SimpleExpSmoothing(train["X1"].values).fit(smoothing_level=0.15)
y_hat_avg["exp_smooth_forecast1"] = model1.forecast(len(test))
 
model2 = SimpleExpSmoothing(train["X1"].values).fit(smoothing_level=0.5)
y_hat_avg["exp_smooth_forecast2"] = model2.forecast(len(test))
 
## 可视化出预测结果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["exp_smooth_forecast1"].plot(style="g--o", lw=2,
                                      label="smoothing_level=0.15")
y_hat_avg["exp_smooth_forecast2"].plot(style="g--s", lw=2,
                                      label="smoothing_level=0.5")
plt.legend()
plt.grid()
plt.title("简单指数平滑预测")
plt.show()
 
## 计算预测结果和真实值的误差
print("smoothing_level=0.15,预测绝对值误差:",
      mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast1"]))
print("smoothing_level=0.5,预测绝对值误差:",
      mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast2"]))

 运行结果如下：

smoothing_level=0.15,预测绝对值误差: 81.10115706423566
smoothing_level=0.5,预测绝对值误差: 106.813228720506

图6-9 简单指数平滑法预测结果

从输出结果和图6-9中可以发现,参数smoothing_level=0.15获得的模型预测效果,比参数 smoothing_level=0.5 获得的模型预测效果更好。但是使用指数平滑法获得的模型,预测效果仍然 较差。

2.3 霍尔特线性趋势法

霍尔特(Holt)线性趋势法是扩展了的简单指数平滑法,其允许有趋势变化的数据预测,所以 对于有趋势变化的序列可能会获得更好的预测结果。Python中可以使用Holt()函数对时间序列进行霍尔特(Holt)线性趋势法的建模和预测,并且可以使用 smoothing_level 和 smoothing_slope 两个参数控制模型的拟合情况。对切分后的序列进行预测的程序如下,程序中分别训练获得了两个 霍尔特(Holt)线性趋势法模型,对应的参数有 smoothing_level=0.15, smoothing_slope=0.05 和smoothing_level=0.15、 smoothing_slope=0.25。程序中还将训练集、测试集和预测数据进 行了可视化,程序运行后的结果如图6-10所示。

## 数据准备
y_hat_avg = test.copy(deep = False)
## 模型构建
model1 = Holt(train["X1"].values).fit(smoothing_level=0.1,
                                      smoothing_slope = 0.05)
y_hat_avg["holt_forecast1"] = model1.forecast(len(test))
 
model2 = Holt(train["X1"].values).fit(smoothing_level=0.1,
                                      smoothing_slope = 0.25)
y_hat_avg["holt