[시계열 분석] 10. Phillips-Perron Test(검정) with Python

이번 포스팅은 단위근 검정의 한 방법으로 (Augmented) Dickey-Fuller(ADF)와 쌍벽을 이루는 Phillips-Perron Test(검정)에 대해서 알아본다. 그리고 파이썬을 이용하여 샘플 데이터를 통한 예제도 알아보고자 한다.

이번 포스팅에서 다룰 내용은 다음과 같다.

1. Phillips-Perron Test(검정)

2. Python 예제 - Phillips-Perron Test(검정)

본 포스팅에서는 수식을 포함하고 있습니다.

티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니

PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다.

1. Phillips-Perron Test(검정)

1.1 정의

다음과 같은 시계열이 있다고 하자.

$\Delta X_t = c_t + (\rho - 1)X_{t-1} + W_t, t=1,2, \ldots, T$

$W_t$ 는 보통은 백색잡음이지만 꼭 그럴 필요 없다. 여기서 $W_t$ 는 등분산성은 만족하지 않을 수도 있고 약한 자기 상관성을 가질 수도 있다. 그리고 $c_t$ 는 다음과 같이 3종류 형태가 있는 트렌드 함수이다.

$\begin{align} (A) \text{ } c_t &= 0 \text{ no time trend, no intercept} \\ (B) \text{ } c_t &= a \text{ intercept but no time trend} \\ (C) \text{ }c_t &= a_0+a_1t \text{ intercept with time trend} \end{align}$

Phillips-Perron Test(검정)의 귀무가설과 대립 가설은 다음과 같다.

$H_0 : \rho = 1 \text{ vs } H_1 : \rho < 1$

귀무가설은 해당 시계열이 단위근이 있다는 것을 의미하며 이는 비정상 시계열임을 뜻한다. 그리고 대립 가설은 부분 정상성(Local Stationary)을 만족하는 시계열이다.

Phillips-Perron Test(검정)은 최소 제곱법을 이용하여 추정한 $\hat{\rho}$ 과 그에 대응하는 $t$ 통계량 $\tau = \frac{\hat{\rho}-1}{SE(\hat{\rho})}$

을 이용하여 검정한다. 하지만 $\hat{\rho}$ 과 $\tau$ 를 바로 이용하는 것이 아니라 $Z$ 통계량이라는 것을 사용한 통계량을 검정 통계량으로 사용한다. 그 이유를 설명하기 위하여 다음을 정의한다.

$\sigma^2 = \lim_{t\rightarrow\infty}E(T^{-1}S_T^2), S_t = \sum_{j=1}^tW_j\tag{1.1}$

$\sigma_W^2 = \lim_{t\rightarrow\infty}E(W_t^2)\tag{1.2}$

$Z$ 통계량을 사용한 이유를 알아보자. 기존 $\hat{\rho}$ 과 $\tau$ 의 기각값을 계산하기 위하여 귀무 가설이 참일 때 극한 분포를 계산해야한다. 근데 이 극한 분포가 nuisance 모수( $\sigma^2, \sigma_W^2$ )에 의존하게 된다. 이렇게 되면 안그래도 샘플이 큰 경우에 그나마 정확한 기각값을 구할 수 있는데 nuisance 모수가 포함되어 있어 이를 추정해서 플러그인(Plug-In)하면 기각값의 정확성이 떨어지게 된다. 또한 nuisance 모수가 있으면 기각값의 테이블을 만들수 없다. 왜냐하면 nuisance 모수의 추정값을 미리 알 수 없기 때문이다. 이러한 연유로 Phillips와 Perron은 이러한 nuisance 모수에 의존하지 않는 변환을 고안하여 $Z$ 통계량을 만든 것이다. 이게 Phillips-Perron Test(검정)의 핵심이다.

1.2 검정 방법

Phillips-Perron 검정 방법은 다음과 같다.

1) Phillips-Perron 검정에서 검정 통계량을 계산한다.

검정 통계량은 다음과 같다.

