8.9 계절성 ARIMA 모델들

지금까지, 비-계절성 데이터와 비-계절성 ARIMA 모델에만 관심을 두었습니다. 하지만, ARIMA 모델로 다양한 종류의 계절성 데이터를 모델링 할 수도 있습니다.

계절성 ARIMA 모델은 지금까지 살펴본 ARIMA 모델에 추가적인 계절성 항을 포함하여 구성됩니다. 계절성 ARIMA 모델을 다음과 같이 쓸 수 있습니다:

ARIMA \(\underbrace{(p, d, q)}\) \(\underbrace{(P, D, Q)_{m}}\)
\({\uparrow}\) \({\uparrow}\)
모델의 모델의
비 계절성 부분 계절성 부분

여기에서 \(m\)은 매년 관측값의 개수입니다. 모델의 계절성 부분에 대문자 기호를 사용하고, 모델의 비-계절성 부분에서는 소문자 기호를 사용하겠습니다.

모델의 계절성 부분은 모델의 비-계절성 성분과 비슷한 항으로 구성됩니다만, 계절성 주기의 후방이동(backshift)을 포함합니다. 예를 들어, (상수가 없는) ARIMA(1,1,1)(1,1,1)\(_{4}\) 모델은 분기별 데이터(\(m=4\))에 대한 경우이고 다음과 같이 쓸 수 있습니다. \[ (1 - \phi_{1}B)~(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} = (1 + \theta_{1}B)~ (1 + \Theta_{1}B^{4})\varepsilon_{t}. \]

단순히 추가적인 계절성 항을 비-계절성 항에 곱합니다.

ACF/PACF

AR이나 MA 모델의 계절성 부분은 PACF와 ACF의 계절성 시차에서 볼 수 있을 것입니다. 예를 들면, ARIMA(0,0,0)(0,0,1)\(_{12}\) 모델은 다음과 같은 모습을 나타낼 것입니다:

  • ACF에 시차 12에서 나타나는 뾰족한 막대는 있지만 다른 유의미한 뾰족한 막대는 없는 모습;
  • PACF의 계절성 시차가 지수적으로 감소하는 모습(즉, 시차 12, 24, 36, …).

비슷하게, ARIMA(0,0,0)(1,0,0)\(_{12}\) 모델은 다음과 같은 모습을 나타낼 것입니다:

  • ACF의 계절성 시차가 지수적으로 감소하는 모습;
  • PACF에 시차 12에서 나타나는 유의미한 뾰족한 막대가 있는 모습.

계절성 ARIMA 모델에서 적절한 계절성 차수를 고려할 때, 계절성 시차를 고려하는 것에 주의하시길 바랍니다.

모델링 과정은 모델의 비-계절성 성분과 계절성 AR과 MA 항을 골라야 하는 것을 제외하고는 비-계절성 데이터에 대한 경우와 거의 같습니다. 예제를 통해 이러한 과정을 잘 익힐 수 있습니다.

예제: 유럽 분기별 소매 거래

1996년부터 2011년까지 분기별 유럽 소매 거래 데이터를 이용하여 계절성 ARIMA 모델링 과정을 설명하겠습니다. 그림 8.17에 데이터를 나타냈습니다.

autoplot(euretail) + ylab("소매 지수") + xlab("연도")
유로 지역 (17개 국가) 1996--2011 기간 동안 도매와 소매 거래를 포함하고, 자동차와 오토바이 수리에 대한 분기별 소매 거래 지수. (지수: 2005 = 100).

Figure 8.17: 유로 지역 (17개 국가) 1996–2011 기간 동안 도매와 소매 거래를 포함하고, 자동차와 오토바이 수리에 대한 분기별 소매 거래 지수. (지수: 2005 = 100).

데이터가 정상성을 나타내지 않는다(non-stationary)는 것을 분명하게 알 수 있고 약간의 계절성도 보입니다. 따라서 먼저 계절성 차분을 구할 것입니다. 계절성으로 차분을 구한 데이터를 그림 8.18에 나타냈습니다. 이것도 정상성을 나타내지 않는 것 같기 때문에, 1차 차분을 한 번 더 구하여 그림 8.19에 나타냈습니다.

