자기상관분석
관련된 내용
상관성(correlation)은 두 변수간의 관계를 나타내는 것에 반해 자기상관(autocorrelation)은 한 변수내에 시간 차이에 따른 값들 사이의 관계를 파악하는 것입니다. 다시말하면 행렬의 형태로 표현되는 자료에서 자기상관성은 하나의 열 내에 존재하는 값들 사이의 관계를 나타냅니다. 반면에 상관성은 열(column)들 사이의 관계를 의미합니다. 자기상관의 정도는 식 1과 같이 자기상관 계수(Rh)로 나타냅니다.
| \begin{align}R_h& =\frac{ \text{Autocovariance}}{\text{Variance}}\\ &=\frac{\sum^{N-h}_{t=1} (x_t-\bar{x})(x_{t+h}-\bar{x})}{\sum^N_{i=1}(x_t-\bar{x})^2}\\& n\,:\;\text{자료의 크기}\\& h\,:\;\text{시차(lag time)}\end{align} | (식 1) |
자기상관계수는 함수 pandas.Series.autocorr(lag=1)에 의해 계산됩니다. 이 함수는 Series 객체 즉, 1개의 열 또는 행으로 이루어진 1차원의 벡터에서만 적용할수 있으며 지정된 lag의 차이로 두 그룹을 분리하여 pearson 상관계수를 계산하는 것입니다. 또는 statsmodels.tsa.stattools.acf() 함수를 적용하여 보다 상세한 정보를 확인할 수 있습니다.
예 1)
kospi 지수의 일일 주가 자료중 시가(Open)을 설명변수로 하여 종가(Close)를 추정하는 회귀모델을 작성하고 오차에 대한 자기상관성을 조사합니다.
| Open | Close | |
|---|---|---|
| 0 | 2874.50 | 2944.45 |
| 1 | 2943.67 | 2990.57 |
| 2 | 2993.34 | 2968.21 |
| ⋮ | ⋮ | ⋮ |
위 자료는 다음 코드로 호출한 것입니다.
import numpy as np import pandas as pd import matplotlib.pyplot as plt import FinanceDataReader as fdr
st=pd.Timestamp(2021,1, 1)
et=pd.Timestamp(2024, 5, 10)
kos=fdr.DataReader('KS11',st, et)[["Open","Close"]]
kos.index=range(len(kos)
위 자료를 표준화 합니다.
X=kos.values[:,0].reshape(-1,1) y=kos.values[:,1].reshape(-1,1) from sklearn.preprocessing import StandardScaler #독립변수 정규화(표준화) xScaler=StandardScaler().fit(X) X_n=xScaler.transform(X) #반응변수 정규화(표준화) yScaler=StandardScaler().fit(y) y_n=yScaler.transform(y)
sklearn.liniear_model.linearRegression(fit_intercept=True) 클래스를 사용하여 회귀모델(mod)를 생성합니다.
from sklearn.linear_model import LinearRegression
mod = LinearRegression()
mod.fit(X_n, y_n)
print(f'계수:{np.around(mod.coef_,3)}, \n편차 : {np.around(mod.intercept_,3)}\nR2:{np.around(mod.score(X_n, y_n),3)}')
계수:[[0.997]], 편차 : [0.] R2:0.994
위에서 생성한 모델 mod에 대한 오차에 대해 자기상관계수를 계산합니다.
pre=mod.predict(X_n)
error=y_n-pre
mu_e=error.mean()
r=np.sum((error[:-1]-mu_e)*(error[1:]-mu_e))/(np.sum((error-mu_e)**2))
print(f'R_h: {round(r, 3)}')
R_h: -0.024
자기상관계수는 함수 pandas.Series.autocorr(lag=1)에 의해 계산됩니다. 이 함수는 Series 객체 즉, 1개의 열 또는 행으로 이루어진 1차원의 벡터에서만 적용할수 있으며 지정된 lag의 차이로 두 그룹을 분리하여 pearson 상관계수를 계산하는 것입니다. 또는 statsmodels.tsa.stattools.acf() 함수를 적용하여 보다 상세한 정보를 확인할 수 있습니다.
r_h=pd.Series(error.flatten()).autocorr(lag=1)
print(f'R_h: {round(r_h, 3)}')
R_h: -0.024
자기상관계수의 값은 [-1,1]의 범위에 존재하며 0인 경우 자기상관이 없음을 나타냅니다. 자기상관계수의 신뢰구간은 식 2와 같이 계산됩니다.
| \begin{align}P\left(-z_{1-\frac{\alpha}{2}}\lt \frac{R_h-\overline{R_h}}{\frac{\sigma}{\sqrt{N}}}\lt z_{1-\frac{\alpha}{2}}\right) & =1-\alpha \\ \Rightarrow P\left(-\frac{z_{1-\frac{\alpha}{2}}}{\sqrt{N}}\lt R_h\lt \frac{z_{1-\frac{\alpha}{2}}}{\sqrt{N}}\right) & =1-\alpha\\ \because\; \overline{R_h}=0,\quad \sigma=1&\\\alpha,\, N: \text{유의수준, 샘플수} \end{align} | (식 2) |
다음은 유의수준 0.05에서의 신뢰구간은 statsmodels.tsa.stattools.acf() 함수를 사용하여 확인할 수 있습니다.
from statsmodels.tsa.stattools import acf
autoRe=acf(error, adjusted=True, nlags=1, qstat=True, alpha=0.05)
print(f"acf: {autoRe[0].round(3)}\nconfint: {autoRe[1].round(3)}\nq-state:{autoRe[2].round(3)}\np-value:{autoRe[3].round(3)}")
acf: [ 1. -0.024] confint: [[ 1. 1. ] [-0.092 0.044]] q-state:[0.488] p-value:[0.485]
위 결과는 귀무가설(H0: R_h=0)을 검정하는 것으로 다음을 의미합니다.
- acf: 시차(nlag) 0과 1에서 자기상관계수
- confint: 시차 0과 1에서의 신뢰구간
- q-state: 통계량
- p-value: 유의확률
그러므로 시차 1의 경우 자기상관성이 없다고 하는 귀무가설을 기각할 수 없습니다.
위 검정 방법과 함께 자기상관성은 식 3으로 계산되는 통계량을 기준으로 상관성의 여부를 판단할 수 있습니다. 이 검정 방법을 Durbin_Watson test 라고 합니다.
| $$\text{DW}=\frac{\sum^T_{t=2}(\text{error}_t-\text{error}_{t-1})^2}{\sum^T_{t=2}\text{error}_t^2}$$ | (식 3) |
위 검정 통계량은 대략적으로 2(1 − Rh)과 같으며 [0, 4]범위에 존재합니다. Rh = 0일 경우 이 검정통계량은 2를 나타내므로 0에 근접할수록 양의 강한 자기상관, 4에 근접할수록 강한 음의 자기상관을 나타냅니다. 일반적으로 (1.5, 2.5) 구간의 값은 자기상관을 무시할 수 있습니다. 통계량 DW는 statsmodel패키지의 stats.durbin_watson() 함수를 적용할 수 있습니다
from statsmodels.stats.stattools import durbin_watson
dw1=durbin_watson(error)
print(f'DW 통계량: {np.around(dw1, 3)}')
DW 통계량: [2.036]
위 결과는 statsmodel.api.OLS() 클래스를 적용하여 생성한 회귀모델의 객체의 메서드 summary()를 사용하여 확인할 수 있습니다. 다음 코드는 회귀모델의 메서드(reg.summary())에 의해 생성된 객체의 속성 tables를 사용한 것입니다. reg.summary()은 3개의 표들을 반환하며 속성 tables를 적용하여 각 표를 분리하여 나타낼 수 있습니다.
from statsmodels.api import add_constant, OLS X_n0=add_constant(X_n) reg=OLS(y_n, X_n0).fit() reg.summary().tables[2]
| Omnibus: | 54.386 | Durbin-Watson: | 2.036 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 226.583 |
| Skew: | 0.020 | Prob(JB): | 6.28e-50 |
| Kurtosis: | 5.564 | Cond. No. | 1.00 |
위에서 적용한 acf() 함수는 다양한 시차에서의 자기상관계수($R_h$)를 반환합니다.
acf_v=acf(error, nlags=30) print(acf_v. round(3))
[ 1. -0.024 -0.032 -0.036 -0.039 -0.018 -0.029 0.084 0.009 0.067 -0.025 -0.065 0.009 0.003 0.042 -0.062 0.015 0.029 -0.008 0.026 -0.046 -0.033 -0.072 0.007 -0.016 0.044 -0.003 -0.053 0.027 -0.032 -0.011]
시차가 0인 경우 $R_h = 1$은 자신은 자신을 참조하는 것으로서 당연한 것입니다. 그러나 그 이후는 매우 낮은 자기상관계수들을 보입니다. 이 결과는 statsmodels.graphics.tsaplots모듈의 plot_acf(x)를 사용하여 그래프로 나타낼 수 있습니다.
from statsmodels.graphics.tsaplots import plot_acf )plot_acf(error, lags=30)
댓글
댓글 쓰기