시계열 프로세스의 sample PACF 계산

Summary

본 내용은 코세라 BAYES4 Lesson 1.2.2 내용의 일부 내용정리입니다. 강의 내용 중 Autocovarinace function, Autocorrelation function 에 대한 부분을 R 코드로 전환했습니다. 그리고 특별히 Partial Autocorrelation function(이하 PACF) 은 아래 소개한 상지대학교 강의자료를 바탕으로 회귀모형 계수 추정하는 방식을 사용합니다. (자료 [1])

예제 데이터셋은 상지대학교 IT통계학과 김승구 교수님의 <2016년도 시계열해석-상관계수와 자기상관계수> 강의자료 내용을 사용했습니다.

본 내용에는 다음을 포함합니다.

• 자기공분산함수(Autocovariance function), 자기상관함수(Autocorrelation function), 편자기상관함수(PACF) 계산 방법을 R의 내장함수외 비교해서 다양한 방식으로 알아본다.

• 특히 PACF 같은 경우, 웹 검색을 통한 자료 조사를 해보면 대부분 역행렬을 이용한 방식으로 구현되었다. 여기서는 회귀모형으로 계수를 추정하고, 잔차에 대한 상관관계를 분석하는 방식으로 구현한다.

절차

데이터셋을 생성하고, R의 내장함수 acf() 를 통해서 자기상관함수 값을 계산한다. 데이터셋은 어느 농가에서 12년 동안의 쌀 수확량을 기록한 자료이다.

ts() 함수를 사용하면 자료를 시계열 개체로 변환한다.

1) R 내장 함수 acf 함수 사용

# library(dplyr)

## 1) 데이터셋 생성
yt <- c(12, 11, 15, 14, 14, 18, 16, 19, 19, 25, 22, 27)
yt <- ts(yt, start = c(1990), frequency = 1) # 시계열 개체로 변환

# 자료 개수 지정
T <- length(yt)
T

## [1] 12

# R 의 내장함수 사용해서 ACF 값을 계산한다.
a <- acf(yt, lag.max = 5, type = 'correlation')

a

## 
## Autocorrelations of series 'yt', by lag
## 
##      0      1      2      3      4      5 
##  1.000  0.573  0.493  0.188  0.076 -0.145

2) 수작업 Autocovariance, Autocorrelation 계산

지연 1(h=1), 지연 2(h=2)에 대해서 자기공분산함수, 자기상관함수를수작업으로 구하고 값을 위에서 구한 값과 서로 비교한다.

## 수작업으로 acf 값을 구한다.
# 평균을 구한다.
y_mean <- (1/T) * sum(yt)
y_mean 

## [1] 17.66667

# 자기공분산 autocovariance 계산
gamma_hat <- function (yt, y_mean, h) {
  # 데이터셋의 처음부터 (마지막 - 지연 개수) 까지의 데이터 셋과
  # 지연 개수 위치부터 마지막 데이터셋 까지의 데이터에 대해 공분산을 구한다.
  # 이때 평균은 샘플의 전체 데이터를 대상으로 구한다.
  (1/T) * sum((yt[(h+1):T] - y_mean) * (yt[1:(T-h)] - y_mean))
} 

# 지연 h = 1 에 대한 공분산 계산
# lag 명시
h <-1 
g_h <- gamma_hat (yt, y_mean, h)
g_h

## [1] 13.21296

## 자기상관을 구해서, R의 내장 함수 acf 와 결과를 비교한다. => 서로 일치함
# 자기상관 autocorrelation 계산
# 지연 h = 1
h <-1 
rho_hat <- gamma_hat (yt, y_mean, h) / gamma_hat(yt, y_mean, 0) 
rho_hat

## [1] 0.5730924

# 또는
# sum( (yt[1:(T-1)] - y_mean) * (yt[2:T] - y_mean ) ) / sum( (yt - y_mean)^2 )


# h = 2
h <-2 
rho_hat <- gamma_hat (yt, y_mean, h) / gamma_hat(yt, y_mean, 0) 
rho_hat

## [1] 0.4931727

# 또는
# rho_hat <- sum( (yt[1:(T-2)] - y_mean) * (yt[3:T] - y_mean ) ) / sum( (yt - y_mean)^2 )
# rho_hat

3) R의 내장된 pacf 사용

R 에 내장된 함수 pacf()를 사용해서 편자기상관함수 값을 계산한다.

# 편자기함수 계산
a <- pacf(yt, lag.max = 5)

a

## 
## Partial autocorrelations of series 'yt', by lag
## 
##      1      2      3      4      5 
##  0.573  0.245 -0.263 -0.075 -0.170