euretail %>% diff(lag=4) %>% ggtsdisplay()
계절성 차분을 구한 유럽 소매 거래 지수..

Figure 8.18: 계절성 차분을 구한 유럽 소매 거래 지수..

euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()
 두 번 차분을 구한 유럽 소매 거래 지수.

Figure 8.19: 두 번 차분을 구한 유럽 소매 거래 지수.

이제 그림 8.19에 나타낸 ACF와 PACF에 기초하여 적절한 ARIMA를 찾아봅시다. ACF에서 시차 1의 유의미하게 뾰족한 막대가 비-계절성 MA(1) 성분을 암시하고, ACF에서 시차 4의 유의미하게 뾰족한 막대는 계절성 MA(1) 성분을 암시합니다. 이를 통해, 1차 차분과 계절성 차분을 나타내는 ARIMA(0,1,1)(0,1,1)\(_4\) 모델과 비-계절성과 계절성 MA(1) 성분을 가지고 시작합니다. 그림 8.20에 맞춘 모델에 대한 잔차를 나타내었습니다. (비슷한 방법으로 PACF를 적용하면, ARIMA(1,1,0)(1,1,0)\(_4\) 모델로 시작할 수 있습니다.)

euretail %>%
  Arima(order=c(0,1,1), seasonal=c(0,1,1)) %>%
  residuals() %>% ggtsdisplay()
유럽 소매 거래 지수 데이터에 대해 맞춘 ARIMA(0,1,1)(0,1,1)$_4$ 모델의 잔차.

Figure 8.20: 유럽 소매 거래 지수 데이터에 대해 맞춘 ARIMA(0,1,1)(0,1,1)\(_4\) 모델의 잔차.

ACF와 PACF 둘 다 시차 2에서 유의미하게 뾰족한 막대가 나타나고 시차 3에서 덜 유의미하지만 뾰족한 막대가 나타나는 것을 보여줍니다. 이것은 몇몇 추가적인 비-계절성 항이 모델에 추가되어야 한다는 것을 의미합니다. ARIMA(0,1,3)(0,1,1)\(_4\) 모델의 AICc 값은 68.53인데, ARIMA(0,1,2)(0,1,1)\(_4\) 모델의 AICc는 74.36입니다. AR 항이 있는 다른 모델도 시도해봤지만, AICc 값이 더 작은 경우는 없었습니다. 결론적으로, ARIMA(0,1,3)(0,1,1)\(_4\) 모델을 골랐습니다. 그림 8.21에 잔차를 나타냈습니다. 모든 뾰족한 막대가 이제 유의미한 범위 안에 들어오고, 따라서 잔차는 백색잡음처럼 보입니다. 융-박스(Ljung-Box) 검정 결과도 잔차에 자기상관값이 남아있지 않다는 것을 보여줍니다.

fit3 <- Arima(euretail, order=c(0,1,3), seasonal=c(0,1,1))
checkresiduals(fit3)
유럽 소매 거래 지수 데이터에 대해 맞춘 ARIMA(0,1,3)(0,1,1)$_4$ 모델의 잔차.

Figure 8.21: 유럽 소매 거래 지수 데이터에 대해 맞춘 ARIMA(0,1,3)(0,1,1)\(_4\) 모델의 잔차.

#> 
#> 	Ljung-Box test
#> 
#> data:  Residuals from ARIMA(0,1,3)(0,1,1)[4]
#> Q* = 0.51, df = 4, p-value = 1
#> 
#> Model df: 4.   Total lags used: 8

따라서, 필요한 확인 절차를 거친 계절성 ARIMA 모델을 얻었고 이것으로 예측할 준비가 됐습니다. 모델로 얻은 다음 3년에 대한 예측값은 그림 8.22에 나타냈습니다. 두 번 차분을 구한 것 때문에 예측값이 데이터의 최근 추세를 따라갑니다. 크고 급격하게 증가하는 예측 구간은 소매 거래 지수가 어느때나 증가하기 또는 감소하기 시작할 수 있다는 것을 나타냅니다 — 반면에, 점 예측값은 하향세를 나타내지만, 예측구간에서 보면 데이터가 예측 기간 동안에 상향 추세를 갖는 것도 가능합니다.

