Является ли этот временной ряд стационарным или нет?
Я хочу проверить стационарные данные временного ряда, сохраненные в TS.csv.
Тем не менее, R's tseries::adf.test()
и питона statsmodels.tsa.stattools.adfuller()
дать совершенно разные результаты.
adf.test()
показывает, что это стационарно (р < 0,05), в то время как adfuller()
показывает, что он нестационарный (р> 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 различны, и как сделать их одинаковыми, изменив параметры.
Для второго вопроса выше: "Как правильно тестировать стационарные временные ряды в R и Python?". Я хотел бы предоставить некоторые детали:
При прогнозировании временного ряда модели ARIMA необходимо, чтобы входной временной ряд был стационарным. Если вход не стационарный, он должен быть log()
ред или diff()
чтобы сделать его неподвижным, затем поместите его в модель.
Таким образом, проблема заключается в следующем: должен ли я считать, что вход является стационарным (с параметрами по умолчанию R), и вписать его непосредственно в модель ARIMA, или считать его нестационарным (с параметрами по умолчанию в Python), и сделать его стационарным с дополнительными функциями (такими как log()
или же diff()
)?
2 ответа
Результаты отличаются, потому что модели подходят немного отличаются, и потому, что порядок отставания моделей совершенно разные. Тест на питоне включает в себя постоянный термин "дрейф" (постоянная оценивается, таким образом, центрируя временной ряд в ноль), но тест R включает в себя как постоянный, так и линейный тренд. Это может быть указано в коде Python с аргументом regression = 'ct'
,
Задержка по умолчанию в r
nlag = trunc((length(x)-1)^(1/3))
Длина лага по умолчанию в питоне
12*(nobs/100)^(1/4)
Когда вы запустили код Python, вы сказали функции выбирать оптимальную длину лага по критериям AIC. Если мы скажем python запустить центрированную и отклоненную модель и скажем, что она использует критерии длины лага R, мы получим:
In [5]: adfuller(ts, regression="ct", maxlag = 4)[1]
Out[5]: 3.6892966741832268e-09
Трудно понять, совпадает ли это с результатом R, поскольку R округляет свое значение p до 0,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
И в питоне:
In [6]: adfuller(ts, regression="ct")[1]
Out[6]: 0.58756464088883864
Не идеально, но довольно близко.
Замечания:
Фактический тест-статистика Дики-Фуллера для модели питона
In [8]: adfuller(ts, regression="ct")[0]
Out[8]: -2.025340637385288
который идентичен результату R. Тесты, вероятно, используют разные способы вычисления p-значения из статистики.
P-значения теста Аугментированного Дики-Фуллера довольно чувствительны к выбору порядка запаздывания. Например, вот тот же тест в 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 по-разному определяются в двух пакетах.