계층적 시계열 분석

❗️블로그 옮김:  https://www.taemobang.com

방태모

최근 데이터 수집과 저장 기술의 발전으로 계층 구조 또는 그룹화 구조를 갖는 많은 다중 시계열 자료를 접할 수 있게되었다. 이러한 자료의 각 계열을 단일 시계열 예측 모형(e.g. ARIMA 모형)을 이용해 각 계열을 예측하는 경우 예측값들이 계층 구조(또는 그룹화 구조)를 만족하지 않는 문제가 발생한다. 이러한 예측값에 계층적 시계열 예측 방법을 적용할 경우, 예측값들이 계층 구조를 만족할 뿐만아니라, 예측 성능의 개선에도 도움을 주게 된다.

 

계층적 또는 그룹화 구조로 표현되는 시계열을 각각 계층적 시계열, 그룹화 시계열이라고 하며, 다음과 같이 정의한다.

 

 ○ 계층적 시계열(Hierarchical time series)

 계층적 구조가 고유하게 존재하는 시계열을 말한다. 예를 들면, 광역 철도의 수송량은 철도 노선 별로 구분되며, 이는 다시 역별로 세분화될 수 있다. 

 

 ○ 그룹화 시계열(Grouped time series)

 고유한 계층 구조가 없는 시계열로, 그룹화할 수 있는 순서가 고유하지 않다. 예를 들어, 한 나라의 사망률은 지역으로 분해된 다음 성별로 분해되거나, 그 반대의 순서로도 분해가 가능하다.

 

출처 - G. Athanasopoulos et al. (2009)

예를 들어, Fig. 1은 2-수준 계층 구조를 갖는 계층적 시계열을 도식화하여 표현한 것이다. 계층적 시계열 예측에서는 계층 구조를 구성하고 있는 각 계열들을 노드라고 표현한다. 위 그림의 Total을 총 계열이라 하며, 해당 층을 최상위 수준(level) 또는 수준 0이라고 한다. 그 다음으로 분해되는 층을 수준 1이라고 하며, 그림의 마지막 층을 최하위 수준이라고 표현한다. 즉, 최상위 수준을 제외한 2개의 수준이 존재하므로, 2-수준 계층 구조를 갖는 시계열이라고 표현한다.

 

이 글에서는 4가지 계층적 시계열 예측 방법에 대해 소개한다. 본격적인 방법론의 설명에 앞서, 계층적 시계열을 다룰 때 알아두어야할 기본적인 표기법과 개념에 대해 정리한다. 먼저, 기본적인 표기법을 정리하자.

 

· $Y_{X, t}$ : 노드 $X$의 $t$시점 관측치

· $Y_t$ : 총 계열의 $t$시점 관측치

· $m_i$ : 수준 $i$에 속하는 노드의 갯수. 즉, $K$개 계층이 존재하는 계층 구조 전체 시계열의 수를 $m$이라 하면 다음의 식이 성립한다.

 

$m = m_0 + m_1 + \cdots + m_k$

 

· $\boldsymbol{Y}_{i, t}$ : 수준 $i$의 모든 노드들의 $t$ 관측치를 요소로 갖는 벡터.

 

위 표기법을 이용해 모든 노드의 $t$시점 관측치를 요소로 갖는 벡터 $\boldsymbol{Y}_t$는 다음과 같이 쓸 수 있다.

 

$\boldsymbol{Y}_t = {[Y_t, {\boldsymbol{Y}}'_{1, t}, \cdots, {\boldsymbol{Y}}'_{K,t}]}'$

 

또한, $\boldsymbol{Y}_t$는 합산 행렬(summing matrix) $\boldsymbol{S}$를 이용하면 다음과 같이 최하위 수준 노드들의 합으로 표현할 수 있다.

 

$\boldsymbol{Y}_t = \boldsymbol{S} \boldsymbol{Y}_{K,t}$  $\cdots$  (Eq 1)

 

예를 들어, (Eq 1)을 이용하여 Fig. 1을 표현하면 다음과 같다.

 

출처 - G. Athanasopoulos et al. (2009)

계층적 시계열 분석에서 우리가 결국 관심이 있는 것은 계층 구조를 반영한 예측이다.

 

· $\hat{Y}_{X,n}(h)$ : $n$개의 관측치를 기반으로한 노드 $X$의 $h$시점 이후 예측값

· $\hat{\boldsymbol{Y}}_{i, n}(h)$ : $n$개의 관측치를 기반으로한 수준 $i$의 모든 노드들의 $h$시점 이후 예측값을 요소로 가지는 벡터

· $\hat{\boldsymbol{Y}}_n (h)$ : $n$개의 관측치를 기반으로한 모든 노드들의 $h$시점 이후 예측값을 요소로 가지는 벡터

 