$Z_{\hat{\rho}} = T(\hat{\rho}-1) - \frac{1}{2}\frac{T^2 \cdot SE(\hat{\rho})}{\hat{s}^2}(\hat{\sigma}^2-\hat{\sigma}_W^2)\tag{1.3}$

$Z_{\tau} = \left(\frac{\hat{\sigma}_W^2}{\hat{\sigma}^2}\right)^{1/2}\tau- \frac{T(\hat{\sigma}^2-\hat{\sigma}_W^2)SE(\hat{\rho})}{2(\hat{\sigma}^2\hat{s}^2)^{1/2}} \tag{1.4}$

여기서 $\hat{\sigma}^2, \hat{\sigma}_W^2$ 은 식 (1.1), (1.2)에서 $\sigma^2, \sigma_W^2$ 의 추정량으로 다음과 같이 계산한다.

$\hat{\sigma}^2 = T^{-1}\sum_{t=1}^T\hat{w}_t^2+2T^{-1}\sum_{s=1}^la_{sl}\sum_{t=s+1}^T\hat{w}_t\hat{w}_s\tag{1.5}$

이때 $\hat{w}_t, t=1, 2, \ldots, T$ 는 회귀 모형으로 얻어진 잔차이며

$a_{sl} = 1-s/(l+1)$

이다. 그리고 $\hat{\sigma}_W^2 = \sum_{t=1}^T\hat{w}_t^2$ , $\hat{s}^2$ 은 평균오차제곱(MSE)이다(자유도를 분모로 한다).

여기서 궁금증이 생긴다. 왜 통계량 두 개를 만든 걸까?

논문을 봤는데 정확한 이유는 찾지 못했으나 개인적인 생각으로 기존 단위근 검정(Said, Dickey) 보다 좋다는 것을 실험하기 위하여 하나보다는 두 개가 나으니까 그런 것 같다. 실제로 Phillips와 Perron 논문에서 어느 하나가 안 좋으면 다른 하나는 좋다는 결과를 보여준다.

2) 기각 여부 결정

이제 검정 통계량을 계산했으면 테이블을 이용하여 (A)~(C)에 대응하는 기각값을 선택한다.

테이블은 다음과 같다.

pp_test_table.xlsx

0.01MB

기각값을 찾았으면 검정 통계량과 비교한다. 이때 검정 통계량이 기각값보다 작으면 귀무가설을 기각한다.

1.3 장단점

Phillips-Perron의 장점은 다음과 같다.

1) 장점

a. 좀 더 폭넓은 오차항 시계열을 다룰 수 있다.

Phillips-Perron 검정은 앞서 이야기했듯이 오차항은 반드시 백색 잡음 또는 정상성을 가질 필요는 없으며 자기 상관성을 갖거나 시간에 따른 등분산성을 만족하지 않아도 사용할 수 있는 검정방법이다. 또한 등분산성을 나타내는 수학적인 형태의 Robust 한 방법이다.

b. 자기 회귀 차수를 신경 쓰지 않아도 된다.

Phillips-Perron 검정은 Augmented Dickey Fuller과 달리 자기 회귀 차수를 신경쓰지 않아도 된다.

2) 단점

a. 자기 상관 차수를 정해야 한다.

Phillips-Perron 검정은 $\sigma^2$ 추정 시 자기 상관 차수 $l$ 을 정해줘야 한다. 주로

$l \approx 4\left(\frac{T}{100} \right)^{1/4}$

또는

$l \approx 12\left(\frac{T}{100} \right)^{1/4}$

를 주로 쓴다고 한다.

b. 소표본일 경우의 검정의 정확성이 매우 안 좋다.

소표본일 경우 검정의 정확성이 Augmented Dickey-Fuller 검정보다 매우 안 좋다고 한다.

2. Python 예제 - Phillips-Perron Test(검정)

이번엔 파이썬을 이용하여 Phillips-Perron 검정을 구현해보자.

2.1 함수 만들어보기

먼저 Phillips-Perron 검정을 수행하는 함수를 만들어보자. 여기서 만든 함수는 R 패키지 aTSA의 pp.test 함수를 참고했다.

먼저 필요한 모듈을 임포트한다.

 import pandas as pd
