目录
时间序列数据是常见的数据类型之一,时间序列分析基于随机过程理论和数理统计学方法,研究时间序列数据所遵从的统计规律,常用于系统描述、系统分析、预测未来等。
时间序列数据主要是根据时间先后,对同样的对象按照等时间间隔收集的数据,比如每日的平均气温、每天的销售额、每月的降水量等。虽然有些序列所描述的内容取值是连续的,比如气温的变化可能是连续的,但是由于观察的时间段并不是连续的,所以可以认为是离散的时间序列数据。
一般地,对任何变量做定期记录就能构成一个时间序列。根据所研究序列数量的不同,可以将时间序列数据分为一元时间序列数据和多元时间序列数据。
时间序列的变化可能受一个或多个因素的影响,导致它在不同时间的取值有差异,这些影响因素分别是长期趋势、季节变动、循环波动(周期波动)和不规则波动(随机波动)。时间序列分析主要有确定性变化分析和随机性变化分析。确定性变化分析包括趋势变化分析、周期变化分析、循环变化分析。随机性变化分析主要有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)利用已通过检验的模型进行预测。
对于时间序列数据,最重要的检验就是时间序列数据是否为白噪声数据、时间序列数据是否平稳,以及对时间序列数据的自相关系数和偏自相关系数进行分析。如果时间序列数据是白噪声数据, 说明其没有任何有用的信息。针对时间序列数据的很多分析方法,都要求所研究的时间序列数据是平稳的,所以判断时间序列数据是否平稳,以及如何将非平稳的时间序列数据转化为平稳序列数据,对时间序列数据的建模研究是非常重要的。
本节将会利用两个时间序列数据进行相关的检验分析,首先读取数据并使用折线图将两组时间序列进行可视化,运行下面的程序后,结果如图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,说明该序列为白噪声,没有分析价值。
时间序列是否是平稳的,对选择预测的数学模型非常关键。如果一组时间序列数据是平稳的, 就可以直接使用自回归移动平均模型(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。
自相关分析和偏自相关分析,是用来确定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一阶差分后的分解结果