Stationary test (정상성 검정) python

시계열 데이터를 다루게 된다면, 정상성 검정이라는 것을 시행해야 한다.

 

시계열 데이터를 통해 회귀를 하게 된다면, 이는 과거 데이터를 가지고 미래를 예측하는 것과 같다. 

 

따라서 통계적 속성이 일정해야 미래 데이터의 예측에 대한 신뢰성이 보장된다고 할 수 있을 것이다.

 

다음 그래프를 살펴보자

https://en.wikipedia.org/wiki/Stationary_process

1번째 그래프는 정상성이 보장된 그래프이다. 이를 보면, 데이터의 분산이 평균을 중심으로 시간의 흐름에 따라 일정하게 유지되는 것을 확인할 수 있다.

반면 2번째 그래프는 분산이 일정하지 않고 추세도 하향하는 것을 볼 수 있다.

그렇다면, 회귀분석을 통해 미래를 예측하는 것은 아무래도 1번째 그래프가 더 쉬워보이는 것을 육안으로도 확인할 수 있을 것이다.

 

하지만, 표를 그려서 확인하는 것보다는 통계적인 검정을 활용하여 정상성을 검정하는 것이 더 명확한 방법이다.

따라서 이번에는 정상성을 검정하는 방식인 ADF, KPSS 검정을 시행해 보겠다.

 

먼저 데이터는 divvy_data로 다음 사이트에서 불러온 데이터이다.

https://ride.divvybikes.com/system-data

Divvy System Data | Divvy Bikes

divvy_daily.csv

0.05MB

 

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.stattools import adfuller, kpss

data = pd.read_csv('https://raw.githubusercontent.com/signature95/tistory/main/dataset/divvy_daily.csv', index_col = 'date')
data.head()

출력 결과

          rides  temperature       weather
date                                      
1/1/2014     95    19.483158  rain or snow
1/2/2014    111    16.833333  rain or snow
1/3/2014      6    -5.633333         clear
1/4/2014    181    30.007735  rain or snow
1/5/2014     32    16.756250  rain or snow

 

이번에는 다른 featrue 없이 rides 만 활용하여 시계열 분석을 할 것이다. 

그렇다면 도표를 그려보고 rides의 변화를 살펴보자

plt.figure(figsize=(16,9))
data.rides.plot(kind='line')
plt.ylabel('divvy_rider')
plt.show()

육안으로 확인해보면, 계절성이 존재하고 추세도 우상향하는 형태로 보인다. 따라서 이번에는 ADF, KPSS 검정을 통해 통계적으로 확인해보자.

 

 

먼저 adf 검정이다.

# 함수 형성
def adf_test(timeseries, pvalue = .05, regression_option = 'ct'):
    print ('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC', regression = regression_option)
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
       dfoutput['Critical Value (%s)'%key] = value
    print (dfoutput)
    if dfoutput[1] < pvalue:
       print(f"정상시계열이 아니라는 귀무가설을 {pvalue*100}%의 유의수준으로 기각할 수 있으므로 해당 데이터는 정상성이 보장됩니다.")
    else:
       print(f"정상시계열이 아니라는 귀무가설을 {pvalue*100}%의 유의수준으로 기각할 수 없으므로 해당 데이터는 정상성을 보장하지 못합니다.")

# adf 함수를 통해 출력
adf_test(data['rides'])
       
>>>
       
Results of Dickey-Fuller Test:
Test Statistic                   -1.683520
p-value                           0.758045
Lags Used                        22.000000
Number of Observations Used    1435.000000
Critical Value (1%)              -3.965093
Critical Value (5%)              -3.413554
Critical Value (10%)             -3.128854
dtype: float64
정상시계열이 아니라는 귀무가설을 5.0%의 유의수준으로 기각할 수 없으므로 해당 데이터는 정상성을 보장하지 못합니다.

adf에는 여러 파라미터가 존재하는데 먼저 여기서 regression option = 'ct'로 한 부분에 대해서 설명하고자 한다.

먼저 adf에는 4가지의 옵션이 존재한다. c, ct, ctt, n으로 구성되어있고 c가 default로 되어 있다.

• c : constant (상수항만 존재)
• ct : constant and trend (상수항과 추세 존재)
• ctt : constant and linear and quardratic trend (상수항과 1차 2차 추세가 존재)
• n : no constant, no trend (상수항과 추세가 둘다 존재하지 않음)

추세? 상수항? 이를 더 명확하게 알기 위해서는 다음식을 확인해보면 된다.

AR(1)