fit3 %>% forecast(h=12) %>% autoplot() + ggtitle("ARIMA(0,1,3)(0,1,1)으로 얻은 예측값")
ARIMA(0,1,3)(0,1,1)\(_4\) 모델을 이용하여 얻은 유럽 소매 거래 지수 데이터 예측값. 80%와 95% 예측구간을 나타냈습니다.

Figure 8.22: ARIMA(0,1,3)(0,1,1)\(_4\) 모델을 이용하여 얻은 유럽 소매 거래 지수 데이터 예측값. 80%와 95% 예측구간을 나타냈습니다.

이 작업 대부분을 수행하기 위해 auto.arima()을 사용할 수도 있습니다. 이 함수를 사용하면 같은 결과를 얻습니다.

auto.arima(euretail)
#> Series: euretail 
#> ARIMA(0,1,3)(0,1,1)[4] 
#> 
#> Coefficients:
#>         ma1    ma2    ma3    sma1
#>       0.263  0.369  0.420  -0.664
#> s.e.  0.124  0.126  0.129   0.155
#> 
#> sigma^2 estimated as 0.156:  log likelihood=-28.63
#> AIC=67.26   AICc=68.39   BIC=77.65

auto.arima() 함수는 \(D\) (사용할 계절성 차분 횟수)를 결정하기 위해 nsdiffs()를 사용하고, \(d\) (보통 차분의 횟수)를 결정하기 위해 ndiffs()를 사용합니다. 다른 모델 매개변수는 (\(p,q,P,Q\)) 비-계절성 ARIMA 모델처럼 모두 AICc를 최소화하여 결정합니다.

예제: 호주 코르티코 스테로이드 약물 판매량

두 번째 예제는 좀 더 어렵습니다. 호주에서 코르티코 스테로이드 (Corticosteroid) 약물이 매월 팔린 양을 예측하려고 합니다. 해부학적 치료용 화학 약품 분류 체계(Anatomical Therapeutic Chemical classification scheme)에서 H02 약물들로 분류된 것으로 알려져있습니다.

lh02 <- log(h02)
cbind("H02 판매량 (백만 처방전)" = h02,
      "로그 H02 판매량"=lh02) %>%
  autoplot(facets=TRUE) + xlab("연도") + ylab("")
호주 코르티코 스테로이드 약물 판매량 (단위: 백만 처방전/월). 로그를 취한 값은 아래쪽에 나타냈습니다.

Figure 8.23: 호주 코르티코 스테로이드 약물 판매량 (단위: 백만 처방전/월). 로그를 취한 값은 아래쪽에 나타냈습니다.

그림 8.23에 1991년 7월부터 2008년 6월까지의 데이터를 나타냈습니다. 수준의 분산이 조금 증가합니다. 따라서 분산을 안정화시키기 위해 로그를 취하겠습니다.

데이터에는 계절성이 강하게 나타납니다. 그리고 데이터에 정상성이 나타나지 않는다는 것이 분명합니다. 그래서 계절성 차분을 사용할 것입니다. 계절성 차분을 구한 데이터를 그림 8.24에 나타냈습니다. 아직까지는 차분을 더 구해야하는지 여부를 분명하게 알기 어렵습니다. 차분을 더 계산하지 않기로 결정하지만, 확실한 선택은 아닙니다.

마지막 몇몇 관측값은 초반부의 데이터와 다른 것처럼 보입니다(더 변동이 큰 것 같은). 초반 판매량 보고가 늦어져서 데이터가 때때로 수정되는 것 때문에 생기는 현상일 수도 있습니다.

lh02 %>% diff(lag=12) %>%
  ggtsdisplay(xlab="연도",
    main="계절성 차분을 구한 H02 처방전 데이터")
