[시계열 분석] 9. (Augmented) Dickey-Fuller Test(검정) with Python

---

이번 포스팅에서는 단위근 검정을 위한 대표적인 방법으로 Dickey-Fuller Test(검정)과 이를 확장한 Augmented Dickey Fuller Test(검정)에 대한 내용을 알아본다. 또한 Python(파이썬)을 이용하여 (Augmented) Dickey-Fuller Test(검정)을 어떻게 수행하는지 알아본다.

 

본 포스팅을 보기에 앞서 지난 포스팅에서 다룬 단위근 검정과 Dickey-Fuller Test에 대한 기본적인 내용을 읽으면 좋다.

 

여기서 다루는 내용은 다음과 같다.

 

1. Dickey-Fuller Test(검정)

2. Augmented Dickey-Fuller Test(검정)

3. (Augmented) Dickey-Fuller Test(검정) 장단점

4. Python 예제

---

---

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

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

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

---

반응형

---

   1. Dickey-Fuller Test(검정)

1) 정의

Dickey-Fuller Test(검정)은 단위근을 검정하는 하나의 방법이다.

 

먼저 아래와 같이 1차 자기회귀모형 $AR(1)$이 있다고 하자.

$$X_{t} = \rho X_{t-1}+W_t \tag{1.1}$$

이때 시계열 이론에 따르면 $|\rho|<1$인 경우 식(1.1)은 정상 과정(Stationary Process)이며 $\rho=1$인 경우 단위근을 갖고 식 (1.1)은 비정상 과정(Non-Stationary Process)이 된다. 

---

- 참고 -

1. $|\rho|>1$인 경우는?

시계열 이론에 따르면 $|\rho|>1$인 경우도 적절한 변환을 통해 $X_{t-1}$의 계수를 1보다 작게 만들 수 있다고 한다. 따라서 $|\rho|>1$인 경우는 생각하지 않아도 된다.

2. 단위근 존재를 파악해야 하는 이유?

단위근이 있으면 $\rho$ 추정치에 편향(Bias)이 생기며 다변량 시계열 분석 시 두 시계열에 실제로 관계가 없지만 관계가 있는 가짜 회귀 모형을 만들 수 있다고 한다.

---

Dickey-Fuller Test(검정)은 귀무가설 $H_0 : \rho = 1 \text{  vs  } H_1 : \rho < 1$을 검정한다. 만약 귀무가설을 채택한다면 식 (1.1)은 단위근을 가지며(따라서 비정상 과정(Non-Stationary Process) 이다) 기각한다면 식 (1.1)은 정상 과정(Stationary Process)이라 할 수 있다. 사실 절대값이 아닌 $\rho<1$이므로 완전한 의미에 정상 과정은 아니다.

 

식 (1.1)을 변형해보자. 양변에 $X_{t-1}$을 빼준다. 그러면 

$$\Delta X_{t} = X_t-X_{t-1} = (\rho-1)X_{t-1}+W_t = \delta X_{t-1}+W_t \tag{1.2}$$

이 된다. 식 (1.2)를 이용하면

$$\begin{align} X_t \text{ has a unit root} &\Leftrightarrow  X_t = X_{t-1}+W_t \\ &\Leftrightarrow  \Delta X_t = 0\cdot X_{t-1} + W_t \end{align}$$

이 된다. 이는 $X_t$가 단위근을 갖는다는 것은 $\delta = 0$인 것과 같고 따라서 Dickey-Fuller Test(검정)은 귀무가설 $H_0 : \delta = 0 \text{  vs  } H_1 : \delta< 0$을 검정하게 된다.

 

식 (1.2)는 $X_{t-1}$과 $\Delta X_t$의 관계를 통해 단위근을 식별하는 것이라고 생각했다. 여기서 $\delta < 0$이면 즉, $X_{t-1}$이 ($X_{t-2}$보다) 증가할 때(감소할 때) $\Delta X_t$이 감소하면(증가하면) 시계열 $X_{t-1}$은 정상 과정이라는 것과 같다.

 

Dickey-Fuller Test(검정)은 아래와 같이 3 버전이 있다.

(1) 단위근 검정

$\Delta X_{t} = \delta X_{t-1}+W_t \tag{1.3}$

(2) 절편항을 포함하는 단위근 검정

$\Delta X_{t} = a_0+\delta X_{t-1}+W_t \tag{1.4}$

(3) 절편항과 선형 추세를 포함하는 단위근 검정

$\Delta X_{t} = a_0+a_1t+\delta X_{t-1}+W_t \tag{1.5}$

 

---파이썬 패키지를 보니 2차 곡선형 추세를 포함하는 단위근 검정도 나온 것 같다. --

 

