시계열 모형(특히 최소 제곱법으로 구한 모형)에서 오차의 독립성이 만족하지 않는다면 모형 파라미터의 정확성(편의 발생)이 떨어지고 예측구간의 신뢰성 또한 보장되지 않는다. 따라서 모형을 추정한 후 오차의 독립성을 만족하는지 확인해봐야할 것이다. 독립성을 만족하지 않는다면 오차는 종속성을 갖는다고 볼 수 있다. 시계열 데이터에서 오차는 종종 자기 상관(Autocorrelation)이라는 형태로 종속성을 갖게된다. 따라서 오차의 자기 상관 여부를 검정하는 방법이 필요하다. 이번 포스팅에서는 오차의 자기 상관 존재 여부를 통계적으로 검정하는 Durbin-Watson 검정을 소개한다.
여기서 다루는 내용은 다음과 같다.
2. Generalized Durbin-Watson 검정
1. Durbin-Watson 검정
- 정의 -
먼저 데이터 $(y_t, \tilde{x}_t), \; (t=1, 2, \ldots, n)$ 가 있다고 하자. 이때 최소 제곱 추정법을 통하여 계산한 적합값을 $\hat{y}_t$라 하고 이에 해당하는 잔차를 $e_t = y_t-\hat{y}_t$라 하자. 이때 잔차 사이에 다음과 같은 관계가 있다고 해보자.
$$e_t = \phi e_{t-1} +v_t$$
여기서 $v_t$는 서로 독립이며 평균이 0이고 분산이 $\sigma^2$인 동일한 분포를 따르는 확률변수이다. $\phi\neq 0$인 경우 이러한 관계를 잔차의 1차 자기 상관이라 한다. Durbin-Watson 검정은 잔차의 1차 자기 상관 여부를 검정한다. 즉, 다음의 가설을 검정한다.
$H_0$ : $\phi \leq 0$ vs $H_a$ : $\phi > 0 $(양의 1차 자기 상관)
또는
$H_0$ : $\phi \geq 0$ vs $H_a$ : $\phi < 0$(음의 1차 자기 상관)
물론 양측 검정도 가능하다.
- 검정 방법 -
잔차의 1차 자기 상관 여부를 검정하기 위해 다음과 같은 검정 통계량 $d$을 사용한다.
$$d = \frac{\sum_{t=2}^n(e_t-e_{t-1})^2}{\sum_{t=1}^ne_t^2}\tag{1-1}$$
이때 $d$값이 작으면 1차 양의 자기 상관이 있다고 판단하며 $d$값이 크면 1차 음의 자기 상관이 있다고 판단한다. 그 이유를 알아보기 위해 $d$를 풀어써보자.
$$\begin{align}d &= \frac{\sum_{t=2}^n(e_t^2+e_{t-1}^2-2e_te_{t-1})}{\sum_{t=1}^ne_t^2} \\&= \frac{\sum_{t=2}^n(e_t^2+e_{t-1}^2)}{\sum_{t=1}^ne_t^2}-2\hat{\rho} \\ &\approx 2-2\hat{\rho} = 2(1-\hat{\rho})\end{align}$$
여기서 $\hat{\rho}$은 표본 1차 자기 상관계수이다. $-1\leq \hat{\rho} \leq 1$이므로 (대략적으로) $0\leq d\leq 4$이다. 따라서 $d$가 작아질수록(0에 가까워질수록) $\hat{\rho}$은 1에 가까워진다. 즉, 1차 양의 자기 상관을 보이게 된다. 반대로 $d$가 클수록(4에 가까워질수록) $\hat{\rho}$은 -1에 가까워진다. 즉, 1차 음의 자기 상관을 보이게 된다.
유의수준 $\alpha$가 주어졌을 때 $d$의 하한 기각값과 상한 기각값을 각각 $d_{L, \alpha}, d_{U, \alpha}$라 하자.
만약 $d < d_{L, \alpha}$라면 잔차는 1차 양의 자기 상관을 갖는다고 판단하며, $d \geq d_{U, \alpha}$이면 1차 양의 자기 상관을 갖지 않는다고 판단한다. 마지막으로 $d_{L, \alpha} < d < d_{U, \alpha}$라면 판단을 보류한다.
$d' = 4-d$라 하자. 만약 $d' < d_{L, \alpha}$라면 잔차는 1차 음의 자기 상관을 갖는다고 판단하며, $d '\geq d_{U, \alpha}$이면 1차 음의 자기 상관을 갖지 않는다고 판단한다. 마지막으로 $d_{L, \alpha} < d' < d_{U, \alpha}$라면 판단을 보류한다.
그렇다면 이제 남은 것은 $d_{L, \alpha}, d_{U, \alpha}$을 찾아야한다. 설명 변수의 개수 $k'$, 데이터 개수$n$에 따른 $d_{L, \alpha}, d_{U, \alpha}$의 값을 여기에서 찾을 수 있다.
- 생각해 볼 점-
여기까지 보았다면 다음과 같은 질문 사항이 나올 수 있다.
1. 그냥 $d$를 이용하여 $p$-value를 이용하면 되지 굳이 $d_{L, \alpha}, d_{U, \alpha}$를 사용하는 이유가 무엇일까?
$d$의 분포는 설명 변수에 의존하는데 주어진 설명 변수마다 $d$의 분포를 안다고 하더라도 그 당시 컴퓨터로는 계산하기 힘들었고(사실 Durbin Watson이 논문을 출판했을 당시에는 $d$의 정확한 분포를 알아내지 못하였다) 그에 따라 확률값도 계산하기 힘들었다. 따라서 Durbin, Watson은 설명 변수에 의존하지 않는 $d_{L, \alpha}$, $d_{U, \alpha}$를 고려하였고 유의수준 $\alpha$, 데이터 개수 $n$ 그리고 설명변수 개수 $k'$별로 그 값을 테이블로 정리하여 사람들이 사용하기 쉽게 한 것이었다.
2. 판단을 보류하게 되는 상황에서는 어떻게해야하나?
판단을 보류하게 되는 이유는 $d$의 정확한 분포를 알지못해 $p$-value를 구하지 못했기 때문이다. 하지만 $d$의 분포를 정확히 알아내었고 이에 대한 확률값을 계산하는 방법이 개발되었기 때문에 판단 보류에 대한 문제는 걱정안해도 된다. 다시말해 정확한 $p$-value를 구할 수 있다는 것이다. 그에 따라 $d_{L, \alpha}, d_{U, \alpha}$가 필요가 없어졌다. 즉, Durbin-Watson 검정 통계량 관측값을 $x$라 한다면 $P(d>x)$ 또는 $P(d<x)$를 계산하여 $\alpha$와 비교해보면 된다. 예를 들어 $P(d<x) < \alpha$인 경우 유의수준 $\alpha$에 대하여 잔차가 1차 양의 자기 상관이 존재한다고 판단한다. 이 확률값은 아래 예제에서 소개하는 함수를 이용하여 계산할 수 있다.
- 예제 -
먼저 데이터를 다운받는다.
이 데이터는 여기에서 손으로 일일이 타이핑해서 만들었다. 데이터에 대한 설명이 없어서 변수가 정확히 무엇을 의미하는지 모르지만 강우량과 온도가 수확량에 얼마나 영향을 미치는지 알아보기위한 데이터라 생각한다.
필요한 모듈을 임포트하고 데이터를 불러온다.
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.integrate import quad
from scipy.optimize import root
df = pd.read_csv('./dataset/rainfall.csv')
이제 Durbin-Watson 검정을 하기 위한 함수를 정의한다. 먼저 검정통계량을 계산하는 함수를 정의한다.
def cal_dw_stat(res,lag):
n = len(res)
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)
dw_stat = res.dot(H.dot(res))/np.sum(np.square(res-np.mean(res)))
return dw_stat
다음으로 Durbin-Watson의 기각값 테이블값을 계산해주는 함수를 정의한다. n은 데이터 개수, k는 설명변수 개수, lag는 자기 상관 차수이며 유의수준 alpha, 유효자리수 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]
return lower_critical_val, upper_critical_val
다음으로 잔차를 구하기 위해 최소 제곱 모형을 적합하고 잔차를 계산한다.
y = df['yield']
X = df[['rainfall','temp']]
X = sm.add_constant(X)
X = X.values
ols_fit = sm.OLS(y,X).fit()
res = ols_fit.resid.values
이제 Durbin-Watson 검정을 하기 위한 준비는 끝났다. 검정해보자. 나는 유의수준을 0.05로 잡았다.
n = len(df)
k = X.shape[1] - 1
lag = 1
test_statistics = cal_dw_stat(res,1)
d_lower, d_upper = test(n,k,lag,alpha=0.05,round_digit = 10)
print(test_statistics, d_lower, d_upper)
왼쪽부터 순서대로 검정통계량 $d$, $d_{L, 0.05}$, $d_{U, 0.05}$이다. 위에서 보는 바와 같이
$d = 0.725950 < 0.757978 = d_{L, 0.05}$이므로 귀무가설을 $H_0 : \phi \leq 0 $ 기각하며 잔차들 사이에는 (강한) 1차 자기 상관이 존재한다고 판단한다.
이번에는 $p$-value를 구하는 함수를 정의한다. 아래 함수는 테스트할 데이터 res, 그리고 절편항이 포함된 설명 변수 행렬 X, 차수 lag, 유효자리수 round_digit를 인자로 받는다. 이 함수가 리턴하는 것은 차수, Durbin-Watson 검정통계량, 왼쪽 꼬리확률 그리고 오른쪽 꼬리확률(2개의 $p$-value)이다. 이 함수는 test 함수와 계산과정이 흡사하므로 설명은 생략한다.
def general_dw_test(res,X,lag=1,round_digit = 10):
n = X.shape[0]
k = X.shape[1]-1
X_tX_inv = np.linalg.inv(X.T.dot(X))
M = (np.identity(n) - X.dot(X_tX_inv.dot(X.T)))
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)
dw_stat = res.dot(H.dot(res))/np.sum(np.square(res-np.mean(res)))
T = H.dot(M)
T_eig_values, _ = np.linalg.eig(T)
T_eig_values = np.sort([round(x, round_digit) for x in T_eig_values if round(x,round_digit) != 0])
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]
return (lag, dw_stat, F(dw_stat,T_eig_values), 1-F(dw_stat,T_eig_values))
$p$-value를 계산해보자.
result = general_dw_test(res,X,lag=1)
print(result)
여기에서 주목해야할 값은 3번째 4번째 값이다. 이는 각각 왼쪽 꼬리확률과 오른쪽 꼬리확률이다. 왼쪽 꼬리확률은 $P(d < 0.72595) = 0.00149$이며 이는 유의수준 0.01보다 굉장히 작은 값이다. 따라서 잔차의 1차 양의 자기 상관을 갖는다고 볼 수 있다. 오른쪽 꼬리확률은 $P(d > 0.72595) = 0.99851$이며 이는 매우 큰 값이므로 1차 음의 자기 상관은 존재하지 않는다고 판단할 수 있다.
2. Generalized Durbin-Watson 검정
Durbin-Watson 검정이 1차 자기 상관을 검정하는 것이었다면 임의의 $p$차 자기 상관 검정도 생각해볼 수 있다. 여기서는 $p$차 자기 상관 여부를 검정해보는 방법을 알아보자.
- 정의 -
잔차 $e_t, \;\;(t=1, \ldots, n)$가 주어졌다고 하자. 그리고 잔차가 다음과 같은식으로 표현된다고 하자.
$$e_t = \phi_p e_{t-p} + v_t$$
여기서 $v_t$는 서로 독립이며 평균이 0이고 분산이 $\sigma^2$인 동일한 분포를 따르는 확률변수이다.
General Durbin-Watson 검정은 주어진 $p$에 대하여 $\phi_1=\phi_2=\cdots=\phi_{p-1}=0$인 경우, $H_0$ : $\phi_p=0$ vs $H_a$ : $\phi_p \neq 0$
을 검정하는 방법이다(여기서는 양측 검정을 나타냈지만 Durbin-Watson 검정에서 살펴본 것처럼 단측 검정도 가능하다). 여기서 분명히 해야할 것은 General Durbin-Watson은 단순히 $\phi_p = 0$인지 아닌지를 보는 것이 아니라 $p$보다 작은 차수의 자기 상관은 존재하지 않는 경우에 대하여 $\phi_p = 0$인지 아닌지를 테스트한다는 것이다.
- 검정 방법 -
$\phi_1=\phi_2=\cdots=\phi_{p-1}=0$인 상황에서 $p$차 자기 상관 여부를 테스트하기 위해 다음과 같은 통계량 $d_p$를 계산한다.
$$d_p = \frac{\sum_{t=p+1}^n(e_t-e_{t-p})^2}{\sum_{t=1}^ne_t^2}$$
이때 $d_p$값이 작으면 $p$차 양의 자기 상관이 있다고 판단하며 $d_p$값이 크면 $p$차 음의 자기 상관이 있다고 판단한다.
유의수준 $\alpha$과 $d_p$의 관측값이 $x$로 주어졌다고하자.
먼저 단측 검정의 경우를 살펴보자. 대립가설이 $H_a$ : $\phi_p > 0$이고 $P(d_p<x) < \alpha$ 라면 귀무가설을 기각하며 대립가설이 $H_a$ : $\phi_p < 0$인 이고 $P(d_p>x) < \alpha$라면 귀무가설을 기각한다. 양측 검정인 경우 즉, $H_a$ : $\phi_p \neq 0$ 일때 앞의 두개의 확률값 중 적어도 하나의 확률값이 $\alpha/2$ 보다 작다면 기각한다.
이때 $P(d_p<x)$ 또는 $P(d_p>x)$ 즉, $p$-value 값을 구해야하는데 Vinod(1973)가 제안한 방법을 이용하여 구해보고자 한다. 이는 아래 예제에서 확인할 수 있으며 Vinod의 논문에 대한 내용과 그가 제안한 방법을 어떻게 구현하는지 궁금하다면 여기를 참고하기 바란다.
- 예제 -
먼저 데이터를 다운 받는다.
필요한 모듈을 임포트하고 데이터를 불러온다. 해당 csv 파일은 구분자가 공백문자이므로 read_csv에서 sep 인자에 '\s+'를 넣어준다.
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.integrate import quad
df = pd.read_csv('./dataset/expenditure_stock.csv',sep='\s+')
나는 Stock 변수가 Expend에 미치는 영향을 알아보기 위하여 Stock을 설명변수, Expend를 종속변수로 하는 회귀 모형을 적합할 것이다. 그러고나서 잔차를 계산해준다.
y = df['Expend']
X = df['Stock']
X = sm.add_constant(X)
X = X.values
ols_fit = sm.OLS(y,X).fit()
res = ols_fit.resid.values
이제 Generalized Durbin-Watson 검정을 수행하는 함수를 정의한다. 이 함수는 Durbin-Watson 검정 부분에서 이미 소개하였다.
def general_dw_test(res,X,lag=1,round_digit = 10):
n = X.shape[0]
k = X.shape[1]-1
X_tX_inv = np.linalg.inv(X.T.dot(X))
M = (np.identity(n) - X.dot(X_tX_inv.dot(X.T)))
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)
dw_stat = res.dot(H.dot(res))/np.sum(np.square(res-np.mean(res)))
T = H.dot(M)
T_eig_values, _ = np.linalg.eig(T)
T_eig_values = np.sort([round(x, round_digit) for x in T_eig_values if round(x,round_digit) != 0])
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]
return (lag, dw_stat, F(dw_stat,T_eig_values), 1-F(dw_stat,T_eig_values))
나는 1차, 2차 자기 상관 여부를 테스트해보았다.
for l in [1,2]:
result = general_dw_test(res,X,lag=l)
print(result)
위 결과는 왼쪽부터 차례대로 차수, Generalized Durbin-Watson 검정통계량, 왼쪽 꼬리확률, 오른쪽 꼬리확률을 나타낸다. 위에서 보는 바와 같이 1차, 2차 양의 자기 상관이 유의한 것을 알 수 있다. 하지만 1차 양의 자기 상관을 고려해야한다. 왜냐하면 2차 양의 자기 상관 결과는 1차 자기 상관이 없다는 가정하에서 유의하기 때문이다.
1차 자기 상관은 없지만 2차 자기 상관이 있는 데이터를 찾지 못하였다... 이 데이터를 사용한 이유는 그저 2차 자기 상관 여부를 보여주기 위한 것이었다.
잔차의 자기 상관이 존재한다고 판단되면 그 다음에는 어떻게 해야할까?
이 경우에는 오차에 자기 상관 모형(Autoregressive Model)을 적용해볼 수 있다. 이에 대해서는 다음 포스팅에서 소개하겠다.
원래는 Ljung-Box 검정과 Box-Pierce 검정도 하려고 했지만 이는 ARIMA 모형의 적합성 정도를 보는 것이기 때문에 ARIMA 모형을 공부하고 나서 소개하려고 한다.
참고자료
Vinod - Generalization of the durbin-watson statistic for higher order autoregressive processes
Durbin-Watson Test - Real Statistics Using Excel
Statistical Consulting - stats.idre.ucla.edu/sas/examples/chp/regression-analysis-by-example-by-chatterjee-hadi-and-pricechapter-8-the-problem-of-correlated-errors/
'통계 > 시계열 모형' 카테고리의 다른 글
| [시계열 분석] 6. 이동 평균 모형(Moving Average Model) 적합하기 with Python (0) | 2021.08.20 |
|---|---|
| [시계열 분석] 5. 자기 상관 함수(Autocorrelation Function : ACF)과 부분 자기 상관 함수(Partial Autocorrelation : PACF) with Python (2) | 2021.07.31 |
| [시계열 분석] 4. 자기 회귀 모형(Autoregressive Model) 적합하기 with Python (2) | 2021.04.09 |
| [시계열 분석] 2. 최소 제곱법을 이용한 시계열 분석 with Python (1) | 2021.02.24 |
| [시계열 분석] 1. 시계열 데이터와 정상 과정(Stationary Process) (2) | 2021.02.21 |
댓글