계층적 시계열 분석에서는 특히 $\hat{Y}$을 특히 기저 예측값(Base forecasts)이라 부른다. 이는 모든 계층적 예측 방법에서 필요로 되는 값으로, 각 시계열의 독립적인 예측을 통해 구해진다(e.g. 지수 평활법(Exponential smoothing), ARIMA models, etc.). 

 

모든 계층적 시계열 예측 방법은 다음과 같은 식으로 표현된다.

 

$\tilde{\boldsymbol{Y}}_n(h) = \boldsymbol{SP}\hat{\boldsymbol{Y}}_n(h)$  $\cdots$  (Eq 2)

 

$\tilde{\boldsymbol{Y}}_n(h)$이 임의의 계층적 시계열 예측 방법에 의해 만들어지는 최종적인 예측값이며, 이를 계층적 시계열 분석에서는 최종 수정된 예측값(The final revised forecasts)이라 한다. 즉, 최종 수정된 예측값들은 계층 구조를 만족할 것이다. 또한, (Eq 2)를 통해 결국 최종 수정된 예측값은 독립적인 기저 예측값들의 선형결합의 결과로 표현됨을 알 수 있다. $\boldsymbol{S}$는 (Eq 1)에서 봤던 합산 행렬이며, 해당 노드가 갖는 계층 구조에 따라 달라진다. 즉, 합산 행렬은 해당 자료의 계층 구조에 대한 정보를 갖는 행렬이라 할 수 있다. 그리고 $\boldsymbol{P}$는 계층적 시계열 예측 방법에 따라 달라진다. 이제 4가지 계층적 시계열 예측 방법에 대해 본격적으로 소개한다.

 

1. 상향식 방법(Bottom-up approach)

상향식 방법은 가장 간단한 계층적 시계열 예측 방법으로, 최하위 수준 노드들의 기저 예측값들의 합산을 통해 상위 수준의 최종 수정된 예측값을 생성한다.

 

$\boldsymbol{P} = [\boldsymbol{0}_{m_k \times (m-m_k)} | \boldsymbol{I}_{m_k}]$

 

여기서 $\boldsymbol{P}$의 역할은 최하위 수준의 기저 예측값들을 추출하는 것으로, 이들을 $\boldsymbol{S}$를 통해 합산하여 상위 수준의 노드들의 최종 수정된 예측값을 생성한다. 이 방법의 가장 큰 이점은 시계열이 최대로 분해된 지점인 최하위 수준으로부터 데이터를 모델링함으로써, 합산으로 인해 정보를 잃게되는 경우를 방지할 수 있다는 것이다. 이러한 이유로 계층 구조내의 계열들이 가지는 개개의 특성을 잘 잡아낼 수 있다. 그러나, 최하위 수준의 각 노드에 noise가 많은 경우, 상향식 방법은 좋지않은 결과를 초래할 수 있다.

 

2. 하향식 방법(Top-down approach)

하향식 방법은 총 계열의 기저 예측값을 일정 비율로 최하위 수준의 노드들로 분배하고, 합산을 통해 상위 수준 노드들의 최종 수정된 예측값을 생성한다.

 

$\boldsymbol{P} = [\boldsymbol{p} | \boldsymbol{0}_{m_k \times (m-1)}]$

 

$\boldsymbol{p} = {[p_1, p_2, \cdots, p_{m_k}]}'$로 총 계열의 기저 예측값을 최하위 수준의 각 노드로 분배할 때 쓰이는 일정 비율의 집합에 해당한다. 즉, $\boldsymbol{p}$의 원소는 최하위 수준 노드들의 갯수만큼 필요할 것이며, 여기서 $\boldsymbol{P}$의 역할은 총계열의 기저 예측값을 최하위 수준의 각 노드로 분배하는 것이다. Top-down approach는 이 일정 비율인 $p_j$를 계산하는 방식에 따라 세 가지 방법으로 나뉜다.

 

 (1) 과거 비율의 평균(The average historical proportions) $\rightarrow$ 과거 비율의 평균을 이용한 하향식 방법(Top-down HP1 method)

 

$p_j = \sum_{t=1}^n \frac{Y_{j, t}}{Y_t} / n$ for $j = 1, 2, \cdots, m_k$

 

 과거 비율의 평균은 총 계열의 과거값 대비 최하위 수준의 각 노드들의 과거값이 차지하는 비율로 총 계열의 기저 예측값을 분배한다. 주어진 과거 관측값에 의해 $p_j$가 결정되므로 이름에 과거(historical)라는 표현이 들어간다.

 

 (2) 과거 평균의 비율(proportions of the historical averages) $\rightarrow$ 과거 평균의 비율을 이용한 하향식 방법(Top-down HP2 method)

 