Dickey-Fuller Test(검정)은 3가지 추세에 대하여 더 정확하게 검정할 수 있도록 만들어 놓았으며 각 경우에 대하여 귀무가설이 참일 때 검정 통계량의 점근 분포가 달라진다. 즉, 각 경우에 기각값이 달라진다.

---

이때 궁금증이 생겼다. 이론적으로는 알겠으나 $\delta < 0$인 것이 정상성(Stationary)과 어떻게 직관적으로 연결될까? $\delta < 0$이라는 것은 아까도 말했지만 $X_{t-1}$이 ($X_{t-2}$보다) 증가할 때(감소할 때) $\Delta X_t$이 감소하면(증가하면) 시계열 $X_{t-1}$은 정상 과정이라는 것과 같다.

 

정상 과정(Stationary Process)의 경우 아래와 같이 결정적 추세선 주위로 왔다 갔다(?)할 것이다.

 

정상 과정(Stationary Process) - Dickey-Fuller Test(검정) 설명1

먼저 $X_{t-1}$값이 $X_{t-2}$보다 작은 경우를 살펴보자. 만약 $\delta <0$이면 $X_t$값은 추세선 쪽으로 더 가까워지거나(case 1) 멀어지더라도 조금만 멀어질 것이다(case 2).

왜냐하면 $\delta$는 $X_{t-2}$에서 $X_{t-1}$로 변할 때 $\Delta X_{t-1}$에서 $\Delta X_{t}$로 변하는 평균적인 변화량이므로 근사적으로

$$\frac{\Delta X_t - \Delta X_{t-1}}{X_{t-1}-X_{t-2}}$$ 라고 볼 수 있고 이 값이 음수인 경우는 case 1, case 2 밖에 없기 때문이다.

 

따라서 $\delta < 0$이 의미하는 바는 시계열 데이터가 추세선 근방으로 왔다 갔다 한다고 볼 수 있고 이는 곧 해당 시계열이 정상 과정(Stationary Process) 임을 암시한다.

 

이번엔 $X_{t-1}$값이 $X_{t-2}$보다 큰 경우를 살펴보자. 이때에도 $\delta <0$이면 $X_t$값은 추세선 쪽으로 더 가까워지거나(case 1) 멀어지더라도 조금만 멀어질 것이다(case 2).

앞에서 설명했던 것과 동일하게 이 경우에도 $\delta <0$이면 해당 시계열이 정상 과정(Stationary Process) 임을 암시한다.

 

따라서 $\delta < 0$이라는 것은 해당 시계열이 정상 과정임을 직관적으로 이해할 수 있다. 식 (1.1)을 식 (1.2)로 바꾸는 이유가 아마도 대립 가설($H_1 : \delta < 0$) 을 바꿔줌으로써 시계열의 정상성을 좀 더 직관적으로 이해할 수 있게 함이었다는 생각이 들었다.

---

아직 끝나지 않았다...

 

데이터를 살펴보고 모델링을 하겠지만 안타깝게도 우리는 모형 (1.3)~(1.5) 중에서 어떤 모형이 진짜인지 알 수 없다. 따라서 $\delta$ 뿐만 아니라 $a_0, a_1$에 대해서도 검정을 해야 한다. Dickey, Fuller 1981년 논문 "Likelihood Ratio Statistics for Autoregressive Time Series with a Unit Root"에서 선형 추세 계수 $a_0, a_1$과 $\delta$를 동시에 검정할 수 있는 통계량과 그에 대한 기각값을 시뮬레이션을 통해 테이블로 제시했다.

 

절편이 없는 단위근 검정(1.3)의 경우는 하나만 검정하면 된다.

$$(\tau 1) \text{ } H_0 : \delta = 0 \text{ vs } H_1 : \delta < 0 \tag{1.6}$$

 

다음으로 절편이 있는 단위근 검정(1.4)은 아래와 같이 귀무가설 두 개를 검정한다. 

$$(\phi 1) \text{ } H_0 : \Delta X_t = W_t \text{ (Reduced Model) } \\ \text{ vs } H_1 : \Delta X_t = a_0+ \delta X_{t-1} + W_t \text{ (Full Model) }$$

$$(\tau 2) \text{ } H_0 : \delta = 0 \text{ vs } H_1 : \delta < 0 \tag{1.7}$$

 

마지막으로 절편과 선형 추세가 있는 단위근 검정(1.5)은 세 개의 귀무가설을 검정한다. 

$$(\phi 2) \text{ } H_0 : \Delta X_t = W_t \text{ (Reduced Model) } \\ \text{ vs } H_1 : \Delta X_t = a_0+a_1t+ \delta X_{t-1} + W_t \text{ (Full Model) }$$