계절성 차분을 구한 호주 코르티코 스테로이드 약물 판매량 (단위: 백만 처방전/월).

Figure 8.24: 계절성 차분을 구한 호주 코르티코 스테로이드 약물 판매량 (단위: 백만 처방전/월).

계절성 차분을 구한 데이터를 나타내는 그래프에서, PACF를 보니 시차 12와 시차 24에서 뾰족한 막대가 있습니다만, ACF를 보니 다른 계절성 시차에서는 나타나지 않습니다. 이 사실은 계절성 AR(2) 항을 암시합니다. 비-계절성 시차에서는, PACF를 보면 유의미하게 뾰족한 막대가 3개 있습니다. 이것은 AR(3) 항이 가능하다는 것을 암시합니다. ACF에서 나타나는 패턴은 어떠한 단순한 모델을 가리키지는 않습니다.

결과적으로, 초기 분석은 이러한 데이터에 가능한 모델이 ARIMA(3,0,0)(2,1,0)\(_{12}\)라는 것을 시사합니다. 이 모델에 몇몇 변형을 준 것과 함께 이 모델을 맞추고, AICc 값을 계산하여 다음의 표에 나타냈습니다.

Model AICc
ARIMA(3,0,1)(0,1,2)\(_{12}\) -485.5
ARIMA(3,0,1)(1,1,1)\(_{12}\) -484.2
ARIMA(3,0,1)(0,1,1)\(_{12}\) -483.7
ARIMA(3,0,1)(2,1,0)\(_{12}\) -476.3
ARIMA(3,0,0)(2,1,0)\(_{12}\) -475.1
ARIMA(3,0,2)(2,1,0)\(_{12}\) -474.9
ARIMA(3,0,1)(1,1,0)\(_{12}\) -463.4

이러한 모델 중에서 ARIMA(3,0,1)(0,1,2)\(_{12}\) 모델이 가장 좋습니다(즉, 가장 작은 AICc값을 갖습니다).

(fit <- Arima(h02, order=c(3,0,1), seasonal=c(0,1,2),
  lambda=0))
#> Series: h02 
#> ARIMA(3,0,1)(0,1,2)[12] 
#> Box Cox transformation: lambda= 0 
#> 
#> Coefficients:
#>          ar1    ar2    ar3    ma1    sma1    sma2
#>       -0.160  0.548  0.568  0.383  -0.522  -0.177
#> s.e.   0.164  0.088  0.094  0.190   0.086   0.087
#> 
#> sigma^2 estimated as 0.00428:  log likelihood=250
#> AIC=-486.1   AICc=-485.5   BIC=-463.3
checkresiduals(fit, lag=36)
ARIMA(3,0,1)(0,1,2)$_{12}$ 모델을 H02 월별 처방전 판매량 데이터에 적용한 결과의 잔차.

Figure 8.25: ARIMA(3,0,1)(0,1,2)\(_{12}\) 모델을 H02 월별 처방전 판매량 데이터에 적용한 결과의 잔차.

#> 
#> 	Ljung-Box test
#> 
#> data:  Residuals from ARIMA(3,0,1)(0,1,2)[12]
#> Q* = 51, df = 30, p-value = 0.01
#> 
#> Model df: 6.   Total lags used: 36

이 모델에서 잔차는 그림 8.25에 나타냈습니다. ACF를 보면 유의미하게 뾰족한 막대가 몇 개 있고, 모델은 융-박스(Ljung-Box) 검정을 통과하지 못했습니다. 여전히 모델을 예측에 사용할 수 있지만, 보정된 잔차 때문에 예측구간이 정확하지 않을 수 있습니다.

다음은 자동 ARIMA 알고리즘을 사용하여 시도해볼 것입니다. 모든 입력값을 기본 설정값으로 두고 auto.arima()를 돌리면 ARIMA(2,1,1)(0,1,2)\(_{12}\) 모델을 얻습니다. 하지만, 이 모델은 여전히 시차 36에 대해 융-박스(Ljung-Box) 검정을 통과하지 못합니다. 때때로 모든 검정을 통과하는 모델을 찾지 못할 수도 있습니다.

