10.7 최적 조정 접근법

일관된 예측값에서 예측 오차를 최소화하는 \(\bm{G}\) 행렬을 찾을 수 있다면, 최적 예측 조정(optimal forecast reconciliation)이 나타날 것입니다. 여기에서는 이 접근 방식을 간단하게 요약하여 설명하겠습니다. 자세한 내용은 Wickramasuriya et al. (2019) 에 있습니다.

(10.6)을 이용하여 일관된 예측값을 내는 상황을 생각해봅시다. 편의상 여기에서 이 상황을 반복하였습니다:

\[ \tilde{\bm{y}}_h=\bm{S}\bm{P}\hat{\bm{y}}_h. \] 먼저 편향되지 않은 예측값을 얻었는지 확인할 필요가 있습니다. 기준 예측값 \(\hat{\bm{y}}_h\)이 편향되지 않았다면, 일관된 예측값(coherent forecast) \(\tilde{\bm{y}}_h\)이 편향되지 않을 것이고 \(\bm{S}\bm{G}\bm{S}=\bm{S}\) 로 주어집니다 (Hyndman et al., 2011). 이것이 행렬 \(\bm{G}\)에 대한 제약 조건을 줍니다. 재미있게도, 모든 하향식 기법(top-down)은 이 제약 조건을 만족하지 않아서, 모든 하향식 기법은 편향됩니다.

다음, 예측값에서 오차를 찾을 필요가 있습니다. Wickramasuriya et al. (2019) 에서는 \(h\)-단계-앞 일관된 예측 오차(coherent forecast error)의 분산-공분산(variance-covariance) 행렬이 다음과 같이 주어지는 것을 보였습니다. \[\begin{equation*} \bm{V}_h = \text{Var}[\bm{y}_{T+h}-\tilde{\bm{y}}_h]=\bm{S}\bm{G}\bm{W}_h\bm{G}'\bm{S}' \end{equation*}\] 여기에서 \(\bm{W}_h=\text{Var}[(\bm{y}_{T+h}-\hat{\bm{y}}_h)]\)는 관련된 기준 예측 오차(base forecast error)의 분산-공분산 행렬입니다.

