Summary

본 내용은 코세라 강의 PTSA 의 Lesson 4.4.1, 4.4.2 내용을 그대로 옮겨서 정리했습니다. 자세한 내용은 해당 강의를 참고하세요.

강의 주소

astsa 패키지에 있는 시계열 데이터를 이용해서, 정상성 시계열 데이터로 변환 후 AR(p) 모델 피팅하는 연습을 해본다. 1번 예제의 데이터셋은 Recruitment(number of new fish) 으로 1950부터 1987년까지 453개월 동안의 새로운 물고기 개수가 된다. 2번 예제의 데이터셋은 마찬가지로 astsa 패키지의 Johnson & Johnson 의 1960부터 1980년 까지의 주가에 대한 분기별 수익을 나타낸다.

정상성 시계열 데이터에 대한 전체적인 모델 피팅 과정의 아이디어는 다음 2개이다.

  1. 평균이 0 인 정상성 시계열로 변환한다.
  2. 행렬형식의 Yule-Walker 방정식을 사용해서 모델 계수 phi 추정

모델 피팅하는 과정을 좀더 자세히 살펴보면 다음과 같다.

  1. 평균이 0 인 정상성 시계열 데이터셋에 대해서, AR(p) 모델을 피팅하고자 한다.
  2. 먼저, R 로 시계열 데이터셋에 대한 플롯, acf, pacf 를 본다.
  3. pacf 가 지연 p 차 다음에 유의미한 값을 가지지 않는다면, p 차로 모델을 피팅한다.
  4. R 로 데이터 셋에 대한 자기상관함수값 r 을 구한다. p 차까지.
  5. 행렬 R 을 정의한다.
  6. 컬럼 벡터 b 는 r 의 전치행렬로 구한다.
  7. 추정된 계수 phi.hat 은 solve(R, b) 를 통해 구한다.
  8. 노이즈에 대한 분산 추정은 자기공분산(0) * (1 - sum(r * phi.hat)) 과 같은 형태로 구한다.

1번 예제 - AR(2) 모델 피팅

원래 모델은 \(X_t = \phi_0 + \phi_1 X_{t-1} + \phi_2 X_{t-2} + Z_t, Z_t \sim N(0, \sigma^2)\)

AR(2) 모델로 피팅해보자. 먼저, 데이터에 대한 시계열 데이터를 평균이 0인 시계열 데이터로 변환한다. 그리고 플롯, acf, pacf 를 그려본다. 아래 그림을 보면, acf 는 교대하면서 점차 감소하는 경향을 보이고 있고, pacf 는 지연 1, 2는 유의미한 값을 가지지만 지연 3 이후로는 유의미한 값을 가지지 않는다.

따라서 해당 시계열 데이터는 AR(p) 모델에서 차수 p 를 2 로 유추할 수 있다.

차수 p = 2 로 설정하고, 샘플 자기상관함수 r 을 구한다. 자기상관함수값의 지연 0 에서의 값은 항상 1이 된다. p 차이기 때문에 2개의 2개의 자기상관함수값을 구하는데 지연 1 에서의 값이 r1 이 되고 지연 2 에서의 값이 r2 가 된다.

p <- 2 # 차수는 2로 설정

# TEST
# result <- acf(ar.process, plot=F)
# summary(result)
# result$acf[2:3]

# sample autocorrelation function 의 처음 2개 값을 추출한다.
r <- acf(ar.process, plot=F)$acf[2:(p+1)] # XXX: 괄호 사용에 주의하자. : 연산자가 + 보다 우선순위가 높다.
r
## [1] 0.9218042 0.7829182

Yule-Waler 방정식: \(R \hat{\phi} = \hat{b}\)

\(\begin{bmatrix} r_0 & \cdots & r_{p-1} \\ \vdots & \ddots & \vdots \\ r_{p-1} & \cdots & r_0 \end{bmatrix} \begin{bmatrix} \hat{\phi_1} \\ \hat{\phi_2} \end{bmatrix} = \begin{bmatrix} \hat{r_1} \\ \hat{r_2} \end{bmatrix}\)

여기서 행렬 R 과 컬럼 벡터 b 를 정의한다.