테스트 데이터 평가:

데이터의 마지막 2년으로 구성된 테스트 데이터를 이용하여 지금까지 맞춘 모델 중 몇 가지를 비교하려고 합니다. 따라서, 1991년 7월부터 2006년 6월까지의 데이터를 이용하여 모델을 맞추고, 2006년 7월부터 2008년 6월까지의 판매량을 예측하겠습니다. 결과는 표 8.2에 요약하였습니다.

Table 8.2: 다양한 ARIMA 모델을 H02 월별 처방전 판매량 데이터에 적용한 결과에 대한 RMSE 값.
모델 RMSE
ARIMA(3,0,1)(0,1,2)\(_{12}\) 0.0622
ARIMA(3,0,1)(1,1,1)\(_{12}\) 0.0630
ARIMA(2,1,3)(0,1,1)\(_{12}\) 0.0634
ARIMA(2,1,1)(0,1,2)\(_{12}\) 0.0634
ARIMA(2,1,2)(0,1,2)\(_{12}\) 0.0635
ARIMA(3,0,3)(0,1,1)\(_{12}\) 0.0637
ARIMA(3,0,1)(0,1,1)\(_{12}\) 0.0644
ARIMA(3,0,2)(0,1,1)\(_{12}\) 0.0644
ARIMA(3,0,2)(2,1,0)\(_{12}\) 0.0645
ARIMA(3,0,1)(2,1,0)\(_{12}\) 0.0646
ARIMA(4,0,2)(0,1,1)\(_{12}\) 0.0648
ARIMA(4,0,3)(0,1,1)\(_{12}\) 0.0648
ARIMA(3,0,0)(2,1,0)\(_{12}\) 0.0661
ARIMA(3,0,1)(1,1,0)\(_{12}\) 0.0679

RMSE 값을 기준으로 볼 때 직접 고른 모델과 auto.arima() 둘 다 상위 4개의 모델에 포함됩니다.

AICc 값으로 모델을 비교할 때, 모든 모델을 같은 차수로 차분을 구한 것이 중요합니다. 하지만, 테스트 데이터을 사용하여 모델을 비교할 때, 예측값이 어떻게 나왔는지가 중요하지 않습니다 — 이러한 비교가 항상 유효합니다. 결과적으로, 위의 표에서, 계절성 차분만 고려한 몇몇 모델과 1차 차분과 계절성 차분 둘 다 고려한 몇몇 모델을 포함시킬 수 있습니다. 반면에, AICc 값을 포함하는 이전의 표에서는 1차 차분은 비교 대상에 포함하지 않고 계절성 차분을 고려한 모델만 비교했습니다.

살펴본 모델 중에서 어떠한 것도 잔차 검정을 전부 통과하지 못했습니다. 실제 상황에서, 보통은 모델이 모든 검정을 통과하지 못하더라도 찾을 수 있었던 가장 좋은 모델을 사용합니다.

(테스트 데이터에 대해 가장 낮은 RMSE 값과 계절성 차분을 고려한 모델 중에서 가장 좋은 AICc 값을 갖는) ARIMA(3,0,1)(0,1,2)\(_{12}\) 모델로 얻은 예측값을 그림 8.26에 나타냈습니다.

h02 %>%
  Arima(order=c(3,0,1), seasonal=c(0,1,2), lambda=0) %>%
  forecast() %>%
  autoplot() +
    ylab("H02 판매량 (백만 처방전)") + xlab("연도") +
    ggtitle("ARIMA(3,0,1)(0,1,2)로 얻은 예측값")
ARIMA(3,0,1)(0,1,2)$_{12}$ 모델을 H02 월별 처방전 판매량 데이터에 적용한 결과로부터 얻은 예측값.

Figure 8.26: ARIMA(3,0,1)(0,1,2)\(_{12}\) 모델을 H02 월별 처방전 판매량 데이터에 적용한 결과로부터 얻은 예측값.