时间序列平稳性相关检验方法

理解平稳性

一般来说，平稳时间序列是指随着时间的推移具有相当稳定的统计特性的时间序列，特别是在均值和方差方面。平稳性可能是一个比较模糊的概念，将序列排除为不平稳可能比说序列是平稳的更容易。通常不平稳序列有几个特征：

平均值随时间推移发生变化，不保持稳定；
每个周期的波峰和波谷之间的距离在增长（缩短），也就是序列的方差随时间推移而增加（减小）；
整个序列表现出强烈的季节性行为，与平稳性相反。

如何判断序列是否平稳？

时间序列平稳性检验方法，一般可分为两类：

图形分析方法
假设检验方法

1.图形分析方法

这是最简单的方法，即可视化时间序列数据，或可视化时间序列的统计特征，进行肉眼判断。

1.1 可视化数据

import numpy as np
import pandas as pd
import akshare as ak
import matplotlib.pyplot as plt
import seaborn as sns
import warnings

sns.set(style='darkgrid')
warnings.filterwarnings('ignore')

np.random.seed(123)

# -------------- 准备数据 --------------
# 白噪声
white_noise = np.random.standard_normal(size=1000)

# 随机游走
x = np.random.standard_normal(size=1000)
random_walk = np.cumsum(x)

# GDP
df = ak.macro_china_gdp()
df = df.set_index('季度')
# df.index = pd.to_datetime(df.index)
gdp = df['国内生产总值-绝对值'][::-1].astype('float')

# GDP DIFF
gdp_diff = gdp.diff(4).dropna()


# -------------- 绘制图形 --------------
fig, ax = plt.subplots(2, 2, figsize=(22,16), dpi=300)

ax[0][0].plot(white_noise)
ax[0][0].set_title('white_noise')
ax[0][1].plot(random_walk)
ax[0][1].set_title('random_walk')

ax[1][0].plot(gdp)
ax[1][0].set_title('gdp')
ax[1][1].plot(gdp_diff)
ax[1][1].set_title('gdp_diff')

plt.show()

在这里插入图片描述

白噪声：曲线围绕0值上下波动，波动幅度前后、上下一致，为平稳序列；
随机游走：曲线无确定趋势，均值、方差波动较大，非平稳序列；
GDP数据趋势上升，均值随时间增加，非平稳序列；
GDP季节差分后数据，曲线大致在一条水平线上下波动，波动幅度前后变化较小，可认为是平稳的。

1.2 可视化统计特征

可视化统计特征，是指绘制时间序列自相关图和偏自相关图，根据自相关图的表现来判断序列是否平稳。

自相关，也叫序列相关，是一个信号与自身不同时间点的相关度，或者说与自身的延迟拷贝（或滞后性）的相关性，是延迟的函数。不同滞后期得到的自相关系数，叫自相关图。（这里有一个默认假设是序列是平稳序列，平稳序列的自相关性只和时间间隔k有关，不随时间t的变化而变化，因而可以称自相关函数是延迟k的函数）

平稳序列通常具有短期自相关性，对于平稳的时间序列，自相关系数往往会迅速退化到零（滞后期约短相关性越高，滞后期为0，相关性为1）；而对于非平稳的数据，退化会发生得更慢，或存在先减后增或者周期性波动等变动。

自相关的计算公式为根据滞后期 $k$ 将序列拆成等长的两个序列，计算两个序列的相关性得到滞后期为 $k$ 时的自相关性。举例：
$X=[2,3,4,3,8,7],\quad A=[2,3,4,3,8], \quad B=[3,4,3,8,7]$
$\bar{X} = \frac{1}{6} \sum_{i=1}^{6}X_i = 4.5$
$s^2(X)=\frac{1}{6} \sum_{i=1}^{6}(X_i- \bar{X})(X_i- \bar{X})=4.916$
$r(1)=\frac{1}{5}(A_i-\bar{X})(B_i-\bar{X})=1.75$
$ACF(1)=\frac{r(1)}{s^2(X)}=0.356$

import statsmodels.api as sm
X = [2,3,4,3,8,7]
sm.tsa.stattools.acf(X, nlags=1, adjusted=True)    # 返回值：[滞后期为0时的自相关性，滞后期为1的自相关系数]

array([1.       , 0.3559322])

