这个时间序列是否平稳?

Is this time series stationary or not?

我想检查 TS.csv 中保存的时间序列数据的平稳性。

但是,R 的 tseries::adf.test() 和 Python 的 statsmodels.tsa.stattools.adfuller() 给出完全不同的结果。

adf.test() 表示它是平稳的 (p < 0.05),而 adfuller() 表示它是非平稳的 (p > 0.05)。

下面的代码有没有问题?

在 R 和 Python 中测试时间序列的平稳性的正确过程是什么?

谢谢。

R 代码:

> rd <- read.table('Data/TS.csv', sep = ',', header = TRUE)
> inp <- ts(rd$Sales, frequency = 12, start = c(1965, 1))
> inp
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1965 154  96  73  49  36  59  95 169 210 278 298 245
1966 200 118  90  79  78  91 167 169 289 347 375 203
1967 223 104 107  85  75  99 135 211 335 460 488 326
1968 346 261 224 141 148 145 223 272 445 560 612 467
1969 518 404 300 210 196 186 247 343 464 680 711 610
1970 613 392 273 322 189 257 324 404 677 858 895 664
1971 628 308 324 248 272
> library(tseries)
> adf.test(inp)

    Augmented Dickey-Fuller Test

data:  inp
Dickey-Fuller = -7.2564, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Python 代码(来自 Time_Series.ipynb):

import pandas as pd
from statsmodels.tsa.stattools import adfuller
df = pd.read_csv('Data/TS.csv')
ts = pd.Series(list(df['Sales']), index=pd.to_datetime(df['Month'],format='%Y-%m'))
s_test = adfuller(ts, autolag='AIC')
print("p value > 0.05 means data is non-stationary: ", s_test[1])
# output: p value > 0.05 means data is non-stationary:  0.988889420517

更新

@gfgm 很好地解释了为什么 R 和 Python 的结果不同,以及如何通过更改参数使它们相同。

对于上面的第二个问题:"What's the right process to test stationary of a time series in R and Python?"。 我想提供一些细节:

在预测时间序列时,ARIMA 模型需要输入的时间序列是平稳的。 如果输入不是静止的,应该 log()ed 或 diff()ed 使其静止, 然后将其拟合到模型中。

所以问题是: 我是否应该认为输入是固定的(使用 R 的默认参数)并将其直接拟合到 ARIMA 模型中, 或者认为它是非固定的(使用 Python 的默认参数), 并使用额外的功能使其静止不动(如 log()diff())?

结果不同是因为拟合的模型略有不同,而且模型的滞后阶数完全不同。 python 测试包括常数 'drift' 项(估计常数,从而使时间序列以零为中心),但 R 测试包括常数和线性趋势项。这可以在 python 代码中使用参数 regression = 'ct'.

指定

r 中的默认滞后长度

nlag = trunc((length(x)-1)^(1/3))

python

中的默认滞后长度

12*(nobs/100)^(1/4)

当您 运行 python 代码时,您告诉函数根据 AIC 标准选择最佳 lag-length。如果我们告诉 python 到 运行 一个居中和去趋势的模型,并且我们告诉它使用 R lag-length 标准,我们得到:

In [5]: adfuller(ts, regression="ct", maxlag = 4)[1]
Out[5]: 3.6892966741832268e-09

很难看出这是否与 R 的结果相同,因为 R 将其 p-value 舍入为 .01,但我们可以告诉 R 使用 python 的滞后长度,并且 python 使用 R 的模型(我不能用这个函数改变 R 中的模型)。我们得到:

adf.test(inp, k = ceiling(12*(length(inp)/100)^(1/4)))

    Augmented Dickey-Fuller Test

data:  inp
Dickey-Fuller = -2.0253, Lag order = 12, p-value = 0.5652
alternative hypothesis: stationary

在python中:

In [6]: adfuller(ts, regression="ct")[1]
Out[6]: 0.58756464088883864

不完美,但非常接近。

注:

python 模型的实际 Dickey-Fuller test-statistic 是

In [8]: adfuller(ts, regression="ct")[0]
Out[8]: -2.025340637385288

这与 R 结果相同。这些测试可能使用不同的方式从统计数据计算 p-value。

A​​ugmented Dickey-Fuller 检验的 p-values 对滞后阶数的选择相当敏感。例如,这是 R 中具有更高滞后阶数的相同测试:

> adf.test(rd$Sales, k=9)

    Augmented Dickey-Fuller Test

data:  rd$Sales
Dickey-Fuller = -2.9186, Lag order = 9,
p-value = 0.2004
alternative hypothesis: stationary

adf.test 的文档说它使用具有恒定线性趋势的回归。我们应该将参数 regression = 'ct' 传递给 adfuller 以使用相同的回归方法。

目前我的机器上的 statsmodels 遇到了一些问题,但我建议你尝试以下参数,看看你是否得到更接近的对应关系:

adfuller(a, maxlag=9, autolag=None, regression='ct')

您要查找的是两者是否显示相同的测试统计信息,因为 p-values 在两个包之间的确定方式不同。