import numpy as np
import math
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')

먼저 시계열 $X_t$ $X_{t-j}, j=1, 2, \ldots$ 를 한 매트릭스에 모아주는 embed 함수를 먼저 만들었다.

 def embed(z, lags):
    res = []
    for i in range(lags):
        temp_z = np.expand_dims(z[lags-1-i:len(z)-i], axis=1)
        res.append(temp_z)
    res = np.concatenate(res, axis=1)
    return res

이제 주인공이 나올 시간이다. Phillips-Perron 검정을 수행하는 함수다. 먼저 검정 통계량 $Z_{\hat{\rho}}, Z_\tau$ 를 계산하는 stat 함수, 테이블을 참조하여 주어진 샘플 수에 대하여 적절히 보간법(Interpolation)을 적용하여 p value를 계산하는 get_pval 함수를 내부에서 정의한다.

그러고 나서 (A)~(C)인 경우의 회귀 모형을 적합하고 이를 이용하여 각 경우의 검정 통계량과 p value를 계산하고 그 결과를 튜플 형식으로 반환한다.

 def phillips_perron_test(x, statistic_type = 'z_rho', lag_short = True):
    def stat(model, m):
        idx = 1
        if m > 1:
            idx = 2
        res1 = model.resid
        rho_hat = model.params[idx-1]
        sig_hat = model.bse[idx-1]
        s2 = np.sum(np.square(res1))/(n-m)
        gamma = []
        for i in range(1, q+2):
            u = embed(res1, i)
            gamma_val = np.sum(u[:, 0]*u[:, i-1])/n
            gamma.append(gamma_val)
 
        lambda2 = gamma[0]+2*np.sum((1-np.array(range(1, q+1))/(q+1))*gamma[1:])
        z_r = n*(rho_hat-1) - ((n**2)*(sig_hat**2)/s2)*(lambda2-gamma[0])/2
        z_t = np.sqrt(gamma[0]/lambda2)*(rho_hat-1)/sig_hat-\
            (lambda2-gamma[0])*n*sig_hat/(2*np.sqrt(lambda2*s2))
        return (z_r, z_t)
 
    def get_pval(table, stat):
        from scipy import interpolate
        table = np.array(table)
        ncol = table.shape[1]
        size = (25, 50, 100, 250, 500, 1000)
        percnt = (0.01, 0.025, 0.05, 0.1, 0.5, 0.9, 0.95, 0.975, 0.99)
        interpol_size = []
        for j in range(ncol):
            f_temp = interpolate.interp1d(size, table[:, j], fill_value='extrapolate')
            interpol_size.append(f_temp(n).item())
 
        f = interpolate.interp1d(interpol_size, percnt, fill_value='extrapolate')
        return f(stat).item()
 
    
    assert statistic_type in ['z_rho', 'z_tau']
    z = embed(x, 2)
    yt = z[:,0]
    yt1 = z[:,1]
    n = len(yt)
    assert n >= 1
    t = range(1, n+1)
    lag = 4
    if not lag_short:
        lag = 12
 
    q = math.ceil(lag*(n/100)**(0.25))
 
    table_r1 = ((-11.8, -9.3, -7.3, -5.3, -0.82, 1.01, 1.41, 1.78, 2.28), 
                (-12.8, -9.9, -7.7, -5.5, -0.84, 0.97, 1.34, 1.69, 2.16), 
                (-13.3, -10.2, -7.9, -5.6, -0.85, 0.95, 1.31, 1.65, 2.09), 
                (-13.6, -10.4, -8, -5.7, -0.86, 0.94, 1.29, 1.62, 2.05), 
                (-13.7, -10.4, -8, -5.7, -0.86, 0.93, 1.29, 1.61, 2.04), 
                (-13.7, -10.5, -8.1, -5.7, -0.86, 0.93, 1.28, 1.6, 2.03))
 
    table_r2 = ((-17.2, -14.6, -12.5, -10.2, -4.22, -0.76, 0, 0.64, 1.39), 
                (-18.9, -15.7, -13.3, -10.7, -4.29, -0.81, -0.07, 0.53, 1.22), 
                (-19.8, -16.3, -13.7, -11, -4.32, -0.83, -0.11, 0.47, 1.13), 
                (-20.3, -16.7, -13.9, -11.1, -4.34, -0.84, -0.13, 0.44, 1.08), 
                (-20.5, -16.8, -14, -11.2, -4.35, -0.85, -0.14, 0.42, 1.07), 
                (-20.6, -16.9, -14.1, -11.3, -4.36, -0.85, -0.14, 0.41, 1.05))
 
    table_r3 = ((-22.5, -20, -17.9, -15.6, -8.49, -3.65, -2.51, -1.53, -0.46),
                (-25.8, -22.4, -19.7, -16.8, -8.8, -3.71, -2.6, -1.67, -0.67), 
                (-27.4, -23.7, -20.6, -17.5, -8.96, -3.74, -2.63, -1.74, -0.76), 
                (-28.5, -24.4, -21.3, -17.9, -9.05, -3.76, -2.65, -1.79, -0.83), 
                (-28.9, -24.7, -21.5, -18.1, -9.08, -3.76, -2.66, -1.8, -0.86), 
                (-29.4, -25, -21.7, -18.3, -9.11, -3.77, -2.67, -1.81, -0.88))
 
    table_t1 = ((-2.65, -2.26, -1.95, -1.6, -0.47, 0.92, 1.33, 1.7, 2.15), 
                 (-2.62, -2.25, -1.95, -1.61, -0.49, 0.91, 1.31, 1.66, 2.08), 
                 (-2.6, -2.24, -1.95, -1.61, -0.5, 0.9, 1.29, 1.64, 2.04), 
                 (-2.58, -2.24, -1.95, -1.62, -0.5, 0.89, 1.28, 1.63, 2.02), 
                 (-2.58, -2.23, -1.95, -1.62, -0.5, 0.89, 1.28, 1.62, 2.01), 
                 (-2.58, -2.23, -1.95, -1.62, -0.51, 0.89, 1.28, 1.62, 2.01))
 
    table_t2 = ((-3.75, -3.33, -2.99, -2.64, -1.53, -0.37, 0, 0.34, 0.71), 
                (-3.59, -3.23, -2.93, -2.6, -1.55, -0.41, -0.04, 0.28, 0.66),
                (-3.5, -3.17, -2.9, -2.59, -1.56, -0.42, -0.06, 0.26, 0.63), 
                (-3.45, -3.14, -2.88, -2.58, -1.56, -0.42, -0.07, 0.24, 0.62), 
                (-3.44, -3.13, -2.87, -2.57, -1.57, -0.44, -0.07, 0.24, 0.61), 
                (-3.42, -3.12, -2.86, -2.57, -1.57, -0.44, -0.08, 0.23, 0.6))
 
    table_t3 = ((-4.38, -3.95, -3.6, -3.24, -2.14, -1.14, -0.81, -0.5, -0.15),
                (-4.16, -3.8, -3.5, -3.18, -2.16, -1.19, -0.87, -0.58, -0.24), 
                (-4.05, -3.73, -3.45, -3.15, -2.17, -1.22, -0.9, -0.62, -0.28), 
                (-3.98, -3.69, -3.42, -3.13, -2.18, -1.23, -0.92, -0.64, -0.31), 
                (-3.97, -3.67, -3.42, -3.13, -2.18, -1.24, -0.93, -0.65, -0.32), 
                (-3.96, -3.67, -3.41, -3.13, -2.18, -1.25, -0.94, -0.66, -0.32))
 
    X = np.expand_dims(yt1, axis=1)
    model1 = sm.OLS(yt, X).fit()
    X = sm.add_constant(X)
    model2 = sm.OLS(yt, X).fit()
    X = np.concatenate([X, np.expand_dims(t, axis=1)], axis=1)
    model3 = sm.OLS(yt, X).fit()
 
    stats = np.array([np.array(stat(model1, 1)), np.array(stat(model2, 2)), np.array(stat(model3, 3))])
 
    pvalue_delta = [get_pval(table_r1, stats[0,0]), get_pval(table_r2, stats[1,0]),
                    get_pval(table_r3, stats[2,0])]
 
    pvalue_t_delta = [get_pval(table_t1, stats[0,1]), get_pval(table_t2, stats[1,1]),
                    get_pval(table_t3, stats[2,1])]
 
    if statistic_type == 'z_rho':
        statistic = stats[:, 0]
        p_value = pvalue_delta
    else:
        statistic = stats[:, 1]
        p_value = pvalue_t_delta
        
    return (q, statistic_type, list(statistic), p_value)