根据 $A CF$ 求出滞后 $k$ 自相关系数时，实际上得到并不是 $X (t)$ 与 $X (t - k)$ 之间单纯的相关关系。因为 $X (t)$ 还会受到中间 $k - 1$ 个随机变量的影响，而这 $k - 1$ 个随机变量又都和 $X (t - k)$ 具有相关关系，如果提出中间 $k - 1$ 个随机变量的影响，所计算得到的是偏自相关系数，计算较为复杂。

# 数据生成过程在第一个代码块中
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, ax = plt.subplots(4, 2, figsize=(22, 16), dpi=300)
fig.subplots_adjust(hspace=1.2)

plot_acf(white_noise, ax=ax[0][0])
ax[0][0].set_title('ACF(white_noise)')
plot_pacf(white_noise, ax=ax[0][1])
ax[0][1].set_title('PACF(white_noise)')

plot_acf(random_walk, ax=ax[1][0])
ax[1][0].set_title('ACF(random_walk)')
plot_pacf(random_walk, ax=ax[1][1])
ax[1][1].set_title('PACF(random_walk)')

plot_acf(gdp, ax=ax[2][0])
ax[2][0].set_title('ACF(gdp)')
plot_pacf(gdp, ax=ax[2][1])
ax[2][1].set_title('PACF(gdp)')

plot_acf(gdp_diff, ax=ax[3][0])
ax[3][0].set_title('ACF(gdp_diff)')
plot_pacf(gdp_diff, ax=ax[3][1])
ax[3][1].set_title('PACF(gdp_diff)')

plt.show()

在这里插入图片描述

白噪音的自相关系数很快就衰减到0附近，是明显的平稳序列。滞后期为0时的自相关系数和偏自相关系数其实就是序列自己和自己的相关性，故为1；滞后期为1时，自相关系数为0，表示白噪声无自相关性；
随机游走，自相关系数下降非常缓慢，故为非平稳序列；另从偏自相关系数中可以看到随机游走只和前一项有关；
GDP数据的自相关图中也可以看到存在一定的周期性，滞后4、8、12等自相关系数较大下降较慢，差分后下降多一些起到一定效果，认为差分后序列是平稳的。同可视化数据一样，直观判断带有较强主观性，但能让我们对数据有更直观的认识。

2.假设检验

平稳性的假设检验方法当前主流为单位根检验，检验序列中是否存在单位根，若存在，则为非平稳序列，不存在则为平稳序列。

2.1 DF检验

$A D F$ 检验(Augmented Dickey-Fuller Testing)是最常用的单位根检验方法之一，通过检验序列是否存在单位根来判断序列是否是平稳的。 $A D F$ 是 $D F$ 检验的增强版。

迪基（Dickey）和弗勒（Fuller）1979年基于非平稳序列的基本特征将其大致归为三类并提出DF检验：

当序列基本走势呈现无规则上升或下降并反复时，将其归为无漂移项自回归过程；
当序列基本走势呈现明显的随时间递增或递减且趋势并不太陡峭时，将其归为带漂移项自回归过程；
当序列基本走势随时间快速递增时，则将其归为带趋势项回归过程。

对应检验回归式为：

无漂移项自回归过程： $Y_t=\rho Y_{t-1} + \varepsilon_t, \quad Y_0=0$
带漂移项自回归过程： $Y_t=\mu + \rho Y_{t-1} + \varepsilon_t, \quad Y_0=0$
带漂移项和趋势项自回归过程： $Y_t=\mu + \beta_t + \rho Y_{t-1} + \varepsilon_t, \quad Y_0=0$
其中 $\mu$ 是常数项， $\beta_t$ 是趋势项， $\varepsilon_t$ 为白噪声无自相关性。

原假设 $H_0: \rho = 1$ （存在单位根，为非平稳时间序列）；
备择假设 $H_1: \rho < 1$ （不存在单位根，时间序列是平稳的：不含截距项和趋势项的平稳/含截距项的平稳/含截距项和趋势平稳）；

若检验统计量大于临界值（ $P$ 值大于显著性水平 $\alpha$ ），不能拒绝原假设，序列是非平稳的；
若检验统计量小于临界值（ $P$ 值小于显著性水平 $\alpha$ ），拒绝原假设，序列是平稳的。

2.2 ADF检验

$D F$ 检验公式为一阶自回归过程，为了能适用于高阶自回归过程的平稳性检验，迪基等1984年对 $D F$ 检验进行了一定的修正，引入了更高阶的滞后项， $A D F$ 的检验回归式修正为：