$$(\phi 3) \text{ } H_0 : \Delta X_t = a_0 + W_t \text{ (Reduced Model) } \\ \text{ vs } H_1 : \Delta X_t = a_0+a_1t + \delta X_{t-1} + W_t \text{ (Full Model) }$$

$$(\tau 3) \text{ } H_0 : \delta = 0 \text{ vs } H_1 : \delta < 0 \tag{1.8}$$

2) 검정 통계량

검정 통계량은 계산하기 위해 식 (1.3)~(1.5)을 모델로 하는 선형 회귀 모형을 최소 제곱법으로 추정한다.

예를 들어 식 (1.3)의 경우 $\Delta X_t$를 반응 변수, $X_{t-1}$을 설명 변수로 하는 절편 없는 선형 모형을 최소 제곱법으로 추정한다.  식 (1.5)의 경우 $\Delta X_t$를 반응 변수, 시간 인덱스 $t$와 $X_{t-1}$을 설명 변수로 하는 절편 있는 선형 모형을 최소 제곱법으로 추정한다. 

 

이때 $X_{t-1}$에 대응하는 회귀 계수 추정치를 $\hat{\delta}$라 할 때 $\tau 1, \tau 2, \tau 3$의 검정 통계량 $T$은 다음과 같다.

$$T=\frac{\hat{\delta}}{\text{std}(\hat{\delta})}$$

이때 $\text{std}(\hat{\delta})$는 $\hat{\delta}$의 표준편차 추정치이다. 

그리고 $\phi 1, \phi 2, \phi 3$ 검정에 대한 검정 통계량 $F$는 아래와 같이 계산한다.

$$F = \frac{SSE(\text{ RM }) - SSE(\text{ FM })}{\text{DF of Residual for RM }-\text{DF of Residual for FM }} \div \frac{SSE(\text{ FM })}{\text{DF of Residual for FM }}$$

여기서 $SSE(\cdot)$은 모델에 대한 잔차 제곱합이며 RM은 Reduced Model, FM은 Full Model이며 DF of Residual for RM와 DF of Residual for FM은 각각 Reduced Model, Full Model에 대한 오차 자유도이며 이는 데이터 개수에서 추정하려는 회귀 계수를 빼면 얻을 수 있다.

3) 검정 방법

통계량 $T$와 $F$이 계산되면 다음으로 Dickey가 제시한 기각값 테이블을 비교해야 한다($T$와 $F$에 대한 테이블이 따로 있다). 이때 각 $T$ 통계량과 $F$ 통계량에 대한 기각값을 $T^*, F^*$라 하면

$T < T^*$인 경우 기각,  $F > F^*$인 경우 기각한다.

 

이때 $T$ 통계량에 대한 테이블은 다음과 같다.

dickey_fuller_tau_table.xlsx

0.01MB

(Augmented) Dickey Fuller Test(검정) 기각값 테이블1

다음은 $F$ 통계량에 대한 테이블이다. 

dickey_fuller_phi_table.xlsx

0.01MB

(Augmented) Dickey Fuller Test(검정) 기각값 테이블2

---

   2. Augmented Dickey-Fuller Test(검정)

1) 정의

Dickey-Fuller Test(검정)은 $AR(1)$ 모형에서 단위근 검정이라고 한다면 Augmented Dickey-Fuller Test(검정)은 $AR(p)$ 모형에서의 단위근 검정이다.

 

먼저 $AR(p)$ 모형을 생각해보자.

$$X_t = \phi_1X_{t-1} + \cdots + \phi_pX_{t-p} + W_t \tag{2.1}$$

여기서 $W_t$는 백색 잡음이다.

단위근을 의미하는 귀무가설을 만들고 싶다. 앞에서 Dickey-Fuller Test(검정)처럼 ($H_0 : \rho = 1 \text{ vs } H_1 : \rho < 1$ 임을 상기하자) 간단하게 만들고 싶은데 지금 상태로는 간단하게 정리가 안된다.

 

식 (2.1)을 변형해보자. Back Shift 연산자 $B$를 도입해보자. 이 말인즉슨, $B^pZ_t = Z_{t-p}$이다.

 

식 (2.1)은 다음과 같이 쓸 수 있다.

$$(1-\phi_1B-\cdots-\phi_pB^p) = W_t \tag{2.2}$$

이때 $\rho = \phi_1+\phi_2+\cdots + \phi_p$,

$\beta_j = -(\phi_{j+1}+\phi_{j+2}+\cdots + \phi_p), j=1, 2, \ldots, p-1$

라고 하면 식 (2.2)는 다음과 같다.