AR은 Autoregressive로 자기상관성을 의미한다. AR(1)에서 1은 t시점 기준으로 직전의 값과 현재의 값이 white noise(백색잡음)으로 이뤄져 있다는 것을 의미한다. 즉 현 시점의 관측치는 직전 시점의 관측치에 현 시점의 에러값(white noise)로 구성되 있는 것이다. 여기서 상수항이 constant를 의미한다. 정리하면 위 모형을 간단하게 random walk 모형이라고도 한다.

 

계수 설명

a : 현 시점에 발생한 새로운 충격 혹은 잡음으로 E(a_t) = 0, COV(a_t, a_s) = 0을 만족한다.
ø : 전기(t-1)의 충격이 현재(t) 시점에 미치는 영향을 계수로 표현한 것.
- ø = 0 인 경우
X_t = a 를 만족한다 
- |ø| <1 인 경우
이를 만족하지 못하면 모형자체가 성립하지 않음. (불만족하는 경우 비정상 시계열임)
ø = 0.5일때는 t-1 기의 충격이 다음기에는 50% 감소하게 됨.
일반적으로 전기의 충격이 더 커지는 성향은 없다는 전제하에 모형을 가정함. 또한, ø이 1 이상인 경우, 또한, X_t의 분산이 ∞으로 커지게 됨

 

 

즉 우리는 귀무가설을 ø = 1로 설정하고 이를 검정하기 위해 ADF를 진행한다.

AR(1)

여기 식에서 ADF를 진행하게 된다면, regression option을 'c'로 설정해야 한다. trend로 표현된 것이 없기 때문이다.

그렇다면 추세가 어느정도 존재하는 것처럼 보인다면 식은 다음과 같이 쓸 수 있다.

ßt 로 표현된 값이 trend이다. 이경우는 detrend를 통해 time을 제거하면 되는데 이는 다음에 다루도록 한다.

 

다시 돌아와서 regression option을 'ct'로 설정하였고, 그에대해 검정해 보았는데 결과는 귀무가설을 기각하지 못했기에 해당 데이터가 non-stationary하다는 사실을 알 수 있었다. 

이번에는 KPSS 검정을 해보겠다.

 

def kpss_test(timeseries, pvalue = .05, regression_option = 'ct'):
    print ('Results of KPSS Test:')
    kpsstest = kpss(timeseries, regression= regression_option)
    kpss_output = pd.Series(kpsstest[0:3], index=['Test Statistic','p-value','Lags Used'])
    for key,value in kpsstest[3].items():
        kpss_output['Critical Value (%s)'%key] = value
    print (kpss_output)
    if kpss_output[1] < pvalue:
        print(f"정상시계열이 맞다는 귀무가설을 {pvalue*100}%의 유의수준으로 기각할 수 있으므로 해당 데이터는 정상성을 보장하지 못합니다.")
    else:
        print(f"정상시계열이 맞다는 귀무가설을 {pvalue*100}%의 유의수준으로 기각할 수 없으므로 해당 데이터는 정상성이 보장됩니다.")
        

kpss_test(data['rides'])

>>>

Results of KPSS Test:
Test Statistic            0.103919
p-value                   0.100000
Lags Used                24.000000
Critical Value (10%)      0.119000
Critical Value (5%)       0.146000
Critical Value (2.5%)     0.176000
Critical Value (1%)       0.216000
dtype: float64
정상시계열이 맞다는 귀무가설을 5.0%의 유의수준으로 기각할 수 없으므로 해당 데이터는 정상성이 보장됩니다.
/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/statsmodels/tsa/stattools.py:2015: InterpolationWarning: The test statistic is outside of the range of p-values available in the
look-up table. The actual p-value is greater than the p-value returned.

statsmodel 라이브러리에서 지원하는 KPSS검정에서는 P값이 출력될때 최대 10%에서 최소 1%로 출력이된다.

그렇기 때문에 실제 출력 결과에서 warning 메세지가 뜨고 실제 P값은 출력된 P값보다 더 큰 값이다라는 정보를 주게 된다.

통계치를 확인해보더라도 10%의 유의수준 값이 0.119지만 검정 통계치는 0.103으로 10% 유의수준 통계치보다 작은 것을 확인할 수 있다.

 

여기서도 regression option을 ct로 지정했는데, 원래 KPSS 검정의 기본 값은 c로 설정되어 있다. (하지만 ADF와 다르게 'c', 'ct'의 옵션만 존재한다)

 

이번에 실시한 KPSS의 결과는 ADF와 상반되는 결과를 보여주었다. ADF는 정상성을 보장하지 못하고 KPSS는 정상성을 보장하는 결과를 가져왔는데, 이를 더 확실하게 하기 위해서는 차분(differential)을 사용해야 한다.

 

다음글에서 이어집니다.

2022.01.19 - [공부/모델링] - Differential (차분) python

Differential (차분) python