4) 회귀모형의 계수를 이용해서 sample PACF 계산

회귀모형을 이용한 편자기상관계수 계산법 평균 \(\mu = 0\) 이고 분산 \(\sigma^2 = 1\) 인 확률과정 \(Y_1, Y_2,,,\) 에 대해 현시점의 관측 \(Y_t\) 를 종속변수로 하고 \(k\) 이전 시점까지의 관측 \(Y_{t-1}, Y_{t-2},,, Y_{t-k}\) 를 독립변수로 하는 중회귀모형 \[Y_t = \beta_1 Y_{t-1} + \beta_2 Y_{t-2} + ... + \beta_k Y_{t-k} + \epsilon_t\] 에서 \(Y_{t-k}\) 의 계수 \(\beta_k\) 는 k 차 편자기상관계수가 된다. (인용: 자료 [1])

자료 [1] 을 참조해서 회귀모형을 구하고, 계수의 값을 확인한다. 그리고 R의 내장함수 pacf 를 이용해서 구한 값과 비교를 한다.

# 기존 자료를 평균이 0 이고 분산이 1 인 자료로 변환한다.
zt <- scale(yt)

zt_0 <- c(zt, 0, 0) # y_t
zt_1 <- c(0, zt, 0) # y_{t-1}
zt_2 <- c(0, 0, zt) # y_{t-2}

# dplyr 패키지의 lag 함수를 이용할 경우...
# lag(zt, 1, default = 0)

# 1차 편자기상관계수: x1 의 계수를 살펴보자.
mod <- lm (y ~ ., data = data.frame(y = zt_0, x1 = zt_1))
round(mod$coefficients, 3)

## (Intercept)          x1 
##       0.000       0.573

# 2차 편자기상관계수: x2 의 계수를 살펴보자.
mod <- lm (y ~ ., data = data.frame(y = zt_0, x1 = zt_1, x2 = zt_2))
round(mod$coefficients, 3)

## (Intercept)          x1          x2 
##       0.000       0.433       0.245

5) 회귀모형으로 잔차를 계산해서 sample PACF 계산 [자료 2]

자료 [2]를 참조해서 회귀모형을 구하고 잔차의 상관관계를 구한다.

1. 회귀모형을 구하고

2. 잔차를 계산하고

3. 잔차간의 상관관계를 구한다.

## lag, h = 2 에 대한 잔차 계산

# 1) 회귀모형을 구한다.
mod1 <- lm (y ~ ., data = data.frame(y = zt_0, x1 = zt_1))
mod2 <- lm (y ~ ., data = data.frame(y = zt_2, x1 = zt_1))

# 2) 각 회귀모형의 잔차를 계산한다.
res1 <- resid(mod1)
res2 <- resid(mod2)

# 3) 잔차에 대한 편자기상관함수값을 계산한다.
cor(res1, res2)

## [1] 0.2453043

# lag, h = 3
zt_0 <- c(zt, 0, 0, 0) # y_t
zt_1 <- c(0, zt, 0, 0) # y_{t-1}
zt_2 <- c(0, 0, zt, 0) # y_{t-2}
zt_3 <- c(0, 0, 0, zt) # y_{t-3}

# 1) 회귀모형을 구한다.
mod1 <- lm (y ~ ., data = data.frame(y = zt_0, x1 = zt_1, x2 = zt_2))
mod2 <- lm (y ~ ., data = data.frame(y = zt_3, x1 = zt_2, x2 = zt_1))

# 2) 잔차를 계산한다.
res1 <- resid(mod1)
res2 <- resid(mod2)

# 3) 잔차간의 상관관계를 구한다.
cor(res1, res2)

## [1] -0.2629033

# model1과 model2의 계수 비교
round(mod1$coefficients, 3)

## (Intercept)          x1          x2 
##       0.000       0.433       0.245

round(mod2$coefficients, 3)

## (Intercept)          x1          x2 
##       0.000       0.433       0.245

# 잔차 플롯
plot(res1, res2)

References

[1] 상지대학교 > IT통계학과 > 강의자료실 > 2016시계열해석-상관계수와 자기상관계수 (*****)
[2] https://timeseriesreasoning.com/contents/partial-auto-correlation/ (*****)
[3] https://www.youtube.com/watch?v=Pyb7Zd-4Qq0 (***)
[4] https://en.wikipedia.org/wiki/Partial_autocorrelation_function
[5] https://direction-f.tistory.com/65
[6] https://zephyrus1111.tistory.com/135