无漂移项自回归过程： $Y_t=\rho Y_{t-1} + \sum_{i=1}^{k}C_i\bigtriangleup Y_{t-i} + \varepsilon_t, \quad Y_0=0$
带漂移项自回归过程： $Y_t=\mu + \rho Y_{t-1} + \sum_{i=1}^{k}C_i\bigtriangleup Y_{t-i} + \varepsilon_t, \quad Y_0=0$
带漂移项和趋势项自回归过程： $Y_t=\mu + \beta_t + \rho Y_{t-1} + \sum_{i=1}^{k}C_i\bigtriangleup Y_{t-i} + \varepsilon_t, \quad Y_0=0$
其中 $\mu$ 是常数项， $\beta_t$ 是趋势项， $\varepsilon_t$ 为白噪声无自相关性。

假设条件不变：

原假设 $H_0: \rho = 1$ （存在单位根，为非平稳时间序列）；
备择假设 $H_1: \rho < 1$ （不存在单位根，时间序列是平稳的：不含截距项和趋势项的平稳/含截距项的平稳/含截距项和趋势平稳）；

检验流程同 $D F$ 检验一致。若要严格判断序列是否是宽平稳的，可以直接检验是否不含截距项和趋势项平稳；若不能拒绝原假设（如 $p > 0.05$ ），序列非平稳，此时仍有必要检验序列是否是趋势平稳的。非平稳且非趋势平稳，可以使用一阶差分等平稳化方法处理后再做检验，弱势趋势平稳，困于过度差分则不宜使用差分方式平稳化。

import numpy as np
from matplotlib import pyplot as plt

np.random.seed(123)

y = np.random.standard_normal(size=100)
for i in range(1, len(y)):
    y[i] = 1 + 0.1*i + y[i]

plt.figure(figsize=(12, 6), dpi=300)
plt.plot(y)
plt.show()

在这里插入图片描述

from arch.unitroot import ADF
adf = ADF(y)
print(adf.summary().as_text())

adf = ADF(y, trend='ct')
print(adf.summary().as_text())

   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -0.739
P-value                         0.836
Lags                                5
-------------------------------------

Trend: Constant
Critical Values: -3.50 (1%), -2.89 (5%), -2.58 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -9.963
P-value                         0.000
Lags                                0
-------------------------------------

Trend: Constant and Linear Time Trend
Critical Values: -4.05 (1%), -3.46 (5%), -3.15 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

arch包中 $A D F$ 检验可指定trend = "n"（不含截距项和时间趋势项）、trend = "c"（含截距项）、trend = "ct"（含截距项和时间趋势项）、trend = "ctt"（含截距项、时间趋势项和二次型时间趋势项），分别对应不同平稳类型的检验。

再来看看GDP季节差分前后数据是否为平稳：

# 数据在第一个代码块中
from arch.unitroot import ADF
adf = ADF(gdp)
print(adf.summary().as_text())

adf = ADF(gdp_diff)            # 差分后平稳
print(adf.summary().as_text())

   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                  2.606
P-value                         0.999
Lags                               12
-------------------------------------

Trend: Constant
Critical Values: -3.55 (1%), -2.91 (5%), -2.59 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -3.651
P-value                         0.005
Lags                                8
-------------------------------------

Trend: Constant
Critical Values: -3.55 (1%), -2.91 (5%), -2.59 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

# 数据在第一个代码块中
from arch.unitroot import ADF
adf = ADF(gdp, trend='ct')        # 差分前非趋势平稳
print(adf.summary().as_text())

   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -2.358
P-value                         0.402
Lags                                4
-------------------------------------

Trend: Constant and Linear Time Trend
Critical Values: -4.10 (1%), -3.48 (5%), -3.17 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

2.3 PP检验

Phillips和Perron(1988)提出了一种非参数检验方法，主要是为了解决残差项中潜在的序列相关和异方差问题，其检验统计量的渐进分布和临界值与ADF检验相同。同样出现较早，假设条件一样，用法相似，可作为ADF检验的补充。

假设条件不变：

原假设 $H_0: \rho = 1$ （存在单位根，为非平稳时间序列）；
备择假设 $H_1: \rho < 1$ （不存在单位根，时间序列是平稳的：不含截距项和趋势项的平稳/含截距项的平稳/含截距项和趋势平稳）；

同样构造一个趋势平稳序列，看下PP检验结果：

import numpy as np
from arch.unitroot import PhillipsPerron

np.random.seed(123)

y = np.random.standard_normal(size=100)
for i in range(1, len(y)):
    y[i] = 1 + 0.1*i + y[i]


