« 読了：「コミュニティ 安全と自由の戦場」 | メイン | 読了：Wilbur & Farris (2014) 配荷率増大が売上に与える効果はシェアが大きい製品において大きい、SKUレベルでみてもやっぱりそうだった »
2019年10月 3日 (木)
Hyndman, R.J., Athanasopoulos, G, Shang, H.L. (2013) hts: An R Package for Forecasting Hierarchical or Grouped Time Series.
階層時系列データを分析するためのRパッケージhtsのvignette。中の人は誰あろう、forecastパッケージで知られるHyndmanさんである。
ここで階層的時系列というのは、たとえば国レベルの死亡者数時系列があり、これを地域別の時系列とか性別の時系列とに分けることができる、というようなことを指している。
説明例として2レベルの階層を考える。全体(レベル0)の時系列を$y_t$、レベル1にAとBがあって、その時系列を$y_{A,t}, y_{B,t}$とする。Aはレベル2に子AA, AB, ACを持ち、Bは子BA, BBを持ち、それらの時系列を$y_{AA,t}, \ldots, y_{BB, t}$とする。
時系列変数のベクトルを考えよう。$\mathbf{y}_t$をこの8本の時系列変数のベクトルとする。$\mathbf{y}_{K,t}$をレベル2の5本の時系列変数のベクトルとする。すると、
$\mathbf{y}_t = \mathbf{S} \mathbf{y}_{K,t}$
と書ける。ここで$\mathbf{S}$は$8 \times 5$のダミー行列で、1行目($y_t$に対応)はすべて1となり, 2行目($y_{A,t}$)は左から順に1, 1, 1, 0, 0となり、下部の5行だけ取り出すと単位行列になる。なるほどね？
で、ここでの問題は、8本の時系列すべてに対して関心を持っているという点である。
以下では、なんらかの方法によって得た、$j(=1,\ldots,8$)本目の時系列の$h$期先予測を$\hat{y}_{j,h}$と書く(ほんとは$h$じゃなくて$T+h|T$という気持ちだが略記する)。これをベース予測と呼ぼう。
ベース予測は、階層構造に照らしてつじつまが合うとは限らない。つじつまがあうように修正した予測を$\tilde{y}_{j,h}$と書き、これを修正予測と呼ぼう。
さて、階層的時系列の予測には、ボトムアップ法、トップダウン法、ミドルアウト法がある。さらに我々の考えた最適予測結合法も紹介しよう。
まずはボトムアップ法。
上の例でいうと、まずレベル2で$\hat{y}_{AA,h}, \ldots, \hat{y}_{BB,h}$を得る。これがベース予測。つまり、レベル0-1の3本についてはまだ手に入らない。
で、これを積み上げてレベル1の$\tilde{y}_{A,h}, \tilde{y}_{B,h}$を得て、さらにこれを合計してレベル0の$\tilde{y}_{h}$とする。これが修正予測。つまり、レベル2の5本については修正しない。
良い点は、累積において情報が失われないという点。悪い点は、レベル2のデータはノイズが多いことが多いので、モデリングが大変だという点。
お次はトップダウン法。
最初にレベル0で$\hat{y}_h$を得る。で、これらを割合$p_1, \ldots, p_{m_2}$によってレベル2に配分する。最後にレベル2の時系列を積み上げてレベル1の時系列を得る。
この方法がうまくいくかどうかは次の3つの要因で決まる。(1)レベル1の予測の正確性。(2)配分割合の正確性。(3)ベース予測の正確性。[← ??? (1)と(3)は同じことじゃないの？]
配分割合の決め方にはいろいろある。紹介しよう。
その1, 過去の割合の平均
$p_j = \frac{1}{T} \sum_t^T \frac{y_{j,t}}{y_t}$
を使う。
その2, 過去の平均の割合
$p_j = \sum_t^T \frac{y_{j,t}}{T} / \sum_t^T \frac{y_t}{T}$
を使う。
その3、割合を予測する。まず、すべての時系列についての$h$期先予測を手に入れる。で、レベル1から順に、$h$期先予測に占めるそのノードの割合を求める。これを一番下のレベルまで繰り返す。[数式はめんどくさいので略]
配分割合をどうやって決めるかはともかくとして、トップダウン法の弱点は、仮にベース予測が不偏であったとしても、修正予測は不偏じゃない、という点。
ミドルアウト法。どこかのレベルに注目し、そのレベルより上に向かってはボトムアップ法、そのレベルより下に向かってはトップダウン法を使う。
お待たせしました、我々が考えた最適予測結合法です。
この方法では、まずは全レベルで個別に$h$期先予測値を得る。そのベクトルを$\mathbf{\hat{y}}_h$としよう。で、
$\mathbf{\hat{y}}_h = \mathbf{S} \mathbf{\beta}_h + \mathbf{\epsilon}_h$
というモデルを考える。$\mathbf{\epsilon}_h$は$h$期先予測誤差ではなくて、単にこの回帰式における誤差である。その共分散行列を$\mathbf{\Sigma}_h$としよう。
一番下のレベルの$h$期先予測誤差を$\mathbf{\epsilon}_{K,h}$としよう。$\mathbf{\epsilon}_h \approx \mathbf{S} \mathbf{\epsilon}_{K,h}$と仮定すれば、$\mathbf{\beta}_h$の線形不偏推定量は
$\mathbf{\tilde{\beta}}_h = \mathbf{S (S' S)^{-1} S' \hat{y}_h}$
となるので、これを求めて、上のモデルに放り込めば、修正予測$\mathbf{\tilde{y}}_h$が手に入る。[なるほどー]
ここからはhtsパッケージの紹介。
関数hts()は、階層構造の指定と、一番下のレベルの時系列データを食って、全レベルの時系列を持つオブジェクトを返す。階層がなくて単にグループがあるだけならばgts()でもよい。
関数forecast()はこのオブジェクトを食って$h$期先予測を返す。予測モデルとしてはarima, ets, ランダムウォークをご用意している。で、ボトムアップ予測とか、トップダウン予測とか、いろいろできるらしい。
forecastパッケージのforecast()とは異なり、どうやら説明変数は使えないようだ。自分の仕事のなかでは、利用場面がかなり限られてくる...がっくり...
論文：データ解析(2018-) - 読了：Hyndman et al. (2013) 階層時系列データ分析のためのパッケージ hts