$$[(1-\rho B)-(\beta_1 B + \beta_2 B^2+\cdots+\beta_{p-1})(1-B)]X_t = W_t \\ \Leftrightarrow X_t = \rho X_{t-1}+\beta_1 \Delta X_{t-1}+\beta_2 \Delta X_{t-2}+\cdots +\beta_{p-1} \Delta X_{t-p+1}+W_t \tag{2.3}$$

이제 식 (2.1)을 (2.3)의 우변으로 대체한다.

$$X_t = \rho X_{t-1}+\beta_1 \Delta X_{t-1}+\beta_2 \Delta X_{t-2}+\cdots +\beta_{p-1} \Delta X_{t-p+1}+W_t \tag{2.4}$$

이때 식 (2.4)에서 $\rho=1$면 단위근을 갖게 된다. 따라서 단위근을 의미하는 귀무가설과 대립 가설을 쉽게 정의할 수 있다.

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

이때 $\delta = \rho-1$이라 하고 식 (2.4)에서 양변에 $X_{t-1}$을 빼주면 다음과 같다.

$$\Delta X_t = \delta X_{t-1}+\beta_1 \Delta X_{t-1}+\beta_2 \Delta X_{t-2}+\cdots +\beta_{p-1} \Delta X_{t-p+1}+W_t \tag{2.5}$$

그러면 Augmented Dickey-Fuller Test(검정)은 다음의 가설을 검정하게 된다.

$$H_0 : \delta = 0 \text{ vs } H_1: \delta < 0$$

 

Augmented Dickey-Fuller Test(검정) 또한 세 가지 케이스가 있다.

(1) 단위근 검정

$\Delta X_{t} = \delta X_{t-1}+ \sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} + W_t \tag{2.6}$

(2) 절편항을 포함하는 단위근 검정

$\Delta X_{t} = a_0+\delta X_{t-1}+ \sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} +W_t \tag{2.7}$

(3) 절편항과 선형 추세를 포함하는 단위근 검정

$\Delta X_{t} = a_0+a_1t+\delta X_{t-1}+\sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} +W_t \tag{2.8}$

 

절편이 없는 단위근 검정(2.6)의 경우는 하나만 검정하면 된다.

$$(\tau 1) \text{ } H_0 : \delta = 0 \text{ vs } H_1 : \delta < 0 \tag{2.9}$$

 

다음으로 절편이 있는 단위근 검정(2.7)은 아래와 같이 귀무가설 두 개를 검정한다. 

$$(\phi 1) \text{ } H_0 : \Delta X_t = \sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} + W_t \text{ (Reduced Model) } \\ \text{ vs } H_1 : \Delta X_t = a_0+ \delta X_{t-1} +\sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} + W_t \text{ (Full Model) }$$

$$(\tau 2) \text{ } H_0 : \delta = 0 \text{ vs } H_1 : \delta < 0 \tag{2.10}$$

 

마지막으로 절편과 선형 추세가 있는 단위근 검정(2.8)은 세 개의 귀무가설을 검정한다. 

$$(\phi 2) \text{ } H_0 : \Delta X_t = \sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} +W_t \text{ (Reduced Model) } \\ \text{ vs } H_1 : \Delta X_t = a_0+a_1t+ \delta X_{t-1} +\sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} + W_t \text{ (Full Model) }$$

$$(\phi 3) \text{ } H_0 : \Delta X_t = a_0 +\sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} + W_t \text{ (Reduced Model) } \\ \text{ vs } H_1 : \Delta X_t = a_0+a_1t + \delta X_{t-1} +\sum_{j=1}^{p-1}\beta_j\Delta X_{t-j} + W_t \text{ (Full Model) }$$

$$(\tau 3) \text{ } H_0 : \delta = 0 \text{ vs } H_1 : \delta < 0 \tag{2.11}$$

2) 차수 $p$의 결정

부분 자기 상관 함수(Partial Autocorrelation Function)를 이용할 수도 있으나 여기에서는 AIC(Akaike Information Criterion) 또는 BIC(Bayesian Information Criterion)을 최소화하는 차수를 결정한다.

3) 검정 통계량

이 부분은 앞서 살펴본 Dickey-Fuller Test(검정)과 매우 유사하다. 먼저 검정 통계량은 계산하기 위해 식 (2.6)~(2.8)을 모델로 하는 선형 회귀 모형을 최소 제곱법으로 추정한다.

예를 들어 식 (2.6)의 경우 $\Delta X_t$를 반응 변수, $X_{t-1}, \Delta X_{t-j}, j=1, \ldots, p-1$을 설명 변수로 하는 절편 없는 선형 모형을 최소 제곱법으로 추정한다.  식 (2.8)의 경우 $\Delta X_t$를 반응 변수, 시간 인덱스 $t$와 $X_{t-1}, \Delta X_{t-j}, j=1, \ldots, p-1$을 설명 변수로 하는 절편 있는 선형 모형을 최소 제곱법으로 추정한다. 

 