먼저 행렬 R 을 정의한다. r_0 는 항상 1 이기 때문에…

\(R = \begin{bmatrix} 1 & \cdots & r_{p-1} \\ \vdots & \ddots & \vdots \\ r_{p-1} & \cdots & 1 \end{bmatrix}\)

p = 2 에서

\(R = \begin{bmatrix} 1 & r_1 \\ r_1 & 1 \end{bmatrix}\)

컬럼 벡터 b 는 r 로 구성된다.

\(b = \begin{bmatrix} \hat{r_1} \\ \hat{r_2} \end{bmatrix}\)

# matrix R
R <- matrix(1,p,p) # matrix of dimension p by p, with entries all 1's.

# define non-diagonal entires of R
for(i in 1:p){
  for(j in 1:p){
    if(i!=j)
      R[i,j] <- r[abs(i-j)]
  }
}

R
##           [,1]      [,2]
## [1,] 1.0000000 0.9218042
## [2,] 0.9218042 1.0000000
# 벡터 b 정의
b <- matrix(r, p, 1)
b
##           [,1]
## [1,] 0.9218042
## [2,] 0.7829182

다음으로 모델의 계수인 phi 값을 계산한다. 이때 R 의 solve() 함수를 사용하는데, 일반적인 연랍방정식을 푸는 함수라고 생각하면 된다.

예를 들어,

A %*% x = b 와 같은 연립방정식이 있을 때, x 의 값을 구할 때, solve(A, b) 라고 하면 x 의 값이 구해지는거야.

phi.hat <- solve(R, b)
phi.hat
##            [,1]
## [1,]  1.3315874
## [2,] -0.4445447

이제 AR(2) 모델에서 노이즈에 대한 분산을 추정해보자. 분산을 계산하기 위해서 지연 0 에서의 자기공분산 함수값 c0 가 필요하다.

\(\hat{\sigma^2} = c0 * [1 - \Sigma_i^p \hat{\phi_i} r_i]\)

c0 는 R 의 acf() 함수를 통해서 구한다.

# 샘플 자기공분산함수를 구한다.
c0 <- acf(ar.process, type='covariance', plot = F)$acf[1]
c0
## [1] 780.991
var.hat <- c0 * (1 - sum(phi.hat * r))
var.hat
## [1] 94.17131

원래 모델의 상수항 phi_0 를 추정해보자. 상수항은 다음과 같이 계산된다.

\(\hat{\phi_0} = \bar{x}(1 - \Sigma_{i=1}^p \hat{\phi_i})\)

# https://blog.naver.com/skkong89/223069860164

phi0.hat <- mean(ar.data) * (1 - sum(phi.hat))
phi0.hat
## [1] 7.033036

전체적인 모델은 다음과 같이 피팅된다.

\(X_t = 7.0330363 + 1.3315874 X_{t-1} + -0.4445447 X_{t-2} + Z_t\)

\(Z_t \sim N(0, 94.1713101)\)

실제로 시계열 데이터가 정상성인지에 대한 검정이나, 모델이 데이터를 어느 정도 설명하는지에 대한 평가 등은 본 주제의 범주를 벗어난다.

참고로 ar() 함수를 사용한 계수 추정 결과와 비교해보자.

# 상수항을 구해주지는 않는다.
ar(ar.data, order.max = 2, method = 'yule-walker')
## 
## Call:
## ar(x = ar.data, order.max = 2, method = "yule-walker")
## 
## Coefficients:
##       1        2  
##  1.3316  -0.4445  
## 
## Order selected 2  sigma^2 estimated as  94.8

2번 예제 - AR(4) 모델 피팅

AR(4) 모델은 \(X_t = \phi_0 + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \phi_3 X_{t-3} + \phi_4 X_{t-4} + Z_t, Z_t \sim N(0, \sigma^2)\)

2번째 데이터셋은 Johnson & Johnson 의 분기별 수익에 대한 내용이다. 먼저 시계열 데이터 플롯을 그려보자.

시간이 지남에 따라 평균도 증가하고 있고, 분산도 점점 커지고 있음을 볼 수 있다. 즉, 정상성 시계열이 아니다.

