8.7 R에서 ARIMA 모델링
auto.arima()은 어떻게 작동하는가?
R에서 auto.arima() 함수는 힌드만-칸다카르(Hyndman-Khandakar) 알고리즘 (Hyndman & Khandakar, 2008)의 변형된 형태를 사용합니다. auto.arima() 함수는 ARIMA 모델을 얻기 위해, 단위 근 검정, AICc 최소화, MLE를 결합하여 사용합니다. 여기에서 설명하는 내용은 기본값으로 설정된 작동 방식입니다.
| 자동화된 ARIMA 모델링을 위한 힌드만-칸다카르(Hyndman-Khandakar) 알고리즘 |
|---|
|
|
|
|
|
|
기본적인 과정은 찾는 속도를 높이기 위해 몇 가지 근사 과정을 사용합니다. 입력값 approximation=FALSE으로 이러한 근사를 생략할 수 있습니다. 이러한 근사나 단계별 탐색 과정 때문에 최소 AICc를 갖는 모델을 찾지 못할 수도 있습니다. 입력값 stepwise=FALSE을 사용하면, 훨씬 더 큰 모델의 모음을 검색합니다. 입력값 전체를 설명하는 내용은 도움말 파일을 참조하시길 바랍니다.
여러분만의 모델을 고르기
여러분만의 모델을 고르고 싶다면, R에서 Arima() 함수를 사용하시길 바랍니다. R에서 ARIMA 모델을 맞추는 arima()라는 또 다른 함수도 있습니다. 하지만, 이 함수는 \(d=0\)이 아닐 때 상수 \(c\)를 제외하고, forecast 패키지로 작업할 수 있는 다른 함수에서 필요하는 것을 출력값으로 내지 않습니다. 마지막으로, 이 함수는 새로운 데이터에 적용하기 위한 추정된 모델을 사용할 수 없습니다(이러한 추정된 모델은 예측 정확도를 확인할 때 유용합니다). 결과적으로 대신에 Arima() 함수를 사용하길 추천합니다.
모델링 과정
ARIMA 모델을 (비-계절성) 시계열 데이터 모음에 맞출 때, 다음과 같은 과정은 유용한 일반적인 접근 방식이 될 수 있습니다.
- 데이터를 그래프로 나타내고, 특이한 관측값을 찾습니다.
- 필요하다면, 분산을 안정화하기 위해 (박스-칵스(Box-Cox) 변환을 사용하여) 데이터를 변환합니다.
- 데이터가 정상성을 나타내지 않는다면, 데이터가 정상성을 나타낼 때까지 데이터를 가지고 1차 차분을 계산합니다.
- ACF/PACF를 살펴봅니다: ARIMA(\(p,d,0\))이 적절한가? 아니면 ARIMA(\(0,d,q\)) 모델이 적절한가?
- 여러분이 고른 모델을 사용해보고, 더 나은 모델을 찾기 위해 AICc를 이용합니다.
- 잔차의 ACF를 그리고 잔차의 포트맨토(portmanteau) 검정을 하여 여러분이 고른 모델에서 잔차를 확인합니다. 백색잡음 같이 보이지 않는다면, 수정된 모델을 사용합니다.
- 잔차가 백색잡음처럼 보이면, 예측값을 계산합니다.
힌드만-칸다카르(Hyndman-Khandakar) 알고리즘은 3–5 단계만 처리합니다. 그래서 여러분이 이 알고리즘을 사용하더라도, 여전히 여러분 스스로 나머지 단계를 처리해야할 필요가 있습니다.
그림 8.11 에 과정을 요약하였습니다.
Figure 8.11: ARIMA 모델을 사용하여 예측하는 일반적인 절차.
예제: 계절성으로 조정된 전자 장비 주문량
이 과정을 그림 8.12에 나타낸 계절성으로 조정된 전자 장비 주문량 데이터에 적용할 것입니다.
elecequip %>% stl(s.window='periodic') %>% seasadj -> eeadj
autoplot(eeadj)
Figure 8.12: 계절성 조정된 유로 지역의 전자 장비 주문량 지수.
시간 그래프는 몇몇 갑작스런 변화를 나타냅니다. 특별히 2008/2009년에 크게 하락합니다. 이러한 변화는 전세계적인 경제 상황 때문에 일어났습니다. 이외에는 시간 그래프에서 특별한 점이 없고 데이터 조정이 필요한 것 같이 보이진 않습니다.
분산이 변한다는 증거가 없으니, 박스-칵스(Box-Cox) 변환을 하지 않을 것입니다.
장기적으로 볼 때 시계열이 위아래로 출렁이기 때문에 데이터가 정상성을 나타내지 않는다는 것이 분명합니다. 결과적으로, 데이터에서 1차 차분을 구할 것입니다. 8.13는 차분을 구한 데이터를 나타냅니다. 정상성을 나타내는 것으로 보이고, 차분을 더 구할 필요가 없습니다.
eeadj %>% diff %>% ggtsdisplay(main="")
Figure 8.13: 계절성 조정된 전자 장비 데이터를 차분한 것에 대한 시간 그래프와 ACF와 PACF 그래프.
그림 8.13에 나타낸 PACF는 AR(3) 모델을 암시합니다. 따라서, 초기 후보 모델은 ARIMA(3,1,0)입니다. 다른 분명한 후보 모델은 없습니다.
ARIMA(3,1,0) 모델을 ARIMA(4,1,0), ARIMA(2,1,0), ARIMA(3,1,1) 등과 함께 맞춥니다. 이 중에서 ARIMA(3,1,1)이 살짝 작은 AICc 값을 갖습니다.
(fit <- Arima(eeadj, order=c(3,1,1))) #> Series: eeadj #> ARIMA(3,1,1) #> #> Coefficients: #> ar1 ar2 ar3 ma1 #> 0.004 0.092 0.370 -0.392 #> s.e. 0.220 0.098 0.067 0.243 #> #> sigma^2 = 9.58: log likelihood = -492.7 #> AIC=995.4 AICc=995.7 BIC=1012ARIMA(3,1,1) 모델에서 얻은 잔차(residual)의 ACF 그래프는 모든 자기상관 관계값이 한도 안에 들어오는 것을 나타냅니다. 이것은 잔차가 백색잡음처럼 행동하는 것을 의미합니다. 포트맨토(portmanteau) 검정이 큰 p-값을 돌려주었습니다. 이것도 잔차가 백색잡음(white noise)이라는 것을 암시합니다.
checkresiduals(fit)
Figure 8.14: ARIMA(3,1,1) 모델에 대한 잔차 그래프.
#> #> Ljung-Box test #> #> data: Residuals from ARIMA(3,1,1) #> Q* = 24, df = 20, p-value = 0.2 #> #> Model df: 4. Total lags used: 24그림 8.15은 고른 모델에서 얻은 예측값을 나타냅니다.
autoplot(forecast(fit))
Figure 8.15: 계절성 조정된 전자 주문 지수에 대한 예측값.
기본 설정을 사용한 자동화된 알고리즘을 사용했다면, ARIMA(3,1,0) 모델을 얻을 것이지만, approximation=FALSE을 사용했다면, ARIMA(3,1,1) 모델을 얻을 것입니다.
R에서 상수를 이해하기
비-계절성 ARIMA 모델을 다음과 같이 쓸 수 있습니다. \[\begin{equation} \tag{8.4} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d y_t = c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] 또는 다음의 것과 동일하게 \[\begin{equation} \tag{8.5} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d (y_t - \mu t^d/d!) = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] 여기에서 \(c = \mu(1-\phi_1 - \cdots - \phi_p )\)이고 \(\mu\)는 \((1-B)^d y_t\)의 평균입니다. R에서는 매개변수화된 식 (8.5)을 사용합니다.
따라서, 정상성을 나타내지 않는 ARIMA 모델에서 상수가 포함되는 것은 예측 함수에서 차수 \(d\)의 다항식 추세를 유도하는 것과 같습니다. (상수가 생략되면, 예측 함수는 차수 \(d-1\)의 다항식 추세를 포함합니다.) \(d=0\)일 때는, \(\mu\)가 \(y_t\)의 평균인 특수한 상황이 됩니다.
기본값으로, Arima() 함수는 \(d>0\)일 때 \(c=\mu=0\)로 정하고 \(d=0\)일 때 \(\mu\)의 추정치를 냅니다. 이 값은 시계열의 표본 평균에 가까워질 것이지만, \(p+q>0\)일 때, 표본 평균이 최대 가능도 추정(maximum likelihood estimation)이 아니기 때문에 보통은 정확하게 같지는 않습니다.
입력값 include.mean은 \(d=0\)일 때만 영향을 주고 기본값은 TRUE입니다. include.mean=FALSE로 두면 강제로 \(\mu=c=0\)로 둘 것입니다.
입력값 include.drift은 \(d=1\)일 때 \(\mu\ne0\)을 허용합니다. \(d>1\)일 때는 2차나 이상의 차수 추세가 예측할 때 특별히 위험할 수 있기 때문에 어떠한 상수도 허용하지 않습니다. \(d=1\)일 때 매개변수 \(\mu\)는 R 출력에서 “표류(drift)”라고 부릅니다.
include.constant라는 입력값도 있는데 TRUE이면 \(d=0\)일 때 include.mean=TRUE로 둘 것이고, \(d=1\) 일 때는 include.drift=TRUE로 둡니다. include.constant=FALSE이면, include.mean와 include.drift 둘 다 FALSE로 두게 됩니다. include.constant을 사용하면, include.mean=TRUE와 include.drift=TRUE는 무시됩니다.
auto.arima() 함수는 상수를 포함하는 과정을 자동화합니다. 기본값으로, \(d=0\)나 \(d=1\)에 대해, AICc 값이 좋아지면 상수가 추가될 것입니다. \(d>1\)에 대해, 상수가 항상 생략됩니다. allowdrift=FALSE로 정하면, \(d=0\)일 때만 상수가 허용됩니다.
특성근을 그래프로 나타내기
(이 절의 내용은 좀 더 어렵습니다. 원하면 넘어가도 됩니다.)
식 (8.4) 을 다음과 같이 다시 쓸 수 있습니다. \[\phi(B) (1-B)^d y_t = c + \theta(B) \varepsilon_t\] 여기에서 \(\phi(B)= (1-\phi_1B - \cdots - \phi_p B^p)\)은 \(B\)에 대한 \(p\)차 다항식이고, \(\theta(B) = (1 + \theta_1 B + \cdots + \theta_q B^q)\)은 \(B\)에 대한 \(q\)차 다항식입니다.
모델에 대한 정상성(stationarity) 조건은 단위 원(unit circle) 바깥에 있는 \(\phi(B)\)의 \(p\)개 복소근(complex root)들이고, 역변환 가능성 조건은 단위 원 바깥에 있는 \(\theta(B)\)의 \(q\)개 복소근입니다. 따라서 복소 단위 원(complex unit circle)과 여러 근의 관계를 나타내서 모델이 역변환 가능성과 정상성을 나타내는지 확인할 수 있습니다.
대신 역 근(inverse root)은 모두 단위 원 안에 위치해야하기 때문에 그림으로 나타내는 것이 좀 더 쉽습니다. 이 작업을 R로 쉽게 할 수 있습니다. 계절성으로 조정된 전자 제품 지수를 맞추는 ARIMA(3,1,1) 모델에 대해, 그림 8.16 을 얻게 됩니다.
autoplot(fit)
Figure 8.16: ARIMA(3,1,1) 모델을 계절성 조정된 전자 장비 지수에 맞춘 것에 대한 역 특성근.
왼쪽 그래프에서 3개의 빨간 점은 다항식 \(\phi(B)\)의 근과 관계 있고, 오른쪽 그래프에서 빨간 점은 \(\theta(B)\)의 근과 관계 있습니다. R에서는 적합 모델(fitted model)이 반드시 정상성을 나타내면서 역변환 가능하도록 만들기 때문에 예상한 것처럼 빨간 점 모두 단위 원 안에 있습니다. 단위 원에 가까운 근은 수치적으로 불안정한 것일 수 있고, 관련 모델은 예측하기에 적합하지 않을 것입니다.
Arima() 함수는 단위 원 바깥의 역 근(inverse root)에 해당하는 모델을 절대 결과로 내지 않을 것입니다. auto.arima()은 좀 더 엄격해서, 단위 원에 가까운 근에 해당하는 모델까지도 고르지 않을 것입니다.