Taylor, S.J., & Letham, B. (2017) Forecasting at scale. PeerJ Preprints.

　Facebookが公開している時系列予測ライブラリ Prophet のテクニカルペーパーに相当する資料。仕事の都合で目を通した。
　読み終えてから気が付いたけど、同名の記事がAmerican Statisticianに載っている。たぶん中身は同じだと思う。
　Prophetについては、検索すると日本語で解説しているブログ記事がたくさんヒットするので、そういうので勉強したほうが効率がいいんだろうけど… なんというか、気分の問題です。

　イントロは端折って…
　えーと、題名でいうat scaleというのは、計算機資源のことじゃなくて、(1)時系列分析について詳しく知らないたくさんの人が使える、(2)いろんな予測問題に使える、(3)たくさんの予測ができてその評価も効率的にできる、ということだそうです。

モデル
　時系列\(y(t)\)をトレンド\(g(t)\)、季節性\(s(t)\)、休日\(h(t)\), 誤差項\(\varepsilon_t\)にわける。$$ y(t) = g(t) + s(t) + h(t) + \varepsilon_t $$ 発想はGAM(一般化加法モデル)に近く、右辺の各要素はそれぞれ線形関数だったり非線形関数だったりしうる。
　予測問題を曲線のあてはめとして捉える。つまり、データにおける時間依存性の構造を明示的に説明しようというつもりはない。このアプローチには以下の長所がある。

• 柔軟性。
• 時点の間隔が不均等でも大丈夫。
• フィッティングが速い。
• パラメータを解釈しやすい。

トレンド
　飽和成長モデルと区分線形モデルを考える。

飽和成長モデル
　まずはロジスティック成長モデルを考える: $$ g(t) = \frac{C}{1+\exp(-k(t-m))}$$ \(C\)が環境収容力、\(k\)が成長率、\(m\)はオフセットね。
　これを以下のように拡張する。

• 環境収容力\(C\)は一定でない\(C(t)\)と考える(たとえばFacebookユーザ数の時系列でいえば、トレンド項の環境収容力とはインターネットユーザ数にあたり、それ自体も成長している)。どんな関数にするか自分で考えてもよいし、外からデータを持ってきてもよい。[要は、ここは自分で決めろってことですかね]
• 成長率\(k\)も一定でなく、\(S\)個の変化点\(s_j \ (j=1,\ldots,S)\)を持つと考える[変化点の間では一定ということね]。成長率のベースラインを\(k\)、\(s_j\)における成長率の変化量を\(\delta_j\)として、$$ k(t) = k + \sum_{j:t > s_j} \delta_j$$ とする。もっとかっこよく、\(t \geq s_j\)のときに1、そうでないときに0となるダミー変数を\(a_j(t)\)とし、その長さ\(S\)のベクトルを\(a(t)\)、\(\delta_j\)のベクトルを\(\delta\)として、$$ k(t) = k + a(t)^\top \delta $$ と書いてもよい。
  　成長率を書き換えたからにはオフセットも書き換えないといけないが、下式でつじつまが合う[面倒くさいので、ここはなにも考えず写経します…]:$$ \gamma_j = \left(s_j – m – \sum_{l < j} \gamma_l \right)\left(1- \frac{k + \sum_{k < j} \delta_l}{k + \sum_{k \leq j} \delta_l} \right)$$ としておいて、 $$ g(t) = \frac{C(t)}{1+\exp(-(k+a(t)^\top \delta)(t-(m+a(t)^\top \gamma)))}$$

区分線形モデル
　こいつはシンプルです。単に\(\gamma_j = -s_j \delta_j\)として$$ g(t) = (k+a(t)^\top \delta)t + (m+a(t)^\top \gamma)$$とする。

　[はぁ？ という感じなので、具体例で考えよう。時点数を5、変化点を2,4とする。
　\(\delta\)は長さ2のベクトルで、その要素は\(\delta_1, \delta_2\)。\(a(t)\)は長さ2のベクトルで、要素\(a_1(t)\)は\(a_1(1)=0, a_1(2)=1, a_1(3)=1, a_1(4)=1, a_1(5)=1\)となり、要素\(a_2(t)\)は\(a_2(1)=0, a_2(2)=0, a_2(3)=0, a_2(4)=1, a_2(5)=1\)となる。\(\gamma\)は長さ2のベクトルで、その要素は\(\gamma_1 = -2 \delta_1, \gamma_2 = -4 \delta_2 \)。さあ深呼吸！
$$g(1) = (k + a(1)^\top \delta) 1 + (m + a(1)^\top \gamma) = k + m$$ $$g(2) = (k + a(2)^\top \delta) 2 + (m + a(2)^\top \gamma) = 2(k + \delta_1) + (m – 2 \delta_1) = 2k + m$$ $$g(3) = (k + a(3)^\top \delta) 3 + (m + a(3)^\top \gamma) = 3(k + \delta_1) + (m – 2 \delta_1) = 3k + \delta_1 + m$$ $$g(4) = (k + a(4)^\top \delta) 4 + (m + a(4)^\top \gamma) = 4(k + \delta_1 + \delta_2) + (m – 2 \delta_1 – 4 \delta_2) = 4k + 2\delta_1 + m$$ $$g(5) = (k + a(5)^\top \delta) 5 + (m + a(5)^\top \gamma) = 5(k + \delta_1 + \delta_2) + (m – 2 \delta_1 – 4 \delta_2) = 4k + 3\delta_1 + \delta_2 + m$$ なるほどー、たしかにこの式で区分線形になっている。一本目の傾きが\(k\), 二本目の傾きが\(k+\delta_1\), 三本目の傾きが\(k+\delta_1+\delta_2\)だ]