이때 $X_{t-1}$에 대응하는 회귀 계수 추정치를 $\hat{\delta}$라 할 때 $\tau 1, \tau 2, \tau 3$의 검정 통계량 $T$은 다음과 같다.
$$T=\frac{\hat{\delta}}{\text{std}(\hat{\delta})}$$
이때 $\text{std}(\hat{\delta})$는 $\hat{\delta}$의 표준편차 추정치이다. 
그리고 $\phi 1, \phi 2, \phi 3$ 검정에 대한 검정 통계량 $F$는 아래와 같이 계산한다.
$$F = \frac{SSE(\text{ RM }) - SSE(\text{ FM })}{\text{DF of Residual for RM }-\text{DF of Residual for FM }} \div \frac{SSE(\text{ FM })}{\text{DF of Residual for FM }}$$
여기서 $SSE(\cdot)$은 모델에 대한 잔차 제곱합이며 RM은 Reduced Model, FM은 Full Model이며 DF of Residual for RM와 DF of Residual for FM은 각각 Reduced Model, Full Model에 대한 오차 자유도이며 이는 데이터 개수에서 추정하려는 회귀 계수를 빼면 얻을 수 있다.

4) 검정 방법

이는 앞서 Dickey-Fuller Test(검정)에서 살펴본 것과 같이 테이블을 이용하여 기각 여부를 결정한다.

반응형

---

   3. (Augmented) Dickey-Fuller Test(검정) 장단점

- 장점 -

1. 직관적이다.

귀무가설과 대립 가설의 의미를 직관적으로 이해할 수 있다.

2. 구현이 쉽다.

(Augmented) Dickey-Fuller Test(검정)에서 선형 회귀 모형을 이용한 회귀 계수와 잔차 제곱합 등을 이용하여 검정 통계량이 만들어지고 기각 테이블도 있으므로 구현이 쉽다.

- 단점 -

1. 불완전한 대립 가설

사실 $AR(1)$의 정상성 조건은 $|\rho| < 1$인데 대립 가설이 $\rho < 1$이므로 대립 가설이 완벽하진 않다. 즉, 검정 통계량이 매우 작다면 주어진 시계열이 정상성을 갖는다고 확실히 말할 수 있는지 의문이 든다. 이 부분은 따로 찾아봐야 할 것 같다. 

2. 약한 검정력(Low Power)

(Augmented) Dickey-Fuller Test(검정)의 검정력은 약하다고 알려져 있으며 특히 실제 $\rho$가 1보다 작지만 가까운 경우(예 $\rho = 0.95$) 이를 잘 식별하지 못한다고 한다. 다시 말해 단위근을 갖지 않는 경우이지만 단위근을 갖는다고 결론을 내릴 가능성이 많다고 한다. 하지만 이는 대부분의 단위근 검정들이 갖고 있는 단점이기도 하다.

3. 데이터 개수(Sample Size)가 적은 경우 검정 정확성 문제

실제로 귀무가설이 참일 때 검정 통계량의 분포는 점근 분포(Asymptotic Distribution)라서 샘플 사이즈가 작은 경우 검정 통계량의 소표본 분포가 점근 분포와 매우 달라질 수 있어서 검정의 정확도가 떨어질 수 있다. 다시 말해 샘플 사이즈가 작은 경우 귀무가설이 참일 때 기각하는 제1종 오류를 범하는 비율이 테이블에서 정한 신뢰 수준과 매우 큰 차이가 날 수 있다.

 

---

   4. Python 예제

1) 함수

Augmented Dickey-Fuller Test(검정)을 수행하는 파이썬 함수를 만들어보자. 이 함수는 R에서 "urca" 패키지의 ur.df 함수를 참고했다.

 

먼저 $\Delta X_{t-j}, j=1, \ldots, p-1$ 변수들을 하나의 행렬(Matrix)로 만들어주는 embed 함수와 $F$ 통계량을 계산하는 get_f_stat 함수를 만들었다. get_f_stat 함수는 statsmodels에서 제공하는 회귀 분석 결과 2개(reduced_model, full_model)를 받는다.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
 
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 get_f_stat(reduced_model, full_model):
    df_num = reduced_model.df_resid - full_model.df_resid
    df_denom = full_model.df_resid
    numerator = reduced_model.mse_resid*reduced_model.df_resid - full_model.mse_resid*full_model.df_resid
    denominator = full_model.mse_resid*full_model.df_resid
    f_stat = (numerator/df_num)/(denominator/df_denom)
    return f_stat

 

