Hyndman, R.J., Billah, B. (2003) Unmasking the Theta method. International Journal of Forecasting, 19, 287-290.
仕事の都合で読んだ奴。タイトルのとおり、時系列予測の手法のひとつTheta法についての解説。たったの4p。
第一著者のHyndmanさんはいわずとしれた有名人で、Rのforecastパッケージの中の人であらせられる。
いわく。
Theta法はAssimakopoulos & Nikolopoulos (2000, Int.J.Forecasting)が提案した時系列予測手法で、時系列予測のコンペティションM3 Competitionで謎の好成績を叩き出し注目を集めた。しかし提案者による説明は複雑であり、おそらくは混乱を招く代物だ[←最初の段落でいきなりディスるスタイル。さすがはHyndmanさん]。
というわけで解説をお送りしよう。話を先取りすると、実はTheta法はドリフト付きの単純指数平滑化だ。
観察時系列を \( \{ X_1, \ldots, X_n \} \) とする。ここから時系列 \( \{ Y”_{1,\theta}, \ldots, X”_{n, \theta} \} \) を生成する。ただし、\( t = 3, \ldots n \)について $$ Y”_{t,\theta} = \theta X”_t $$ であり、ダブルプライムは2階階差を表す[ なるほど？ つまり \(Y_{t,\theta}\) はもとの時系列の差分系列の差分系列の\( \theta \)倍の和分系列の和分系列ってことね]。A & Nはこれを「シータ線」と呼んでいる。
この二次差分方程式の解は、\( t=1 \ldots, n \)について $$ Y_{t, \theta} = a_\theta + b_\theta (t-1) + \theta X_t$$ ただし \( a_\theta, b_\theta \)は定数。つまり\( Y_{t,\theta} \) は \(X_t\) に線型トレンドをつけたものだ。[ \( t=1,2 \) についてどう導出するのかちょっとぴんとこないんだけど、まあ俺の学力不足だろう]
A & Nは、2つの時系列の差の平方和を最小化するような \( Y_{1,\theta} \) と \( Y_{2,\theta} – Y_{1,\theta} \)をみつけている。これはどういうことかというと、$$ \sum_{i=1}^t [X_t – Y_{t,\theta}] = \sum_{i=1}^t [(1-\theta) X_t – a_\theta – b_\theta (t-1) ]^2 $$ を最小化する \( a_\theta, b_\theta \)を見つけるという話である。これ、\( (1-\theta) X_t \) を \( t-1 \) に回帰してるよね。ってことは、解は
$$ \hat{b}_\theta = \frac{6 (1-\theta)}{n^2-1} \left( \frac{2}{n} \sum_{t=1}^n t X_t – (n+1) \bar{X} \right) $$ $$ \hat{a}_\theta = (1-\theta) \bar{X} – \hat{b}_\theta (n-1) /2 $$ だよね。これが A&N が示したことなのだ。まあ彼らはこの導出に2頁以上使ってるけどね。[←いちいち辛口ですね]
以下の点にご注目。
- \( \theta = 0 \)のとき、\( \hat{a}_\theta, \hat{b}_\theta \) は \( \{ X_t \} \) にあてはめた線型トレンドのパラメータになっている。
- \( \hat{b}_\theta / (1-\theta) \) は\( \theta \) と独立である。
- \( \bar{Y}_\theta = \hat{a}_\theta + \hat{b}_\theta (n-1)/2 + \theta \bar{X} = \bar{X} \) である。
- \( (1/2) [Y_{t,\theta} + Y_{t,2-\theta} ] = X_t \) となる。
さて。A&Nが提案した予測手法は、異なる \( \theta \) について \( Y_{t,\theta} \) を求めこれを加重平均するというものだ。A&NがM3 Competitionで使ったのは \(\theta = 0\) と \(\theta = 2\) である。
\(\theta = 0\)のときの予測とは、ただの線形補外だ。\( h \)期先予測は $$ \hat{Y}_{n,0}(h) = \hat{a}_0 + \hat{b}_0(n + h – 1)$$ \(\theta = 2\)のときの予測とは、\( \{ Y_{t,2} \} \) の単純指数平滑化だ。つまり$$ \hat{Y}_{n,2}(h) = \alpha \sum_{i=0}^{n-1} (1-\alpha)^t Y_{n-t, 2} + (1-\alpha)^n Y_{1,2} $$ [あっそうなの？ もともとのA&Nの提案からして, 所与の\(\theta\)における予測はただの単純指数平滑化なの？ それともA&Nの提案を読み解くと単純指数平滑化に帰着するってこと？ いや、まあ、どっちでもいいっすけどね]
我々の定式化から、いろんなことが簡単に導ける。
\( \{X_t\} \)の単純指数平滑化による予測を \(\tilde{X}_n(h) \) としよう。[…中略…] \(n\)が大きいとき、Theta法による予測は $$ \hat{X}_n(h) = \tilde(X)_n (h) + \frac{1}{2} \hat{b}_0 (h-1+1/\alpha) $$ となる。つまり、単純指数平滑化による予測に、もとの\( \{X_t\} \)にあてはめたトレンド線の傾きを半分にしたトレンド線をくっつけたものである。
[←へえええ。つまりM3 Competitionで起きたことは、単純指数平滑化よりも、それにちょっと過去トレンドを追加した奴が勝った、ってことなのね。なるほど？]
ついでに、この予測の背後にあるstochasticなモデルをお示ししよう。
\(X_1 = l_1)\)とし、\( t=2 \)以降について$$ X_t = l_{t-1} + b + \epsilon_t $$ $$ l_t = l_{t-1} + b + \alpha \epsilon_t $$ とする。ただし \( \{\epsilon_t\} \)は平均0, 分散\( \sigma^2 \)の正規ホワイトノイズ。
[おおっと、状態空間表現のおでましだ。状態方程式はドリフトつきランダムウォークで、観察方程式は状態に状態方程式のドリフト項とノイズをのっけたものなわけね]
この\( X_t \)はドリフト付き単純指数平滑化による予測と等価である。[えっそうなの?! もう少し親切に教えてよ、先生…]
さらにいうと、これはHolt法で傾きの平滑化パラメータをゼロにしたのと同じである。さらにさらに、これは$$ X_t = X_{t-1} + b + (\alpha – 1)\epsilon_{t-1} + \epsilon $$ つまりドリフト付きのARIMA(0,1,1)である。[へええ… はいはい信じますよ先生]
さて、ここで \(b = \hat{b}_0 / 2 \)とするとTheta法と同じ点予測が得られる。またこの状態空間モデルに基づけば […中略…] \(h\)期先予測の95%信頼区間は$$ \hat{X}_n (h) \pm 1.96 \sigma \sqrt{(h-1) \alpha^2 +1 } $$ である。ただし、(たいていの時系列予測と同じく) 推定誤差は入ってないことに注意。
こうしてみると、\(b = \hat{b}_0 / 2 \)じゃない\(b\)を使えば、M3 Competitionでもっと勝てたのではないか。
M3 Competitionのデータを使い、次の3通りを比べてみた。
- A&Nが出した予測値
- 我々が実装したTheta法で出した予測値
- 上の状態空間モデルで出した予測値。\(b\)も最適化するわけである。ついでに(上の定式化とは違って)\(l_0\)も未知パラメータとして最適化した。
初期化の違いだか推定量の違いだかで、1と2はわずかに異なった。sMAPEで勝ったのは3であった。[\(h\)別に表になっているけど、そんなに大差で勝てているわけでもない。レベルシフトが大きい時系列では負けた旨の注記がある]
… なるほどです。勉強になりましたです。
Hyndmanさんは私のなかの「難しい話を分かりやすく説明する人」ランキングのかなり上位にいる人なんだけど、この論文のせいでさらにランクアップしました。誉めてつかわすぞ。[態度がでかい]