$p_j = \sum_{t=1}^n \frac{Y_{j, t}}{n} / \sum_{t=1}^n \frac{Y_t}{n}$ for $j = 1, 2, \cdots, m_k$

 

 과거 평균의 비율은 총 계열의 과거값의 평균 대비 최하위 수준 각 노드의 과거값 평균이 차지하는 비율을 통해 총 계열의 기저 예측값을 분배한다.

 

Top-down HP1, Top-down HP2는 각 $p_j$의 계산이 과거값에 기반을 두므로, 실제 예측 과정에서는 총 계열의 기저예측값만이 필요로 된다는 것이 장점이다. 두 방식은 합산 수준들(The aggregate levels)에서 꽤 신뢰할만한 예측 결과를 보여주며, Low count 데이터에서도 매우 효과적이다. 반대로, 이 두 방식이 가진 가장 큰 단점은 합산으로 인해 각 노드가 가지는 특성들을 잃을 수 있다는 것이다. 예를 들어, 지역별로 분할되는 어느 나라의 관광 수요 데이터를 떠올려보면, 스키 리조트가 있는 지역과 해수욕장이 있는 지역간의 계절적 패턴은 매우 다를 것이며, 이러한 경우 Top-down HP1과 Top-down HP2는 각 노드가 가지는 특성을 잃은채로 예측을 수행하게 될 것이다. 또 하나의 단점은, 데이터가 가지는 추세(Trends)를 놓치게 된다는 것이다.

 

이제 Top-down HP1과 Top-down HP2를 개선한 방법을 소개한다. 이 방식은 예측 비율(Forecasted proportions)을 기반으로 총 계열의 기저 예측값을 최하위 수준의 각 노드로 분배한다.

 

 (3) 예측 비율(Forecasted proportions) $\rightarrow$ 예측 비율을 이용한 하향식 방법(Top-down FP method)

다음의 식을 통해 계산되는 비율을 우리는 예측 비율이라고 부른다.

 

$p_j = \prod_{l=0}^{K-1} \frac{ {\hat{Y}}_{j,n}^{(l)} (h)}{ {\hat{S}}_{j, n}^{(l+1)} (h)}$  for $j = 1, 2, \cdots, m_k$

 

여기서 $(l)$은 노드 $j$의 상위 $l$수준, ${\hat{S}}_{j, n}^{(l+1)} (h)$는 노드 $j$의 상위 $(l+1)$수준과 연결된 하위 수준 모든 노드의 기저 예측값들의 합, {\hat{Y}}_{j,n}^{(l)} (h)은 노드 $j$의 상위 $(l)$수준 노드의 기저 예측값을 의미한다. 예를 들어 Fig. 1의 총 계열은 다음과 같이 쓸 수 있다.

 

$\hat{S}_{AA, n}^{(2)} (h) = \hat{S}_{Total, n} (h) = \hat{Y}_{A,n} (h) + \hat{Y}_{B,n} (h) + \hat{Y}_{C,n} (h)$

 

즉, 예측 비율은 해당 노드가 속한 수준의 기저 예측값들의 합 대비 해당 노드의 기저 예측값이 갖는 비율로 결정이 된다. 예측 비율을 이용한 상향식 방법은 앞선 두 방식과 달리, $p_j$의 계산에 기저 예측값을 활용하여 좀 더 세부적인 정보를 고려하지만, 여전히 기저 예측값들이 unbiased여도 최종 수정된 예측값은 unbiasedness를 만족하지 못하는 문제점을 안고있다. 이는 추후에 소개할 최적 조합 방법을 제외한 모든 계층적 시계열 예측 방법에서 발생하는 결점이다.

 

3. 미들 아웃 방법(Middle-out approach)

미들 아웃 방법은 중간 수준(middle-level)을 고르고, 위로는 상향식, 아래로는 예측 비율을 이용한 하향식 방법을 이용한다. Fig. 1과 같은 2-수준 계층 구조에 미들 아웃 방법을 적용하는 경우, 수준 1을 중간 수준으로 고르고 위로는 상향식, 아래로는 예측 비율을 이용한 하향식 방법을 사용하게 될 것이다.

 

4. 최적 조합 방법(Optimal combination approach)

최적 조합 방법은 계층적 시계열이 가진 계층 구조에 대한 정보를 가진 합산 행렬 $\boldsymbol{S}$를 예측 변수로 최하위 수준 노드들의 기저예측 값을 회귀적으로 추정하여, 상위 수준의 기저 예측값을 생성한다. 다음의 식이 이를 표현한 것이다.

 