pp = PhillipsPerron(y)             # 10%显著水平下平稳
print(pp.summary().as_text())

pp = PhillipsPerron(y, trend='ct')   # 1%显著水平下趋势平稳
print(pp.summary().as_text())

     Phillips-Perron Test (Z-tau)    
=====================================
Test Statistic                 -2.825
P-value                         0.055
Lags                               12
-------------------------------------

Trend: Constant
Critical Values: -3.50 (1%), -2.89 (5%), -2.58 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
     Phillips-Perron Test (Z-tau)    
=====================================
Test Statistic                -10.009
P-value                         0.000
Lags                               12
-------------------------------------

Trend: Constant and Linear Time Trend
Critical Values: -4.05 (1%), -3.46 (5%), -3.15 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

2.4 DF-GLS检验

DF-GLS检验，是Elliot,Rothenberg,and Stock 1996年提出的一种单位根检验方法，全称Dickey-Fuller Test with GLS Detredding，即“使用广义最小二乘法去除趋势的检验”，是目前最有功效的单位根检验。

DF-GLS检验利用广义最小二乘法，首先对要检验的数据进行一次“准差分”，然后利用准差分的数据对原序列进行去除趋势处理，再利用ADF检验的模型形式对去除趋势后的数据进行单位根检验，但此时ADF检验模型中不再包含常数项或者时间趋势变量。

原假设：序列存在单位根（时间序列是非平稳的）；
备择假设：序列不存在单位根（时间序列是平稳的或趋势平稳的）

同样构造一个趋势平稳序列来看下检验效果（DF-GLS检验trend只能指定为c或ct）：

import numpy as np
from arch.unitroot import DFGLS

np.random.seed(123)

y = np.random.standard_normal(size=100)
for i in range(1, len(y)):
    y[i] = 1 + 0.1*i + y[i]

dfgls = DFGLS(y)
print(dfgls.summary().as_text())

dfgls = DFGLS(y, trend='ct')            # 趋势平稳
print(dfgls.summary().as_text())

      Dickey-Fuller GLS Results      
=====================================
Test Statistic                  1.186
P-value                         0.945
Lags                                4
-------------------------------------

Trend: Constant
Critical Values: -2.77 (1%), -2.15 (5%), -1.83 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
      Dickey-Fuller GLS Results      
=====================================
Test Statistic                 -7.565
P-value                         0.000
Lags                                0
-------------------------------------

Trend: Constant and Linear Time Trend
Critical Values: -3.62 (1%), -3.04 (5%), -2.74 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

2.5 KPSS检验

KPSS检验是由Kwiatkowski, Phillips, and Shin 1992年提出的，与之前的检验方法相比，最大的不同点就是它的原假设是平稳序列或趋势平稳序列，而备择假设是存在单位根。

原假设：序列不存在单位根（时间序列是平稳的或趋势平稳的）
备择假设：序列存在单位根（时间序列是非平稳的）

import numpy as np
from arch.unitroot import KPSS

np.random.seed(123)

y = np.random.standard_normal(size=100)
for i in range(1, len(y)):
    y[i] = 0.1 + y[i-1] + y[i]

kpss = KPSS(y)
print(kpss.summary().as_text())

kpss = KPSS(y, trend='ct')
print(kpss.summary().as_text())

    KPSS Stationarity Test Results   
=====================================
Test Statistic                  1.393
P-value                         0.000
Lags                                5
-------------------------------------

Trend: Constant
Critical Values: 0.74 (1%), 0.46 (5%), 0.35 (10%)
Null Hypothesis: The process is weakly stationary.
Alternative Hypothesis: The process contains a unit root.
    KPSS Stationarity Test Results   
=====================================
Test Statistic                  0.114
P-value                         0.115
Lags                                5
-------------------------------------

Trend: Constant and Linear Time Trend
Critical Values: 0.22 (1%), 0.15 (5%), 0.12 (10%)
Null Hypothesis: The process is weakly stationary.
Alternative Hypothesis: The process contains a unit root.

KPSS检验中原假设为不存在单位根。默认检验趋势类型下P值为0.000，拒绝原假设，存在单位根，序列非平稳。指定trend="ct"后，P值不拒绝原假设，认为序列趋势平稳，检验错误。以上几种检验均不能100%保证检验正确，PP检验可认为是ADF检验的补充，KPSS检验同样也可和其他检验一同使用，当均认为是平稳或趋势平稳时方判定为平稳。