다음은 실제로 Augmented Dickey-Fuller Test(검정)을 수행하는 함수이다. 시계열 데이터 y, 차수는 lags, 검정 케이스는 test_type으로 지정한다. test_type=None이면 그냥 단위근 검정이고 drift는 절편 포함, trend는 절편에 선형 추세를 포함하는 모형의 단위근 검정을 수행한다. 마지막으로 최적 차수를 계산하고 싶다면 selectlags에 'AIC' 또는 'BIC'를 지정하며 None인 경우는 lags 지정 값을 차수 $p$로 고정시킨다.

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 get_f_stat(reduced_model, full_model):
    df_num = reduced_model.df_resid - full_model.df_resid
    df_denom = full_model.df_resid
    numerator = reduced_model.mse_resid*reduced_model.df_resid - full_model.mse_resid*full_model.df_resid
    denominator = full_model.mse_resid*full_model.df_resid
    f_stat = (numerator/df_num)/(denominator/df_denom)
    return f_stat
 
def augmented_df_test(y, lags, test_type=None, selectlags=None):
    '''인자 검사'''
    if test_type:
        assert test_type in ('drift', 'trend'), 'test_type must be the one of (None, drift, trend)'
    if selectlags:
        assert selectlags in ('AIC', 'BIC'), 'test_type must be the one of (None, AIC, BIC)'
    lags += 1
    z = np.diff(y)
    n = len(z)
    x = embed(z, lags)
    z_diff = x[:,0]
    z_lag_1 = y[lags-1:n]
    tt = list(range(lags, n+1))
    '''
    1. lags가 1보다 큰 경우 selectlags가 None이 아니면 AIC 또는 BIC 기준에 따라 
    차수 p를 구한다. 
    2. test_type 별로 검정을 수행한다.
    '''
    if lags > 1:
        if selectlags:
            crit_res = []
            for i in range(2, lags+1):
                z_diff_lag = x[:, 1:i]
                if test_type is None:
                    X = np.concatenate([np.expand_dims(z_lag_1,axis=1), z_diff_lag], axis=1)
                    model = sm.OLS(z_diff, X).fit()
                elif test_type == 'drift':
                    X = np.concatenate([np.expand_dims(z_lag_1,axis=1), z_diff_lag], axis=1)
                    X = sm.add_constant(X)
                    model = sm.OLS(z_diff, X).fit()
                else:
                    X = np.concatenate([np.expand_dims(z_lag_1,axis=1), np.expand_dims(tt, axis=1)], axis=1)
                    X = np.concatenate([X, z_diff_lag], axis=1)
                    X = sm.add_constant(X)
                    model = sm.OLS(z_diff, X).fit()
 
                if selectlags == 'AIC':
                    crit_res.append(model.aic)
                else:
                    crit_res.append(model.bic)
            lags = 2+np.argmin(crit_res)
        z_diff_lag = x[:, 1:lags]
        if test_type is None:
            X = np.concatenate([np.expand_dims(z_lag_1,axis=1), z_diff_lag], axis=1)
            model = sm.OLS(z_diff, X).fit()
            tau = model.tvalues[0]
            test_stat = {'tau1':tau}
        elif test_type == 'drift':
            X = np.concatenate([np.expand_dims(z_lag_1,axis=1), z_diff_lag], axis=1)
            X = sm.add_constant(X)
            model = sm.OLS(z_diff, X).fit()
            tau = model.tvalues[1]
            model_phi = sm.OLS(z_diff, z_diff_lag).fit()
 
            phi1 = get_f_stat(model_phi, model)
            test_stat = {'tau2':tau, 'phi1':phi1}
        else:
            X = np.concatenate([np.expand_dims(z_lag_1,axis=1), np.expand_dims(tt, axis=1)], axis=1)
            X = np.concatenate([X, z_diff_lag], axis=1)
            X = sm.add_constant(X)
            model = sm.OLS(z_diff, X).fit()
            tau = model.tvalues[1]
            model_phi2 = sm.OLS(z_diff, z_diff_lag).fit()
            X_temp = sm.add_constant(z_diff_lag)
            model_phi3 = sm.OLS(z_diff, X_temp).fit()
            phi2 = get_f_stat(model_phi2, model)
            phi3 = get_f_stat(model_phi3, model)
            test_stat = {'tau3':tau, 'phi2':phi2, 'phi3':phi3}
    else:
        if test_type is None:
            model = sm.OLS(z_diff, np.expand_dims(z_lag_1,axis=1)).fit()
            tau = model.tvalues[0]
            test_stat = {'tau1':tau}
        elif test_type == 'drift':
            X = np.expand_dims(z_lag_1,axis=1)
            X = sm.add_constant(X)
            model = sm.OLS(z_diff, X).fit()
            tau = model.tvalues[1]
            rss_null = np.sum(np.square(z_diff))
            df_num = len(z_diff) - model.df_resid
            df_denom = model.df_resid
            numorator = rss_null - model.mse_resid*model.df_resid
            denominator = model.mse_resid*model.df_resid
            phi1 = (numorator/df_num)/(denominator/df_denom)
            test_stat = {'tau2':tau, 'phi1':phi1}
        else:
            X = np.concatenate([np.expand_dims(z_lag_1,axis=1), np.expand_dims(tt, axis=1)], axis=1)
            X = sm.add_constant(X)
            model = sm.OLS(z_diff, X).fit()
            tau = model.tvalues[1]
 
            rss_null = np.sum(np.square(z_diff))
            df_num = len(z_diff) - model.df_resid
            df_denom = model.df_resid
            numorator = rss_null - model.mse_resid*model.df_resid
            denominator = model.mse_resid*model.df_resid
            phi2 = (numorator/df_num)/(denominator/df_denom)
            model_phi3 = sm.OLS(z_diff, [1]*len(z_diff)).fit()
            phi3 = get_f_stat(model_phi3, model)
            test_stat = {'tau3':tau, 'phi2':phi2, 'phi3':phi3}
 