이제 샘플 데이터에 대하여 Phillips-Perron 검정을 해보자.

x = [3, 4, 4, 5, 6, 7, 6, 6, 7, 8, 9, 12, 10]

시계열을 그려보자.

 import matplotlib.pyplot as plt
 
fig = plt.figure(figsize=(7,4))
fig.set_facecolor('white')
plt.plot(range(len(x)), x, marker='o')
plt.show()

딱 봐도 정상성을 갖는 시계열은 아니다. 먼저 $Z_\tau$ 검정 통계량을 이용한 Phillips-Perron 검정을 해보자.

 results = phillips_perron_test(x, statistic_type='z_tau', lag_short=False)
print("Phillips-Perron Unit Root Test")
print("alternative: stationary")
lag = results[0]
statistic_type = results[1] 
text_dict = {
    0:"Type 1: no drift no trend",
    1:"Type 2: with drift no trend",
    2:"Type 3: with drift and trend"
}
for i, (statistic, pval) in enumerate(zip(results[2], results[3])):
    print(text_dict[i])
    print('lag : {0}, {1} : {2}, p_value : {3}'.format(lag, statistic_type, statistic, pval) )

모든 경우에 대해서 p value가 0.05보다 크다. 따라서 귀무가설은 기각되며 해당 시계열은 단위근을 갖는다고 할 수 있다.