　変化点はユーザが与えてもいいし、\(\delta\)にスパースな事前分布を与えて自動選択してもよい[おっと、やっぱしベイジアンなのね]。我々がよくやるのは、変化点をたくさん用意しておいて(月一回とか)、事前分布を\(\delta_j \sim Laplace(0, \tau)\)とするというやりかた。\(\tau\)を0に近づければ区分でないロジスティック成長モデルないし線形成長モデルとなる。

　さて、変化点は未来にやってくるかもしれない。\(\tau\)に階層事前分布を与えておいて\(\tau\)の事後分布を得てもいいけど、\(\delta_j\)の推定値を使うという手もある。過去の変化点について推定値の絶対値の平均\(\lambda = \frac{1}{S} \sum_j^S |\delta_j|\)を求めておいて、未来の時点\(j\)の変化量\(\delta_j\)については、確率\((T-S)/T\)で\(\delta_j = 0\)、確率\(S/T\)で\(\delta_j \sim Laplace(0, \lambda)\)とする。「未来においても変化が同じように生じる」という強い仮定を置いているわけで、予想の不確実性を正確に評価できるわけではないけれど、すくなくともオーバーフィッティングの目安にはなる[なるほど… \(\tau\)を大きく見積もるとオーバーフィットするからね]。

季節性
　フーリエ級数を使います。
　\(P\)を規則性を持つ期間の長さとして(たとえば日次データで期間が年なら\(P=365.25\)とか)、$$ s(t) = \sum_n^N \left( a_n \cos\left(\frac{2 \pi n t}{P}\right) + b_n \sin \left(\frac{2 \pi n t}{P} \right)\right)$$ とするわけです。パラメータは\(2N\)個となる。これをベクトル\(\beta\)とし、\(cos(2 \pi n t/ P), sin(2 \pi n t / P)\)の多変量時系列を\(X(t)\)として、$$ s(t) = X(t) \beta $$ と書いてもよい。
　事前分布として\(\beta \sim Normal(0, \sigma^2)\)を与える。\(N\)は年なら10, 週なら3ぐらいがおすすめ。AICとかで選んでもよかろう。

休日・イベント
　ビジネスデータを扱うものたるや、スーパーボールとか感謝祭とかには気を付けなあきませんね。
　ある種類の休日・イベント\(i (=1,\ldots, L)\) について、過去未来の年月日の集合を\(D_i\)、パラメータを\(\kappa_i\)とし、パラメータのベクトルを\(\kappa\)として、$$ Z(t) = [ \mathbf{1}(t \in D_1), \ldots, \mathbf{1}(t \in D_L)] $$ $$ h(t) = Z(t) \kappa$$ とする。事前分布を\(\kappa \sim Normal(0, \nu^2)\)とする。[要はイベント別のフラグを立てまくろうというわけだ。この作戦って週次データだとちょっと困ったことになるよね… 去年の正月三が日と今年の正月三が日とでは週とのダブりが異なるので、同じフラグを使えるのか疑問である]

フィッティング
　上記モデルはStanで簡単に書ける。[メモは省略するけど、事前分布は k ~ normal(0, 5); m ~ normal(0, 5); epsilon ~ normal(0, 0.5), delta ~ double_exonential(0, tau); beta ~ normal(0, sigma); となっている。なお beta は季節性と休日・イベントのパラメータを一緒にしたベクトル]
　Stanコードにおけるtauとsigmaは正則化の程度を表す。過去データからCVで決めるだけの十分な過去データがない場合も多いので、デフォルト値をご用意している。
　
　分析者は往々にして、領域知識は豊かで統計学の知識は豊かでない。Prophetは、その領域知識を統計学ヌキでモデルに取り入れることができるようにしてある。取り入れる方法は次の4つである。

• 収容力。
• 変化点。直接指定できます。
• 休日と季節性。
• スムージングパラメータ \(\tau, \sigma, \nu\)。