#     res = model.resid
    '''샘플 사이즈에 맞는 테이블을 가져온다.'''
    if n < 25:
        rowselec = 0
    elif n < 50:
        rowselec = 1
    elif n < 100:
        rowselec = 2
    elif n < 250:
        rowselec = 3
    elif n < 500:
        rowselec = 4
    else:
        rowselec = 5
 
    if test_type is None:
        cval_tau1 = np.array([(-2.66, -1.95, -1.6),(-2.62, -1.95, -1.61),(-2.6, -1.95, -1.61),
             (-2.58, -1.95, -1.62),(-2.58, -1.95, -1.62),(-2.58, -1.95, -1.62)])
        cval = cval_tau1[rowselec, :]
        cval = np.expand_dims(cval, axis=0)
        test_name = ['tau1']
    elif test_type == 'drift':
        cval_tau2 = np.array([(-3.75, -3, -2.63),(-3.58, -2.93, -2.6),(-3.51, -2.89, -2.58),
                         (-3.46, -2.88, -2.57),(-3.44, -2.87, -2.57),(-3.43, -2.86, -2.57)])
        cval_phi1 = np.array([(7.88, 5.18, 4.12),(7.06, 4.86, 3.94),(6.7, 4.71, 3.86),
                          (6.52, 4.63, 3.81),(6.47, 4.61, 3.79),(6.43, 4.59, 3.78)])
        cval = np.array([cval_tau2[rowselec, :],cval_phi1[rowselec, :]])
        test_name = ['tau2', 'phi1']
    else:
        cval_tau3 = np.array([(-4.38, -3.6, -3.24),(-4.15, -3.5, -3.18),(-4.04, -3.45, -3.15),
                         (-3.99, -3.43, -3.13),(-3.98, -3.42, -3.13),(-3.96, -3.41,-3.12)])
        cval_phi2 = np.array([(8.21, 5.68, 4.67), (7.02, 5.13, 4.31),(6.5, 4.88, 4.16),
                             (6.22, 4.75, 4.07),(6.15, 4.71, 4.05),(6.09, 4.68, 4.03)])
        cval_phi3 = np.array([(10.61, 7.24, 5.91),(9.31, 6.73, 5.61),(8.73, 6.49, 5.47),
                              (8.43, 6.49, 5.47),(8.34, 6.3, 5.36),(8.27, 6.25, 5.34)])
        cval = np.array([cval_tau3[rowselec, :], cval_phi2[rowselec, :], cval_phi3[rowselec, :]])
        test_name = ['tau3', 'phi2', 'phi3']
 
    results = pd.DataFrame(cval)
    results.columns = ('1pct', '5pct', '10pct')
    results.insert(0, 'TEST_NAME', test_name)
    results.insert(1, 'TEST_STAT', [test_stat[x] for x in test_name])
 
    return (results, lags-1)

반응형

---

2) 적용

이제 실제 데이터를 가지고 Augmented Dickey-Fuller Test(검정)을 수행해보자. 먼저 아래의 데이터를 다운받는다.

denmark.csv

0.00MB

데이터의 대한 설명은 다음과 같다.

ENTRY : period Time index from 1974:Q1 until 1987:Q3.
LRM : Logarithm of real money, M2.
LRY : Logarithm of real income.
LPY : Logarithm of price deflator.
IBO : Bond rate.
IDE : Bank deposit rate.

 

이제 데이터를 불러오자.

 

df = pd.read_csv('./dataset/denmark.csv')
df.head()

 