이번엔 $Z_{\hat{\rho}}$ 를 검정 통계량으로 하는 Phillips-Perron 검정을 수행해보자.

 results = phillips_perron_test(x, statistic_type='z_rho', lag_short=False)
print("Phillips-Perron Unit Root Test")
print("alternative: stationary")
lag = results[0]
statistic_type = results[1] 
text_dict = {
    0:"Type 1: no drift no trend",
    1:"Type 2: with drift no trend",
    2:"Type 3: with drift and trend"
}
for i, (statistic, pval) in enumerate(zip(results[2], results[3])):
    print(text_dict[i])
    print('lag : {0}, {1} : {2}, p_value : {3}'.format(lag, statistic_type, statistic, pval) )

마찬가지로 (A)~(C) 모든 경우에 대하여 p value가 0.05보다 크므로 귀무가설은 기각된다. 따라서 해당 시계열은 단위근을 갖는다고 할 수 있다. 앞선 결과의 차이점은 p value가 더 보수적으로 잡혔다는 것이다.

2.2 arch 패키지 이용하기

이번엔 시계열 분석을 위한 arch 패키지를 이용하여 Phillips-Perron 검정을 수행해보자.

사용 방법은 쉽다. PhillipsPerron 클래스를 생성하고 trend 인자에 사전에 정해진 문자열을 입력한다. 여기에는 'nc'(No Trend, No Constant), 'c'(Only Constant), 'ct'(Trend and Constant)를 입력할 수 있다. 그리고 검정 통계량을 $Z_\tau$ 로 할지 $ Z_{\hat{\rho}}$로 할지는 test_type인자를 이용하여 설정한다.

먼저 $Z_\tau$ 검정 통계량을 이용한 Phillips-Perron 검정을 해보자.

 from arch.unitroot import PhillipsPerron
 
for tt in ['nc','c','ct']:
    pp = PhillipsPerron(x, trend=tt, test_type='tau')
    print(pp.summary().as_text())