ar.data <- JohnsonJohnson
plot(ar.data, main='Johnson & Johnson earnings per share')

여기서는 log-returns 라는 변환과정을 사용한다.

\(r_t = log( \frac{X_t}{X_{t-1}} ) = log(X_t) - log(X_{t-1})\)

R 에서는 diff(log(데이터)) 와 같은 형태를 사용한다.

ar.process.log.return <- diff(log(ar.data)) # log-returns 변환 적용
ar.process <- ar.process.log.return - mean(ar.process.log.return) # 평균을 0 으로 이동시킴

par(mfrow = c(2,1))
plot(ar.process.log.return, main='J&J log-returns transform')
plot(ar.process, main='J&J log-returns + mean 0')

## TEST: log-returns 데이터셋을 원래 데이터셋으로 복원
# restored.ar <- exp(cumsum(ar.process.log.return))
# restored.ar <- ts(restored.ar, start = 1960, frequency = 4)
# par(mfrow = c(2,1))
# plot(ar.data, main = 'J&J original dataset')
# plot(restored.ar, main = 'restored dataset from log-returns')

acf, pacf 를 플롯해보자. acf 는 교대로 점차 감소하는 모습을 보이고 있고, pacf 는 지연 4 이후로 유의미한 값을 가지지 못한다.

모델의 차수 p 는 4로 추정한다.

par(mfrow=c(3,1))
plot(ar.process, main='Log-return (mean zero) of Johnson&Johnosn earnings per share')
acf(ar.process, main='ACF of J&J')
pacf(ar.process, main='PACF of J&J')

p = 4 로 추정했다. 다음으로 샘플 자기상관함수를 구한다.

p <- 4

r <- acf(ar.process, plot = F)$acf[2:(p+1)]

행렬 R 과 컬럼 벡터 b 를 구한다.

# matrix R
R <- matrix(1,p,p) # matrix of dimension 4 by 4, with entries all 1's.

# define non-diagonal entires of R
for(i in 1:p){
  for(j in 1:p){
    if(i!=j)
      R[i,j]=r[abs(i-j)]
  }
}
R
##             [,1]        [,2]        [,3]        [,4]
## [1,]  1.00000000 -0.50681760  0.06710084 -0.40283604
## [2,] -0.50681760  1.00000000 -0.50681760  0.06710084
## [3,]  0.06710084 -0.50681760  1.00000000 -0.50681760
## [4,] -0.40283604  0.06710084 -0.50681760  1.00000000
# b-column vector on the right
b <- matrix(r,p,1)# b- column vector with no entries

모델의 계수를 추정한다.

phi.hat <- solve(R, b)
phi.hat
##            [,1]
## [1,] -0.6293492
## [2,] -0.5171526
## [3,] -0.4883374
## [4,]  0.2651266

분산을 추정한다.

# Variance estimation using Yule-Walker Estimator
c0 <- acf(ar.process, type='covariance', plot=F)$acf[1]
c0
## [1] 0.04365692
var.hat <- c0*(1-sum(phi.hat*r))
var.hat
## [1] 0.01419242

상수항 phi0 는 다음과 같다.

# Constant term in the model
phi0.hat <- mean(ar.process)*(1-sum(phi.hat))
phi0.hat
## [1] 4.958957e-18

전체적인 모델은 다음과 같이 피팅된다. 참고로 우리는 X_t 에 대해서 모델 계수를 구한게 아니고, X_t 를 log-returns 변환한 r_t 에 대해서 피팅한거야.

\(r_t = 4.9589573\times 10^{-18} + -0.6293492 r_{t-1} + -0.5171526 r_{t-2} + -0.4883374 r_{t-3} + 0.2651266 r_{t-4}+ Z_t\)

\(Z_t \sim N(0, 0.0141924)\)

ar() 함수의 결과와 비교해보자.

ar(ar.process, order.max = 4)
## 
## Call:
## ar(x = ar.process, order.max = 4)
## 
## Coefficients:
##       1        2        3        4  
## -0.6293  -0.5172  -0.4883   0.2651  
## 
## Order selected 4  sigma^2 estimated as  0.0151