[data analysis] 자기상관분석(Autocorrelation Analysis)

자기상관분석

관련된 내용

• 회귀모델의 오차에 대해
• 오차제곱평균(Mean of Square Error)
• 오차의 분산

상관성(correlation)은 두 변수간의 관계를 나타내는 것에 반해 자기상관(autocorrelation)은 한 변수내에 시간 차이에 따른 값들 사이의 관계를 파악하는 것입니다. 다시말하면 행렬의 형태로 표현되는 자료에서 자기상관성은 하나의 열 내에 존재하는 값들 사이의 관계를 나타냅니다. 반면에 상관성은 열(column)들 사이의 관계를 의미합니다. 자기상관의 정도는 식 1과 같이 자기상관 계수(R_h)로 나타냅니다.

\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 − R_h)과 같으며 [0, 4]범위에 존재합니다. R_h = 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)