결과를 보면 검정 통계량이 앞서서 직접 구현했던 것과 동일하였다. 하지만 p value가 조금 달랐는데 이는 샘플 사이즈에 따른 보간법(Interpolation) 적용 방식이 달라서 그런 것 같다. 이 경우 (A)~(B) 모든 경우에 대하여 p value가 0.05보다 크므로 귀무가설은 기각된다.

이번엔 $Z_{\hat{\rho}}$ 를 검정 통계량으로 하는Phillips-Perron 검정을 수행해보자.

 for tt in ['nc','c','ct']:
    pp = PhillipsPerron(x, trend=tt, test_type='rho')
    print(pp.summary().as_text())

이 경우에도 (A)~(B) 모든 경우에 대하여 p value가 0.05보다 크므로 귀무가설은 기각된다. 따라서 해당 시계열에 단위근이 존재한다고 할 수 있다.

참고자료

Peter C, B. Phillips, Pierre Perron - Testing for a Unit Root in Time Series Regression

https://faculty.washington.edu/ezivot/econ584/notes/unitroot.pdf

https://en.wikipedia.org/wiki/Phillips%E2%80%93Perron_test

저작자표시 비영리 변경금지

'통계 > 시계열 모형' 카테고리의 다른 글

[Change Point Detection] 1. CUSUM(Cumulative Sum) 알고리즘에 대해서 알아보자 with Python (3)	2023.01.02
[시계열 분석] 9. (Augmented) Dickey-Fuller Test(검정) with Python (419)	2022.04.14
[시계열 분석] 8. Random Walk와 단위근 검정(Unit Root Test)에 대해서 알아보자. (839)	2021.09.24
[시계열 분석] 7. ARMA 모형에 대해서 알아보자 with Python (803)	2021.09.18
[시계열 분석] 6. 이동 평균 모형(Moving Average Model) 적합하기 with Python (809)	2021.08.20

내 블로그 - 관리자 홈 전환	`Q` `Q`
새 글 쓰기	`W` `W`

글 수정 (권한 있는 경우)	`E` `E`
댓글 영역으로 이동	`C` `C`

이 페이지의 URL 복사	`S` `S`
맨 위로 이동	`T` `T`
티스토리 홈 이동	`H` `H`
단축키 안내	`Shift` + `/` `⇧` + `/`

부자 되고픈 꽁냥이

[시계열 분석] 10. Phillips-Perron Test(검정) with Python

1. Phillips-Perron Test(검정)

1.1 정의

1.2 검정 방법

1) Phillips-Perron 검정에서 검정 통계량을 계산한다.

2) 기각 여부 결정

1.3 장단점

1) 장점

2) 단점

2. Python 예제 - Phillips-Perron Test(검정)

2.1 함수 만들어보기

2.2 arch 패키지 이용하기

'통계 > 시계열 모형' 카테고리의 다른 글

댓글

티스토리툴바

단축키

내 블로그

블로그 게시글