목표는 수정된 예측값의 오차 분산을 최소화하는 최적 행렬 \(\bm{G}\)을 찾는 것입니다. Wickramasuriya et al. (2019)\(\bm{V}_h\)의 대각합(trace)을 최소화하는 \(\bm{S}\bm{G}\bm{S}=\bm{S}\)인 최적 행렬 \(\bm{G}\)이 다음과 같이 주어지는 것을 보였습니다. \[ \bm{G}=(\bm{S}'\bm{W}_h^{-1}\bm{S})^{-1}\bm{S}'\bm{W}_h^{-1}. \] 따라서 최적으로 조정된 예측값(optimal reconciled forecasts)은 다음과 같이 주어집니다. \[\begin{equation} \tag{10.8} \tilde{\bm{y}}_h=\bm{S}(\bm{S}'\bm{W}_h^{-1}\bm{S})^{-1}\bm{S}'\bm{W}_h^{-1}\hat{\bm{y}}_h. \end{equation}\] 이것을 MinT (또는 최소 대각합-Minimum Trace) 추정값으로 부르겠습니다.

이것을 실제로 사용하기 위해, \(h\)-단계-앞 기준 예측값의 예측 오차 분산인 \(\bm{W}_h\)을 추정해야 합니다. 이것은 어려운 작업일 수 있습니다. 그래서 시뮬레이션과 실제 상황 두 경우에서 모두 잘 작동한다고 알려진 4가지 근사 기법을 소개하려고 합니다.

  1. \(\bm{W}_h=k_h\bm{I}\) \(\forall h\)로 둡시다. 여기에서 \(k_{h} > 0\)입니다.21 이것은 만들기 위해 가장 단순화한 가정입니다. 내용은 다음과 같습니다. \(\bm{P}\)는 데이터와 관계 없다고 하여 상당한 계산량을 줄일 수 있도록 합니다. 하지만, 단점은 구조의 수준 사이의 눈금 차이나 시계열 사이의 관계를 고려하지 못한다는 것입니다. 이 방식은 forecast() 함수에서 method = "comb"weights = "ols"로 두고 구현합니다.

    (10.8)에서 \(\bm{W}_h=k_h\bm{I}\)로 두면 5.7 절에서 \(\bm{X}=\bm{S}\)\(\bm{y}=\hat{\bm{y}}\)로 소개했던 최소 제곱 추정량이 되기 때문에, 여기에서 가중치를 OLS(ordinary least square; 보통의 최소 제곱) 로 부릅니다.

  2. 모든 \(h\)에 대해 \(\bm{W}_{h} = k_{h}\text{diag}(\hat{\bm{W}}_{1})\) 로 둡시다. 여기에서 \(k_{h} > 0\)이고, \[ \hat{\bm{W}}_{1} = \frac{1}{T}\sum_{t=1}^{T}\bm{e}_{t}\bm{e}_{t}', \] 입니다. 그리고 \(\bm{e}_{t}\)는 데이터와 동일한 순서로 누적된 기준 예측값을 낸 모델의 잔차의 \(n\)-차원 벡터입니다. forecast() 함수로 method = "comb"weights = "wls"로 두고 이 방식을 구현합니다.

    잔차의 분산을 가지고 기준 예측값의 눈금을 맞춥니다. 그래서 이것을 분산 눈금 잡기(variance scaling) 방법이라고 부릅니다.

  3. 모든 \(h\)에 대해 \(\bm{W}_{h}=k_{h}\bm{\Lambda}\)로 두고 여기에서 \(k_{h} > 0\)이고, \(\bm{\Lambda}=\text{diag}(\bm{S}\bm{1})\)로 둡시다. \(\bm{1}\)은 차원 \(m\) (밑바닥 수준 시계열의 수)의 단위 벡터입니다. 이 방법에서는 밑바닥 수준(bottom-level) 기준예측(base forecast) 오차 각각이 분산 \(k_{h}\)를 갖고 마디(node) 사이에 상관관계(correlation)가 없다는 것을 가정합니다. 따라서 \(\bm{\Lambda}\) 행렬의 대각성분(diagonal element) 각각은 각 마디에 기여하는 예측 오차 분산의 수를 포함합니다. 이 추정량은 실제 데이터에 의존하지 않고, 합산 구조에만 의존합니다. 따라서 이 방법을 구조적 눈금 잡기(structural scaling)라고 부릅니다. 구조적 눈금 잡기를 적용하면 잔차(residual)를 구할 수 없는 경우에 특히 유용합니다. 따라서 분산 눈금잡기를 사용할 수 없습니다. 판단 예측(4 장)으로 낸 기준 예측값이 이와 같은 상황이 적용되는 예가 되겠습니다. 이 방법은 method = "comb"weights = "nseries"로 두고 구현합니다.

  4. 모든 \(h\)에 대해 \(\bm{W}_h = k_h \bm{W}_1\)로 둡시다. 여기에서 \(k_h>0\)입니다. 여기에서는 오차 공분산 행렬(error covariance matrix)이 서로 비례한다고만 가정합니다. 전체 1단계 공분산 행렬 \(\bm{W}_1\)을 직접 계산합니다. 가장 분명하고 단순한 방법은 표본 공분산(sample covariance)을 사용하는 것입니다. 이 방식은 forecast() 함수에서 method = "comb", weights = "mint", covariance = "sam"로 두고 구현합니다.

    하지만, 밑바닥 수준 시계열의 수 \(m\)이 시계열의 길이 \(T\)에 비해 큰 경우에는, 이것은 좋은 추정량이 아닙니다. 대신에 표본 공분산을 대각행렬로 만드는 축소추정량(shrinkage estimator)을 사용합니다. 이것은 method = "comb", weights = "mint", covariance = "shr"로 두고 구현합니다.

정리하면, 존재하는 다른 기법과는 다르게, 최적 조정 예측값(optimal reconciliation forecast)은 계층적 구조 또는 그룹화된 구조 안에 있는 모든 사용 가능한 정보를 가지고 냅니다. 이 차이는 중요합니다. 특정한 합산 수준이나 그룹이, 사용자가 관심 있어 할 그리고 모델링 과정에서 중요하게 사용될 수 있는 데이터의 특징을 나타낼 수 있기 때문입니다.

예를 들면, 10.1 절에서 소개한 호주 여행자 데이터를 살펴봅시다. 여기에서는 계층 구조가 국가의 지리적인 구분 즉, 주와 지역으로 나타납니다. 몇몇 해안가 지역은 여름 휴양지로 유명해서 몰릴 것이지만, 몇몇 산악 지역은 겨울 휴양지로 몰릴 것입니다. 이러한 차이는 합산하면 국가 수준에서는 묻힐 것입니다.

예제: 호주 수감 인구 예측

10.2 절에서 다룬 호주 수감 인구에 대한 예측값을 계산해보겠습니다. forecast() 함수의 기본 입력값을 사용하여, 분산 눈금을 사용하는 WLS 추정량과 최적 조정 접근법을 이용하여 일관된 예측값을 계산하겠습니다.

prisonfc <- forecast(prison.gts)

합산의 각 수준에 대한 예측값을 얻기 위해, aggts() 함수를 사용할 수 있습니다. 예를 들면, 종합적인 전체 수감 인구와 1개 요인에 대한(주, 성별, 법적 상태) 예측값을 계산하기 위해, 다음을 이용합니다:

fcsts <- aggts(prisonfc, levels=0:3)

다음과 같은 내용을 이용하여 단순한 그래프를 얻습니다.

groups <- aggts(prison.gts, levels=0:3)
autoplot(fcsts) + autolayer(groups)

다음의 코드를 이용하여 더 나은 그래프를 얻을 수 있습니다. 결과는 그림 10.8 에 있습니다. 수직 선은 예측 시기의 시작을 나타냅니다.

prisonfc <- ts(rbind(groups, fcsts),
  start=start(groups), frequency=4)
p1 <- autoplot(prisonfc[,"Total"]) +
  ggtitle("호주 수감 인구") +
  xlab("연도") + ylab("전체 수감자 수 ('000)") +
  geom_vline(xintercept=2017)
cols <- sample(scales::hue_pal(h=c(15,375),
          c=100,l=65,h.start=0,direction = 1)(NCOL(groups)))
p2 <- as_tibble(prisonfc[,-1]) %>%
  gather(Series) %>%
  mutate(Date = rep(time(prisonfc), NCOL(prisonfc)-1),
         Group = str_extract(Series, "([A-Za-z ]*)")) %>%
  ggplot(aes(x=Date, y=value, group=Series, colour=Series)) +
    geom_line() +
    xlab("연도") + ylab("수감자 수 ('000)") +
    scale_colour_manual(values = cols) +
    facet_grid(. ~ Group, scales="free_y") +
    scale_x_continuous(breaks=seq(2006,2018,by=2)) +
    theme(axis.text.x = element_text(angle=90, hjust=1)) +
    geom_vline(xintercept=2017)
gridExtra::grid.arrange(p1, p2, ncol=1)
전체 호주 성인 수감 인구와 인구를 주, 법적 상태, 성별로 그룹화한 것에 대한 일관된 예측값.

Figure 10.8: 전체 호주 성인 수감 인구와 인구를 주, 법적 상태, 성별로 그룹화한 것에 대한 일관된 예측값.

그림 10.9 을 그리는데 비슷한 코드를 사용했습니다. 왼쪽 패널은 주와 성별 사이의 상호작용에 대한 일관된 예측값(coherent forecast)을 나타냅니다. 오른쪽 패널은 밑바닥 수준(bottom-level) 시계열에 대한 예측값을 나타냅니다.

호주 성인 수감 인구를 모든 가능한 속성의 조합으로 그룹화한 것에 대한 일관된 예측값.

Figure 10.9: 호주 성인 수감 인구를 모든 가능한 속성의 조합으로 그룹화한 것에 대한 일관된 예측값.

acurracy() 명령어는 계층적이거나 그룹화된 구조에 걸쳐 예측 정확도를 계산할 때 유용합니다. 다음에 나오는 표에서는 상향식 방법과 최적 조정 기법으로 2015년 1분기부터 2016년 4분기까지 테스트 기간을 예측한 결과를 요약합니다.

결과를 통해 특별히 꼭대기 수준에서 최적 조정 기법이 더 정확한 예측값을 내는 것을 알 수 있습니다. 일반적으로 최적 조정 기법이 구조 안에 있는 모든 수준의 정보를 사용하기 때문에 제한된 정보를 사용하는 다른 전통적인 방법보다 더 정확한 일관된 예측값을 내는 것을 알 수 있습니다.

train <- window(prison.gts, end = c(2014, 4))
test <- window(prison.gts, start = 2015)

fcsts.opt <- forecast(train, h=8, method="comb",
                      weights="wls", fmethod="ets")
fcsts.bu <- forecast(train, h=8, method="bu",
                     fmethod="ets")
tab <- matrix(NA, ncol = 4, nrow = 6)
rownames(tab) <- c("전체", "주", "법적 상태", "성별", "밑바닥", "모든 시계열")
colnames(tab) <- c("MAPE", "MASE", "MAPE", "MASE")
tab[1, ] <- c(
  accuracy(fcsts.bu, test, levels = 0)[c("MAPE", "MASE"), "Total"],
  accuracy(fcsts.opt, test, levels = 0)[c("MAPE", "MASE"), "Total"]
)
j <- 2
for (i in c(1:3, 7)) {
  tab[j, ] <- c(
    mean(accuracy(fcsts.bu, test, levels = i)["MAPE", ]),
    mean(accuracy(fcsts.bu, test, levels = i)["MASE", ]),
    mean(accuracy(fcsts.opt, test, levels = i)["MAPE", ]),
    mean(accuracy(fcsts.opt, test, levels = i)["MASE", ])
  )
  j <- j + 1
}
tab[6, ] <- c(
  mean(accuracy(fcsts.bu, test)["MAPE", ]),
  mean(accuracy(fcsts.bu, test)["MASE", ]),
  mean(accuracy(fcsts.opt, test)["MAPE", ]),
  mean(accuracy(fcsts.opt, test)["MASE", ])
)
out <- knitr::kable(tab, digits=2, booktabs=TRUE,
    format = ifelse(html, "html", "latex"),
    caption="서로 다른 시계열 그룹에 대한 호주 수감 인구 예측값의 정확도.") %>%
    kableExtra::add_header_above(c(" "=1, "Bottom-up" = 2, "Optimal" = 2))
if(!html) {
  out <- gsub("\\[t\\]","\\[b\\]",out)
}
out
Table 10.1: 서로 다른 시계열 그룹에 대한 호주 수감 인구 예측값의 정확도.
Bottom-up
Optimal
MAPE MASE MAPE MASE
전체 5.32 1.84 3.08 1.06
7.59 1.88 7.62 1.85
법적 상태 6.40 1.76 4.32 1.14
성별 8.62 2.68 8.72 2.74
밑바닥 15.82 2.23 15.25 2.16
모든 시계열 12.41 2.16 12.02 2.08

참고 문헌

Hyndman, R. J., Ahmed, R. A., Athanasopoulos, G., & Shang, H. L. (2011). Optimal combination forecasts for hierarchical time series. Computational Statistics and Data Analysis, 55(9), 2579–2589. https://robjhyndman.com/publications/hierarchical/

Wickramasuriya, S. L., Athanasopoulos, G., & Hyndman, R. J. (2019). Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. J American Statistical Association, 114(526), 804–819. https://robjhyndman.com/publications/mint/


  1. \(k_{h}\)는 비례 상수입니다. (10.8) 에서 사라지기 때문에 이 값은 추정하지 않아도 되고 여기에서 정할 필요가 없습니다.↩︎