나는 로그 실질 화폐량(LRM)에 대하여 단위근이 있는지 알아보고 싶었다. 먼저 시계열 그래프를 그려보자.

 

import matplotlib.pyplot as plt
 
df['ENTRY'] = df['ENTRY'].map(lambda x: pd.to_datetime(x, format='%Y:%m'))
fig = plt.figure(figsize=(8, 5))
fig.set_facecolor('white')
plt.scatter(df['ENTRY'], df['LRM'])
plt.show()

 

정상이 아닌 시계열

주기성을 가지면서 선형 추세가 있는 것 같기도 하다. 확실한 건 시계열이 정상은 아니다 ㅎㅎ.

 

이제 Augmented Dickey-Fuller Test(검정)을 해보자. 각 검정 케이스별로 검정을 수행하고 모형 차수 $p$는 최대 4로 설정하고 AIC를 통해 2~4 중에서 최적 차수를 선택하도록 했다.

y = df['LRM'].values
for test_type in (None, 'drift', 'trend'):
    results, lags = augmented_df_test(y, lags=4, test_type=test_type, selectlags='AIC')
    print('Optimal order : ', lags)
    print(results)

 

 

해석을 해보자.

 

먼저 절편이 없는 모형의 경우 최적 차수는 3이 나왔다. $\tau 1$의 검정 통계량은 0.956이며 이는 유의 수준 10% 기각값 -1.61보다 크다. 따라서 귀무가설이 채택되고 단위근이 존재한다고 할 수 있다.

 

절편이 있는 모형의 경우 최적 차수는 5가 나왔다. $\tau 2$의 검정 통계량은 -1.702이며 이는 유의 수준 10% 기각값 -2.58보다 크다. 따라서 귀무가설이 채택되고 단위근이 존재한다고 할 수 있다. 또한 $\phi 1$ 검정 통계량은 1.849이며 이는 유의 수준 10% 기각값 3.86보다 작은 값이다. 따라서 귀무가설은 채택되며 절편항은 유의하지 않다고 판단한다.

 

마찬가지로 절편과 선형 추세가 있는 모형의 경우 $\tau 3, \phi 2, \phi 3$ 모두 채택된다. 따라서 단위근이 존재하고 절편항과 선형 추세항은 유의하지 않다고 할 수 있다.

---

$\tau 1, \tau 2, \tau 3$에 대한 검정은 statsmodels에서 제공하는 adfuller 함수를 이용하여 Augmented Dickey-Fuller Test(검정)을 수행할 수도 있다. 아래 코드는 앞의 $\tau 1, \tau 2, \tau 3$에 대한 검정을 수행한다. adfuller 함수의 regression인자에는 'c', 'ct', 'ctt', 'nc' 값을 가질 수 있는데 각각의 의미는 다음과 같다.

• “c” : 절편항이 있는 단위근 검정
• “ct” : 절편항과 선형 추세가 있는 단위근 검정
• “ctt” : 절편항, 선형 추세 그리고 이차 추세가 있는 단위근 검정
• “n” : 아무것도 없는 단위근 검정

adfuller의 함수는 튜플 타입을 리턴하며 순서대로 검정 통계량, p-value, 최적 차수, 관측치 개수, 기각값, 최적 AIC값을 출력한다.

 

from statsmodels.tsa.stattools import adfuller
 
regression_types = ['nc','c','ct']
tau_list = ['tau 1', 'tau 2', 'tau 3']
for i, regression in enumerate(regression_types):    
    results = adfuller(y, maxlag=4, regression=regression, autolag='AIC')
    print('Optimal order : ', results[2])
    print(f'{tau_list[i]} statistics : {results[0]} ', results[4])

 

 

검정 통계량은 $\tau 1$을 제외하고 동일하며 검정 결과는 모두 동일하다. 또한 최적 차수도 동일하다. 그리고 기각값에서 차이가 났다. statsmodels에서는 샘플 사이즈 값에 따라 적절한 보간법(Interpolation)을 이용하여 기각값을 더 정확하게 계산해주는 것 같다.

 

이번 포스팅을 통해 듣기만 했던 Augmented Dickey-Fuller Test(검정)에 대해서 좀 더 자세히 알 수 있는 기회가 되었다. R 코드를 통해 검정 알고리즘을 수월하게 이해할 수 있었다. 단위근 검정에 대해선 Augmented Dickey-Fuller Test(검정) 말고 여러가지가 더 있는데 기회가 되면 공부하고 포스팅해야겠다.

---

참고자료

Augmented Dickey-Fuller (ADF) Test in R

James D. Hamilton - Time Series Analysis

Ruey S. Tsay - Analysis of Financial Time Series

Peter J. Brockwell, Richard A. Davis - Time Series : Theory and Methods

 

---

---