모든 영역

	import pandas as pd
	import numpy as np
	import math
	import statsmodels.api as sm
	import warnings
	warnings.filterwarnings('ignore')

	def embed(z, lags):
	res = []
	for i in range(lags):
	temp_z = np.expand_dims(z[lags-1-i:len(z)-i], axis=1)
	res.append(temp_z)
	res = np.concatenate(res, axis=1)
	return res

	def phillips_perron_test(x, statistic_type = 'z_rho', lag_short = True):
	def stat(model, m):
	idx = 1
	if m > 1:
	idx = 2
	res1 = model.resid
	rho_hat = model.params[idx-1]
	sig_hat = model.bse[idx-1]
	s2 = np.sum(np.square(res1))/(n-m)
	gamma = []
	for i in range(1, q+2):
	u = embed(res1, i)
	gamma_val = np.sum(u[:, 0]*u[:, i-1])/n
	gamma.append(gamma_val)

	lambda2 = gamma[0]+2np.sum((1-np.array(range(1, q+1))/(q+1))gamma[1:])
	z_r = n(rho_hat-1) - ((n2)(sig_hat*2)/s2)(lambda2-gamma[0])/2
	z_t = np.sqrt(gamma[0]/lambda2)*(rho_hat-1)/sig_hat-\
	(lambda2-gamma[0])nsig_hat/(2np.sqrt(lambda2s2))
	return (z_r, z_t)

	def get_pval(table, stat):
	from scipy import interpolate
	table = np.array(table)
	ncol = table.shape[1]
	size = (25, 50, 100, 250, 500, 1000)
	percnt = (0.01, 0.025, 0.05, 0.1, 0.5, 0.9, 0.95, 0.975, 0.99)
	interpol_size = []
	for j in range(ncol):
	f_temp = interpolate.interp1d(size, table[:, j], fill_value='extrapolate')
	interpol_size.append(f_temp(n).item())

	f = interpolate.interp1d(interpol_size, percnt, fill_value='extrapolate')
	return f(stat).item()


	assert statistic_type in ['z_rho', 'z_tau']
	z = embed(x, 2)
	yt = z[:,0]
	yt1 = z[:,1]
	n = len(yt)
	assert n >= 1
	t = range(1, n+1)
	lag = 4
	if not lag_short:
	lag = 12

	q = math.ceil(lag(n/100)*(0.25))

	table_r1 = ((-11.8, -9.3, -7.3, -5.3, -0.82, 1.01, 1.41, 1.78, 2.28),
	(-12.8, -9.9, -7.7, -5.5, -0.84, 0.97, 1.34, 1.69, 2.16),
	(-13.3, -10.2, -7.9, -5.6, -0.85, 0.95, 1.31, 1.65, 2.09),
	(-13.6, -10.4, -8, -5.7, -0.86, 0.94, 1.29, 1.62, 2.05),
	(-13.7, -10.4, -8, -5.7, -0.86, 0.93, 1.29, 1.61, 2.04),
	(-13.7, -10.5, -8.1, -5.7, -0.86, 0.93, 1.28, 1.6, 2.03))

	table_r2 = ((-17.2, -14.6, -12.5, -10.2, -4.22, -0.76, 0, 0.64, 1.39),
	(-18.9, -15.7, -13.3, -10.7, -4.29, -0.81, -0.07, 0.53, 1.22),
	(-19.8, -16.3, -13.7, -11, -4.32, -0.83, -0.11, 0.47, 1.13),
	(-20.3, -16.7, -13.9, -11.1, -4.34, -0.84, -0.13, 0.44, 1.08),
	(-20.5, -16.8, -14, -11.2, -4.35, -0.85, -0.14, 0.42, 1.07),
	(-20.6, -16.9, -14.1, -11.3, -4.36, -0.85, -0.14, 0.41, 1.05))

	table_r3 = ((-22.5, -20, -17.9, -15.6, -8.49, -3.65, -2.51, -1.53, -0.46),
	(-25.8, -22.4, -19.7, -16.8, -8.8, -3.71, -2.6, -1.67, -0.67),
	(-27.4, -23.7, -20.6, -17.5, -8.96, -3.74, -2.63, -1.74, -0.76),
	(-28.5, -24.4, -21.3, -17.9, -9.05, -3.76, -2.65, -1.79, -0.83),
	(-28.9, -24.7, -21.5, -18.1, -9.08, -3.76, -2.66, -1.8, -0.86),
	(-29.4, -25, -21.7, -18.3, -9.11, -3.77, -2.67, -1.81, -0.88))

	table_t1 = ((-2.65, -2.26, -1.95, -1.6, -0.47, 0.92, 1.33, 1.7, 2.15),
	(-2.62, -2.25, -1.95, -1.61, -0.49, 0.91, 1.31, 1.66, 2.08),
	(-2.6, -2.24, -1.95, -1.61, -0.5, 0.9, 1.29, 1.64, 2.04),
	(-2.58, -2.24, -1.95, -1.62, -0.5, 0.89, 1.28, 1.63, 2.02),
	(-2.58, -2.23, -1.95, -1.62, -0.5, 0.89, 1.28, 1.62, 2.01),
	(-2.58, -2.23, -1.95, -1.62, -0.51, 0.89, 1.28, 1.62, 2.01))

	table_t2 = ((-3.75, -3.33, -2.99, -2.64, -1.53, -0.37, 0, 0.34, 0.71),
	(-3.59, -3.23, -2.93, -2.6, -1.55, -0.41, -0.04, 0.28, 0.66),
	(-3.5, -3.17, -2.9, -2.59, -1.56, -0.42, -0.06, 0.26, 0.63),
	(-3.45, -3.14, -2.88, -2.58, -1.56, -0.42, -0.07, 0.24, 0.62),
	(-3.44, -3.13, -2.87, -2.57, -1.57, -0.44, -0.07, 0.24, 0.61),
	(-3.42, -3.12, -2.86, -2.57, -1.57, -0.44, -0.08, 0.23, 0.6))

	table_t3 = ((-4.38, -3.95, -3.6, -3.24, -2.14, -1.14, -0.81, -0.5, -0.15),
	(-4.16, -3.8, -3.5, -3.18, -2.16, -1.19, -0.87, -0.58, -0.24),
	(-4.05, -3.73, -3.45, -3.15, -2.17, -1.22, -0.9, -0.62, -0.28),
	(-3.98, -3.69, -3.42, -3.13, -2.18, -1.23, -0.92, -0.64, -0.31),
	(-3.97, -3.67, -3.42, -3.13, -2.18, -1.24, -0.93, -0.65, -0.32),
	(-3.96, -3.67, -3.41, -3.13, -2.18, -1.25, -0.94, -0.66, -0.32))

	X = np.expand_dims(yt1, axis=1)
	model1 = sm.OLS(yt, X).fit()
	X = sm.add_constant(X)
	model2 = sm.OLS(yt, X).fit()
	X = np.concatenate([X, np.expand_dims(t, axis=1)], axis=1)
	model3 = sm.OLS(yt, X).fit()

	stats = np.array([np.array(stat(model1, 1)), np.array(stat(model2, 2)), np.array(stat(model3, 3))])

	pvalue_delta = [get_pval(table_r1, stats[0,0]), get_pval(table_r2, stats[1,0]),
	get_pval(table_r3, stats[2,0])]

	pvalue_t_delta = [get_pval(table_t1, stats[0,1]), get_pval(table_t2, stats[1,1]),
	get_pval(table_t3, stats[2,1])]

	if statistic_type == 'z_rho':
	statistic = stats[:, 0]
	p_value = pvalue_delta
	else:
	statistic = stats[:, 1]
	p_value = pvalue_t_delta

	return (q, statistic_type, list(statistic), p_value)

	import matplotlib.pyplot as plt

	fig = plt.figure(figsize=(7,4))
	fig.set_facecolor('white')
	plt.plot(range(len(x)), x, marker='o')
	plt.show()

	results = phillips_perron_test(x, statistic_type='z_tau', lag_short=False)
	print("Phillips-Perron Unit Root Test")
	print("alternative: stationary")
	lag = results[0]
	statistic_type = results[1]
	text_dict = {
	0:"Type 1: no drift no trend",
	1:"Type 2: with drift no trend",
	2:"Type 3: with drift and trend"
	}
	for i, (statistic, pval) in enumerate(zip(results[2], results[3])):
	print(text_dict[i])
	print('lag : {0}, {1} : {2}, p_value : {3}'.format(lag, statistic_type, statistic, pval) )

	from arch.unitroot import PhillipsPerron

	for tt in ['nc','c','ct']:
	pp = PhillipsPerron(x, trend=tt, test_type='tau')
	print(pp.summary().as_text())

[시계열 분석] 10. Phillips-Perron Test(검정) with Python

1. Phillips-Perron Test(검정)

1.1 정의

1.2 검정 방법

1) Phillips-Perron 검정에서 검정 통계량을 계산한다.

2) 기각 여부 결정

1.3 장단점

1) 장점

2) 단점

2. Python 예제 - Phillips-Perron Test(검정)

2.1 함수 만들어보기

2.2 arch 패키지 이용하기

'통계 > 시계열 모형' 카테고리의 다른 글

관련글

댓글

티스토리툴바

단축키

내 블로그

블로그 게시글

모든 영역