Prophetは豊富な視覚化によって分析者の判断を支援いたします。

　我々の目指すところの、視覚化→モデリング→予想の評価→問題の発見というループは、いわゆる統計的予想とjudgmental予想のいいとこどりである。Hadley Wickhamさんいうところのtransform-visualizeループに近い。分析者たちよ、胸にあるものいつか見えなくなるもの、それはそばにいることをいつも思い出して。統計学を超えて行け、過去経験を超えて行け。[すいませんこんなことは書いてありませんでした]

予想の自動的評価
　最後に、予想の成績評価を自動的に行う手順について。
　まずはベースライン予想を作っておくことをお勧めします。最終の値を予想にするとか、標本平均を予想にするとか、ARIMAモデルを使うとか。

　予想ホライゾンを\(H\)とする(たとえば日次で向こう一年を予想するなら\(H=365\))。時点\(T\)までの過去情報を使った時点\(t\)の予想を\(\hat{y}(t|T)\)とする。なんらかの距離\(d(y, y’)\)が定義できるとする(なんでもいいけどMAPEがわかりやすいのでお勧め)。
　\(h \in (0, h]\)における予想のの正確性を$$ \phi(T, h) = d(\hat{y}(T+h|T), y(T+h)) $$ とする。いろんな\(h\)についての誤差の期待値のモデル\(\xi(h) =E[\phi(T, h)] \)について考える。この関数は局所平滑で、\(h\)とともに弱増大だと仮定する。

　さて、具体的にどうするかというと… 予想ホライゾンが過去データに入っているようなカットオフ点をK個設定し、それぞれのカットオフ点までのデータを使って予想を得る。我々はこれをsimulated histrocal forecasts(SHFs)と呼んでいる。いわゆるrolling-origin評価と似ているが、可能なカットオフ点をすべて使うわけではない。目安としては\(H/2\)ごとに予想を得ればよい。で、横軸に\(h\), 縦軸にMAPEをとった散布図を描いて、(たとえば)LOESSでフィッティングする。
　[カットオフ点の決め方はよくわからんが、まあどうでもいいや。とにかく、真面目にrolling-origin評価をしなくてもいいじゃんという主旨である。さすがはIT企業、著者らは暗黙のうちに長めの日次時系列を想定しているのだ。伝統的なマーケターが思い浮かべる月次売上時系列ではなかなかこうはいくまい]
　大きな予想誤差を素早く見つけ出すことが大事である。そのためには、(1)ベースラインと比べる、(2)すべての手法において大きな誤差が出る場所を探す(外れ値と考えられる)、(3)SHF誤差が急に大きくなるカットオフ点を探す(変化点を追加すべきと考えられる)、などがおすすめ。
　云々。
—–
　なるほどね… 時系列データについてその生成モデルを組むつもりはさらさらなくて、単に区分成長曲線(or 区分線形な曲線)＋フーリエ級数＋ダミーフラグにフィッティングしようとしているわけだ。これはこれでひとつの考え方ではある。

　疑問その１。説明変数の時系列を使うという話は全然出てこなかったけど、それはやらないのかな？
　すでに\(s(t), h(t)\)で外生変数を使ってるから、\(y(t)\)に即時的に効き永続性を持たない説明変数であれば、追加するのは容易だろう。飲料売上に対する気温とか。
　しかし、ものの本によれば、ブランド売上時系列の68%が単位根を持つのだそうだ。つまり、外生的ショックが永続的な影響を持っているわけである。たとえば新規ブランドの導入期などには、きっと売上時系列を下から叩く小さなショック(広告やプロモーション)と上から叩く小さなショック(競合の広告やプロモーション)が日々こつこつと続き、それぞれのショックが長期的な影響を持つのだろう。そういうのって、このモデルではどう表現すればいいんだろうか？ \(g(t)\)の分母に外生変数を組み込む方向で拡張したとすると、今度は、説明変数を\(g(t)\)の分母に入れるべきか\(h(t)\)に入れるべきかをユーザが判断しなければならなくなり、統計的知識が必要になり… ううむ、主旨に反するぞ。

　疑問その２。進化的な市場の売上時系列がそうであるように、\(y(t)\)が単位根過程であると分析者が知っている場合(つまり、未知の外生的ショックがあったとしてそれは永続的な影響を持つだろうと思っている場合)、それでも区分成長モデルをあてはめるというのは良い方法なのか。あ、ひょっとして差分系列に区分線形モデルをあてはめればいいの？　いや、おそらく「これは進化的な市場だから云々」という考え方自体が生成モデル的であって、著者らの発想とは異なるのだろう。

　追記: この解説の公開後、prophetはいろいろ機能追加されている模様。それらについての日本語の解説も山のようにある模様。だんだん勉強する気が失せていく… まあいいや、あとは公式ドキュメントをみるのが早かろう。