$\hat{\boldsymbol{Y}}_n(h) = \boldsymbol{S\beta}_h + \epsilon_h$,  where $ \boldsymbol{\beta}_h = E[\hat{\boldsymbol{Y}}_{k, n}(h) | \boldsymbol{Y}_1, \cdots, \boldsymbol{Y}_n]$,  $E(\epsilon_h) = 0$,  $Var(\epsilon_h) = \sum_h$

 

이때 $\boldsymbol{\beta}_h = E[\hat{\boldsymbol{Y}}_{k, n}(h) | \boldsymbol{Y}_1, \cdots, \boldsymbol{Y}_n]$는 최하위 수준 $K$의 알려지지 않은 조건부 기댓값, ${\boldsymbol{S}}$는 $m \times m_k$ 합산 행렬, $\epsilon_{h}$는 선형 회귀에 대한 오차로 평균으로 0, 공분산 행렬로 $\sum_h$를 갖는다. 즉, 이 방식의 아이디어는 $h$시점 이후 기저 예측값들을 선형 회귀모형으로 표현하는 것에 있다. 이를 통해 계열간의 상관(correlations)과 상호작용(Interactions)을 고려할 수 있다. 그러나, 이 방식에는 한 가지 문제점이 존재한다. 일반적으로 실제 문제를 다룰 때에는 $\sum_h$에 대한 정보가 없기 때문에, $\boldsymbol{\beta}_h$의 최소분산 불편추정치(The minimum variance unbiased estimate)를 얻기위해 필요한 GLS(Generalized least squres) 추정을 사용하기 어렵다는 것이다(GLS에는 $\sum_h$의 추정이 필요로 됨). 다행히 이 문제점은 Hyndman et al. (2007)이 제안한 다음의 적절한 가정을 통해 해결된다.

 

$\epsilon_h \approx \boldsymbol{S}\epsilon_{K, h}$

 

이 가정에 따르면, 모든 노드의 $h$시점 이후 기저 예측값이 만들어내는 회귀모형의 잔차는 결국 최하위 수준 노드들의 $h$시점 이후 기저 예측값이 만들어내는 회귀모형 잔차의 합산으로 근사한다. 즉, 회귀모형의 잔차 또한 계층 구조를 만족해야 한다는 가정이다. 이를 이용하면 $\sum_h$에 의존하는 $\boldsymbol{\beta}_h$의 GLS 추정은 다음과 같이 매우 간단해진다:

 

$\hat{\boldsymbol{\beta}}_h = ({\boldsymbol{S}}'\boldsymbol{S})^{-1} {\boldsymbol{S}}'\hat{\boldsymbol{Y}}_n(h)$

 

이는, $\boldsymbol{\beta}_h$에 대한 BLUE(Best linear unibiased estimator)로 OLS 추정이 가능하게끔 해준다. 그러나, $\tilde{\boldsymbol{Y}}_n(h)$의 분산은 여전히 $\sum_h$로 신뢰구간에는 여전히 $\sum_h$에 대한 추정이 필요로 된다. 결론적으로, 최적 조합 방법의 최종 수정된 예측값은 다음과 같이 주어진다.

 

$\tilde{\boldsymbol{Y}}_n(h) = \boldsymbol{S}\hat{\boldsymbol{\beta}}_h = \boldsymbol{S}({\boldsymbol{S}}'\boldsymbol{S})^{-1} {\boldsymbol{S}}'\hat{\boldsymbol{Y}}_n(h) = \boldsymbol{SP}\hat{\boldsymbol{Y}}_n(h)$,  where $\boldsymbol{P} = ({\boldsymbol{S}}'\boldsymbol{S})^{-1} {\boldsymbol{S}}'$

 

즉, 최적 조합 방법에서 $\boldsymbol{SP}$는 선형 회귀의 Hat matrix $\boldsymbol{H}$와 같은 역할을 한다. 이때 $\boldsymbol{SP} = \boldsymbol{S}({\boldsymbol{S}}'\boldsymbol{S})^{-1} {\boldsymbol{S}}'$를 기저예측에 곱해지는 가중치라고 칭하기도 한다. 이 방식이 가지는 또 하나의 큰 장점은 기저 예측값이 unbiased일 경우 최종 수정된 예측값도 unbiasedness를 만족하게된다는 것이다.

 

 

참고 논문

G. Athanasopoulos et al. "Hierarchical forecasts for Australian domestic tourism." International Journal of Forecasting 25 (2009): 146–166

R.J. Hyndman et al. "Optimal combination forecasts for hierarchical time series." Computational Statistics & Data Analysis Volume 55, Issue 9 (2011) : 2579-2589