자기상관함수와 부분자기상관함수[개선 필요]

자기상관함수와 부분자기상관함수

자기상관함수와 부분자기상관함수를 좀 더 쉽게 이해하기 위해 상관계수와 부분상관함수의 개념을 알아보고 오자.

자기상관함수

시계열 데이터 $y_t, \cdots, y_1$이 있다면 이 때 자기상관함수는 아래와 같다. 시간에 따라 변하는 확률변수로 이루어진 상관함수를 의미한다.

$$ ACF(k) =cor(y_{t+k},y_t)= {{cov(y_{t+k}, y_t)}\over{\sqrt{var(y_{t+k})var(y_t)}}}$$

직관적으로 자기상관은 이전시점으로 다음시점을 예측할 수 있는 정도를 보는 것이다.

$$ y_i=\beta_0+ \beta_1 y_{i-1}$$

library(forecast)
acfValue=stats::acf(AirPassengers)
acfValue$acf[,,1]

# 구현
acf=function(x,lags=40){
  # 분산
  var=sum((x-mean(x))^2)/length(x)
  out=c()
  for(lag in 1:lags){
    lag_x=shift(x,lag)
    out=c(out,
      # 공분산
      (sum((x-mean(x))*(lag_x-mean(x)),
      na.rm=T)/length(x))/var)
  }
  # 상관계수(x, lag_x)
  out
}
acf(AirPassengers)

부분자기상관함수

ACF는 $Y_t$와 $Y_{t-2}$가 무관함에도 $Y_t$와 $Y_{t-2}$가 $Y_t$와 $Y_{t-1}$ 사이의 유의미한 관계로 인해 상관있다는 결과가 도출될 수 있다. PACF는 $t$와 $t-k$ 사이에 다른 영향력을 배제한 후에도 남아있는 상관관계를 의미한다. 자료를 찾아보았을 때 정의는 아래와 같았다. 우선 정의는 부분상관함수와 유사한 정의를 가지지만 수식은 조금 차이를 보였다.

$$PACF(k)=cor(z_{t+k}-P_{t,k} (z_{t+k}), z_t - P_{t,k} (z_t)) , \\ for k \ge 2$$

$$ P_{t,k} (x)=\beta_0+\beta_{1} z_{t+1}+\cdots \beta_{k} z_{t+k-1}$$

실제로 코드 구현에서도 차이를 보였으므로, 이는 아래 링크들을 대신하려고 한다.

설명

아래 링크에서는 설명이 잘 되어 있고 출력결과도 동일하게 나타나서 따라서 구현해 보았으나 내 자료에서는 다른 결과를 보였다.

from rpy2.robjects.packages import importr
import rpy2.robjects as robjects
import rpy2.robjects.packages as rpackages
utils = rpackages.importr('utils')
import numpy as np
import pandas as pd
forecast=importr('forecast')
series=pd.Series(
    robjects.r['AirPassengers'])

from sklearn import linear_model
df=pd.DataFrame({'T_i':series})
df['T_(i-1)']=df['T_i'].shift(1)
df['T_(i-2)']=df['T_i'].shift(2)
df=df.drop(df.index[[0,1]])# drop NA

# T_i에 대한 잔차 계산
lm = linear_model.LinearRegression()
df_X = df[['T_(i-1)']] 
df_y = df['T_i'] 
model = lm.fit(df_X, df_y)
df['Predicted_T_i|T_(i-1)'] = \
    lm.predict(df_X)
df['Residual_T_i|T_(i-1)'] = \
    df['T_i'] - 
    df['Predicted_T_i|T_(i-1)']

# T_(i-2)에 대한 잔차 계산
lm = linear_model.LinearRegression()
df_X = df[['T_(i-1)']] 
df_y = df['T_(i-2)']
model = lm.fit(df_X, df_y)
df['Predicted_T_(i-2)|T_(i-1)'] = \
    lm.predict(df_X)
df['Residual_T_(i-2)|T_(i-1)'] = \
    df['T_(i-2)'] - 
    df['Predicted_T_(i-2)|T_(i-1)']

구현

아래 링크 또한 설명이 아주 잘 되있고 구현 또한 잘 되어있지만 내가 읽기에 충분히 어려워 보였다.

내용 중에서 중요한 부분만을 추려보았다.

편의상 $Y_t$가 평균이 0 이라고 가정하고 진행하였으므로, $E(Y_l)=0$ 이므로 $\beta_0$는 0이 된다.

2. $\hat{Y_{k+1}}=\sum^k_{j=1} \phi_{kj} Y_{k+-j}$

아래 해당식의 유도과정이 잘 이해되지는 않으나, 역행렬을 구하여 $\beta_i$들의 값을 찾는 것 같다.

$$\gamma_1 = \beta_1\gamma_0 +\cdots+\beta_{k-1} \gamma_{k-2}\\ \vdots \\ \gamma_{k-1}= \beta_1 \gamma_{k-2} + \cdots + \beta_{k-1} \gamma_0 $$

$\begin{pmatrix} & \gamma_1 \\ & \gamma_2 \\ & \vdots \\& \gamma_{k-1} \end{pmatrix} = \begin{pmatrix} & \gamma_0 & \gamma_1 & \cdots & \gamma_{k-2} \\ & \gamma_1 & \gamma_0 & \cdots & \gamma_{k-3} \\ & \vdots & \vdots & \ddots & \vdots \\ & \gamma_{k-2} & \gamma_{k-3} & \cdots & \gamma_0 \end{pmatrix} \begin{pmatrix} & \beta_1 \\ & \beta_2 \\ & \vdots \\& \beta_{k-1} \end{pmatrix}$

from rpy2.robjects.packages import importr
import rpy2.robjects as robjects
import rpy2.robjects.packages as rpackages
utils = rpackages.importr('utils')
import numpy as np
import pandas as pd
forecast=importr('forecast')
series=pd.Series(robjects.r['AirPassengers'])

def pacf(series, k):
    if k == 0:
        pacf_val = 1
    else:
        gamma_array = 
          np.array([acf(series, k) 
            for k in range(1,k+1)])
 
        gamma_matrix = []
        for i in range(k):
            temp = [0]*k
            temp[i:] = [acf(series, j) 
              for j in range(k-i)]
            gamma_matrix.append(temp)
 
        gamma_matrix = np.array(gamma_matrix)
        gamma_matrix = gamma_matrix + 
          gamma_matrix.T - np.diag(
              gamma_matrix.diagonal())
        pacf_val = 
          np.linalg.inv(
            gamma_matrix).dot(gamma_array)[-1]
    return pacf_val
    
for i in range(0,10):
    print(i,pacf(series,i))