본 포스팅에서는 수식을 포함하고 있습니다.
티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니
PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다.
이번 포스팅에서는 오차의 2차 이상의 자기 상관 여부를 검정할 수 있는 방법을 제시한 Vinod의 논문 "Generalization of the Durbin-Watson Statistic for Higher Order Auturegressive Processes"을 리뷰하고 파이썬으로 구현해보려고 한다.
여기서 다루는 내용은 다음과 같다.
3. Generalization of the Statistic for Higher Order Autoregressive Processes
4. Generalizations of Some Variations fo the Durbin-Watson Theory
Abstract
Durbin-Watson 통계량은 선형 모형에서 잔차의 1차 자기 상관 여부를 검정할 때 사용한다. 본 논문에서는 2차 이상의 자기 상관 여부를 검정할 수 있는 방법을 제안한다.
1. Introduction and Notation
본 논문에서 제안하는 방법으로 1차 자기 상관부터 그 보다 더 높은 차수의 자기 상관 여부를 차례대로(Sequentially) 검정할 수 있다.
먼저 다음의 회귀 모형을 생각해보자.
$$Y = X\beta+u$$
여기서 $Y$는 $n\times 1$ 종속변수 벡터, $X$는 $n \times k$ 설명 변수 행렬(1열은 1로 이루어진 상수 열이다), $\beta$는 $k \times 1$ 벡터 그리고 $u$는 평균이 0이고 분산이 $\sigma^2$인 Stationary 오차 벡터이다. 이때 최소 제곱 추정량은 $\hat{\beta} = (X^tX)^{-1}X^ty$ 이다. 그리고 rank가 $n-k$인 행렬 $M = I-X(X^tX)^{-1}X$라 하자. 그리고 $\hat{u} = Y - X\hat{\beta}$라 하자.
2. Durbin-Watson Theory
Durbin, Watson은 오차에 1차 자기 상관 여부를 검정하기 위해 다음과 같은 통계량 $\hat{d}$를 제안했다.
$$\hat{d} = \sum_{t=2}^n(\hat{u}_t-\hat{u}_{t-1})^2/\sum_{t=1}^n\hat{u}_t^2\tag{2.1}$$
식 (2.1)은 다음과 같이 표현할 수 있다.
$$\hat{d} = \hat{u}^tA\hat{u}/\hat{u}^t\hat{u}\tag{2.2}$$
여기서 $A$는 대칭이고 rank가 $n-1$인 행렬이며 다음과 같이 주어진다.
$$A =
\begin{pmatrix}
1 & -1 & 0 &\cdots & 0 \\
-1 & 2 & 0 &\cdots & 0 \\
\vdots & \vdots & \vdots& \ddots & \vdots \\ 0 & 0 & \cdots & 2 & -1 \\
0 & 0 & \cdots & -1 & 1
\end{pmatrix}$$
이때 0이 아닌 $A$의 eigen value를 다음과 같이 오름차순으로 정렬하자.
$$\lambda_1 < \lambda_2 < \cdots < \lambda_{n-1} $$
행렬 $M$을 이용하면 식 (2.2)는 다음과 같다.
$$\hat{d} = \frac{u^tM^tAMu}{u^tMu}$$
$AM$의 0이 아닌 $n-k$개 eigen value $\delta_i$를 오름차순으로 정렬하자. 이 $\delta_i$는 관측 데이터인 $X$에 의존한다는 사실을 명심하자.
$u_t$가 정규분포를 따른다고 가정하자. Durbin과 Watson에 따르면 다음과 같은 $u$의 Orthogonal Linear Transformation $\xi=T(u)$ 존재한다고 한다.
$$\hat{d} = \sum_{i=1}^{n-k}\delta_i\xi_i^2/\sum_{i=1}^{n-k}\xi_i^2$$
여기서 $\xi_i$는 평균이 0이고 분산이 $\sigma^2$인 정규분포를 따르는 확률 변수이다.
Durbin과 Watson의 논문 "Testing for Serial Correlation in Least Squares Regression I"에 다음과 같은 Lemma가 있다.
$\textbf{Lemma}$
$X$의 $l$개 열이 $A$의 $l$개 eigen vector의 선형 결합으로 표현되고 나머지 $n-l$개의 eigen vector에 대응하는 eigen value $\lambda_j, \;j=1, 2, \ldots, n-l$ 오름차순으로 정렬되어 있다고 하자. 그러면 다음과 같은 부등식을 만족한다.
$$\lambda_i \leq \delta_i \leq \lambda_{i+k-l}, \;\; i=1, 2, \ldots, n-k$$
이 Lemma가 말해주는 것은 $X$의 열과 $A$의 eigen vector 간에 선형 결합 정보를 미리 알 수 있다면 $\delta_i$의 타이트한 바운드를 구할 수 있다는 것이다. 하지만 범주형 변수를 사용하는 경우가 아니라면 $l$을 미리 알 수 없으므로 웬만한 경우 $l=0$이다.
이 Lemma의 따름 정리로 다음과 같은 $d_L$, $d_U$를 생각할 수 있다.
$$d_L \leq \hat{d} \leq d_U \tag{2.3}$$
여기서 $d_L = \sum_{i=1}^{n-k}\lambda_i\xi_i^2/\sum_{i=1}^{n-k}\xi_i^2$, $d_U=\sum_{i=1}^{n-k}\lambda_{i+k-l}\xi_i^2/\sum_{i=1}^{n-k}\xi_i^2$
여기서 질문이 생길 수 있다.
그냥 $\hat{d}$의 분포를 구하면 되지 이 통계량의 하한과 상한을 왜 구하는 것일까?
앞에서 $\delta_i$는 $X$에 의존한다는 사실을 기억하는가? 위 질문에 답은 이러한 사실과 옛날 컴퓨터의 성능이 좋지 않다는 점에 있다. 즉, $\hat{d}$의 분포는 $X$에 의존하는데 주어진 $X$마다 $\hat{d}$의 분포를 안다고 하더라도 계산하기 힘들었고 그에 따라 확률값도 계산하기 힘들었다. 따라서 Durbin, Watson은 $X$에 의존하지 않는 $d_L$, $d_U$를 고려하였고 유의수준 $\alpha$, 데이터 개수 $n$ 그리고 설명변수 개수 $k$별로 그 값을 테이블로 정리하여 사람들이 사용하기 쉽게 한 것이었다.
Durbin과 Watson은 $d_L$과 $d_U$의 정확한 분포는 알아내지 못했고 베타 분포를 이용한 근사적인 분포를 사용하였고 이를 이용하여 주어진 데이터의 개수, 설명 변수의 개수에 대하여 유의 수준 $\alpha$에 대응하는 값을 계산하였다.
Durbin-Watson 통계량 $\hat{d}$는 다음과 같은 가설을 검정한다.
$H_0 : \rho_1 = 0 $ vs $H_a : \rho_1>0$
여기서 $\rho_1 = E(u_tu_{t-1})/\sigma^2$이다(1차 자기 상관 계수). 이때 $\hat{d} \leq d_L$이라면 귀무가설을 기각하고 $\hat{d} \geq d_U$라면 기각하지 않는다. 그리고 $d_L < \hat{d} < d_U$라면 판단을 보류한다. 위에서는 1차 자기 상관계수가 0인지 0보다 큰지를 검정하였는데 0보다 작은지를 검정하고 싶다면 $4-\hat{d}$를 이용한다. 이때 $4-\hat{d} \leq d_L$이면 귀무가설을 기각하고 $4-\hat{d} \geq d_U$이면 기각하지 않는다. 만약 $d_L < 4-\hat{d} < d_U$라면 판단을 보류한다. 이 경우에는 $\hat{d}$의 확률분포를 이용하면 되는데 Vinod가 논문을 쓸 당시 $\hat{d}$의 확률 분포를 정확히 계산할 수 있는 방법이 개발되었다. 이에 대한 내용은 다음 섹션에서 소개한다.
3. Generalization of the Statistic for Higher Order Autoregressive Processes
본 논문에서는 다음과 같은 가설 검정 프로세스를 차례대로 수행한다.
$j (>1)$번째 스텝에서 검정하고자 하는 가설은
$$H_0 : \rho_j = 0 \text{ when } \rho_1 = \rho_2 = \cdots =\rho_{j-1} = 0 \;\;vs\;\; H_a : \rho_j \neq 0 $$
이다.
위와 같은 검정을 하기 위해 Vinod는 다음과 같은 통계량을 제안했다.
$$\hat{d}_j = \frac{\sum_{t=j+1}^n(\hat{u}_t-\hat{u}_{t-j})^2}{\sum_{t=1}^n\hat{u}_t^2}\tag{3.1}$$
이때 $\hat{d}_1$은 (2.1)에서 $\hat{d}$와 같다.
이제 (3.1)을 다음과 같이 행렬로 표현하고자 한다.
$$\hat{d}_j = \hat{u}^tH_j\hat{u}/\hat{u}^t\hat{u}$$
$n\times n$행렬 $H_j$는 양쪽 끝 $j$번째는 1이고 나머지는 2인 대각 원소를 가지고 있으며 각 대각 원소로부터 $j$만큼 왼쪽과 아래쪽 원소는 -1인 행렬이다. 다음은 $H_2$를 나타낸 것이다.
$$H_2 =
\begin{pmatrix}
1&0&-1&0& \cdots &0&0&0&0 \\
0&1&0&-1&\cdots &0&0&0&0\\ -1&0&2&0&\cdots&0&0&0&0\\ 0&-1&0&2&\cdots&0&0&0&0 \\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots \\ 0&0&0&0&\cdots&2&0&-1&0 \\ 0&0&0&0&\cdots&0&2&0&-1 \\ 0&0&0&0&\cdots&-1&0&1&0 \\ 0&0&0&0&\cdots&0&-1&0&1
\end{pmatrix}$$
행렬 $M$을 사용하면 $\hat{d}_j$는 다음과 같이 쓸 수 있다.
$$\hat{d}_j = u^tM^tH_jMu/u^tMu\tag{3.2}$$
행렬 $H_jM$은 rank가 $m'=\min (n-j, n-k)$인 행렬이다. 이때 0이 아닌 $H_jM$의 eigen value를 오름차순으로 정리한 것을 $\delta_{j,i},\;(i=1, 2, \ldots, m')$라 하자. 그리고 $\lambda_{j,i},\;(i=1, 2, \ldots, n-k)$를 오름차순으로 정렬된 $H_j$의 eigen value라 하자. 이때 Durbin-Watson의 Lemma를 이용하면 다음과 같은 사실을 알 수 있다.
$$d_j^L \leq \hat{d}_j \leq d_j^L$$
여기서 $d_j^L = \sum_{i=1}^{n-k}\lambda_{j,i}\xi_i^2/\sum_{i=1}^{n-k}\xi_i^2$, $d_j^U=\sum_{i=1}^{n-k}\lambda_{j,i+k-l}\xi_i^2/\sum_{i=1}^{n-k}\xi_i^2$
여기서도 웬만한 경우 $l=0$로 하자. 근데 논문에서는 더해지는 항이 $n-k$개라 했는데 $m'$인 것 같지만 구현을 해보니 $n-k$가 맞는 것 같다.
Durbin, Watson이 논문을 출판한 이후 Imhof는 $\hat{d}_j$의 정확한 확률분포를 계산하는 방법을 개발하였다. 이 논문에서는 $d_j^U$은 비슷하게 구할 수 있으므로 $d_j^L$의 누적확률분포만을 고려했다. $F$를 $d_j^L$의 누적확률분포라 한다면 Imhof 방법으로 구한 $F$는 다음과 같다.
$$F(d_j^L) = \frac{1}{2}-\frac{1}{\pi}\int_0^{\infty}\frac{\sin \epsilon (u)}{u\gamma (u)}du$$
여기서 $\epsilon (u) = \frac{1}{2}\sum_{i=1}^{n-k}\tan^{-1}(b_iu)$, $\gamma (u) = \prod_{i=1}^{n-k}(1+b_i^2u^2)^{1/4}$ 그리고 $b_i = \lambda_{j,i}-d_j^L$이다.
이제 주어진 데이터 개수 $n$, 설명 변수의 개수 $k$ 그리고 차수 $j$에 대하여 유의 수준 $\alpha$에 대응하는 $d_j^L$과 $d_j^U$를 찾을 수 있다. 찾는 방법은 $F(d_j^L) = \alpha$를 풀면 된다. Vinod는 $k$는 2~6, $j$는 2~4 그리고 $n$은 15~40 그리고 40~100은 5개 단위로 해서 테이블을 만들었다.
4. Generalizations of Some Variations fo the Durbin-Watson Theory
이 섹션에서는 Durbin-Watson 통계량을 이용하여 귀무가설의 기각 여부를 판단할 수 없는 경우에 해결방안에 대해서 소개하고 있다.
Durbin-Watson이 제안한 베타 근사를 이용하여 (2.1)에 있는 $\hat{d}$ 확률 분포의 꼬리값을 계산한다. 이를 이용하여 기각 여부를 판단한다. 이외에도 몇 가지 방법을 소개한다. 하지만 Imhof의 방법으로 $\hat{d}$의 분포를 계산할 수 있으니 이 방법을 쓰자.
5. Implementation with Python
이번엔 파이썬으로 논문의 테이블에 있는 값들을 계산해보자. 먼저 필요한 모듈을 다운 받는다.
import numpy as np
from scipy.integrate import quad
from scipy.optimize import root
아래 코드는 데이터 개수 n, 절편을 제외한 설명변수의 개수 k, 차수 lag 그리고 유의 수준 alpha를 입력받아 $d_j^L$, $d_j^U$를 계산하는 코드이다. 추가적으로 정밀도를 조정하는 round_digit도 입력값에 포함했다.
def test(n,k,lag,alpha=0.05,round_digit = 10):
assert lag < n, 'invalid lag'
## make H matrix
a = np.zeros(n)
a[0] = -1
a[lag] = 1
a = np.expand_dims(a,axis=1)
for i in range(n-lag-1):
temp = np.zeros(n)
temp[i+1] = -1
temp[i+1+lag] = 1
temp = np.expand_dims(temp, axis=1)
a = np.concatenate([a,temp],axis=1)
H = a.dot(a.T)
## get eigenvalue of H in increasing order
eig_value, eig_vector = np.linalg.eig(H)
idx = eig_value.argsort()
eig_value = eig_value[idx]
lower_eig_value = eig_value[1:n-k] ## eigen values for d L
lower_eig_value = [round(x,round_digit) for x in lower_eig_value]
upper_eig_value = eig_value[k+1:] ## eigen values for d U
upper_eig_value = [round(x,round_digit) for x in upper_eig_value]
def epsilon(u,mu,x):
## u scalar, mu vector
return 0.5*np.sum(np.arctan(u*mu))
def gamma(u,mu):
## u scalar, mu vector
return np.power(np.prod(1+np.square(u*mu)),0.25)
def F(x,eig_val):
mu = eig_val-x
f = lambda u: np.sin(epsilon(u,mu,x))/(u*gamma(u,mu))
return 0.5-(1/np.pi)*quad(f,0,np.inf)[0]
def objective(x,alpha,eig_val):
return alpha-F(x,eig_val)
lower_critical_val = root(objective,1,args=(alpha,lower_eig_value)).x[0]
upper_critical_val = root(objective,2,args=(alpha,upper_eig_value)).x[0]
print('n : ',n, 'order : ',lag, ' k : ', k,
' lower : ', lower_critical_val, ' upper : ', upper_critical_val)
line 6~17
행렬 $H_j$를 계산한다.
line 19~27
$H_j$의 eigen value와 $d_j^L, d_j^U$에 대응하는 eigen value를 할당해주었다.
line 29~39
Imhof가 개발한 방법으로 $d_j^L$과 $d_j^L$의 누적 분포 함수를 계산하기 위해 필요한 함수를 정의한다. 이때 quad를 이용하여 무한대를 포함하는 구간의 적분값을 계산할 수 있다.
line 41~42
$F(d) = \alpha$를 만족하는 $d$를 찾아준다.
이제 위의 함수가 잘 작동하는지 테스트해보자. 먼저 $\alpha = 0.05, n=20$으로 고정하고 차수를 2~4, 절편항을 제외한 설명변수의 개수 k를 1~5까지 설정하고 $d_j^L, d_j^U$ 값을 출력했다.
alpha = 0.05 ## significance
n = 20 ## the number of dat
for lag in [2,3,4]:
for k in [1,2,3,4,5]:
test(n,k,lag)
아래 표는 Vinod의 논문에 나와있는 차수는 4일때 유의수준 5%에 대응하는 테이블이다. 위 그림과 아래 그림의 빨간부분을 비교해보니 $d_j^L$의 경우 -0.0002~0.0002 정도의 오차를 보였고 $d_j^U$는 정확하게 나왔다. 잘 구현된 것 같다.
이 함수는 또한 Durbin-Watson의 Table 값도 계산할 수 있다. 차수 lag를 1로 하면 된다.
## Durbin-Watson Table
for n in [15,16,17,18,19,20]:
for k in [1,2,3,4,5]:
test(n,k,lag=1)
아래 그림은 실제 테이블의 일부를 나타낸 것이다. 위 아래 그림에서 n=15인 경우를 빨간 네모로 표시했는데 비교해보니 일치했다. 그 외의 경우에도 모두 일치했다.
Durbin Watson 통계량을 1차 이상의 자기 상관 여부를 파이썬으로 검정하는 방법이 없는 것 같아서 내가 직접 개발해보고자 시도하는 과정에서 이 논문을 읽게 되었다. 글씨체가 좀 옛날 거라서 읽기 힘들었지만 그래도 이해하기는 쉬웠다. 이번 논문 리뷰를 통해 파이썬을 이용하여 1차 이상의 차수에 대해서도 자기 상관 여부를 검정할 수 있는 Durbin Watson Test를 구현했다는 점에서 의미가 있었다.
댓글