線形ガウス状態空間モデルとPOSデータ分析（前半）

■ 前書き

前回は、状態空間モデルの基本形と、「説明変数の効果が時間とともに変わる」状況を潜在状態として表現できることを解説した。本稿ではその発展として前半と後半に分け、読者がやや複雑な線形ガウス状態空間モデルをドメイン知識に合わせてカスタムに設計し、Pythonで実装できる形に落とし込めるようにすることを目標にする。

■ 記法

本稿で扱う線形ガウス状態空間モデルは、以降つぎの記法で統一する。

$\mathbf{x}_t$は各時点の状態、$y_t$は各時点の観測値、$\mathbf{G}_t$は状態遷移行列、$\mathbf{F}_t$は観測行列である。(本稿では、$y_t$がスカラーであるため$\mathbf{F}_t$は行ベクトルになる。)また、$\boldsymbol{w}_t$と$v_t$は互いに独立な正規分布に従い、さらにともに時点間でも独立である。さらに、任意の$t$と$s$について、$\boldsymbol{w}_t$と$v_s$は独立であり、初期状態$\mathbf{x}_0$もこれらのノイズと独立と仮定する。$\mathbf{m}_0,$$\mathbf{P}_0$は初期状態の事前分布の平均ベクトルと共分散行列を表す。

■ 線形ガウス状態空間モデルの利点

線形・ガウスという仮定は一見シンプルで強い制約であるように思われるかもしれないが、実務では強い利点がある。線形・ガウスという仮定の下ではカルマンフィルタにより、潜在状態の推定と尤度計算を複雑な近似計算を必要とせず数値的に扱えるため、モデルの試行錯誤が速い。

さらに、状態空間モデルの強みは、現場の仮説（例：レベル、成長率、曜日パターン、外生要因の効き方）を、状態ベクトル$\mathbf{x}_t$の各要素として定式化し、必要に応じて次元を拡張するだけで、1つのモデルに統合できる点にある。例えば、

• – レベル（全体水準）は競合や露出で少しずつ動く
• – 成長率（傾き）は加速したり鈍化したりする
• – 曜日パターンはあるが、流通や生活様式で少しずつ形が変わる
• – 検索（Googleトレンド）が効いても、季節によって効き方が弱まる時期がある

といった「関係者の直感」を、パラメータの中身として設計できる。

実務で状態空間モデルを活かす時の難所は、カルマンフィルタのアルゴリズムよりも、むしろ 「ドメイン知識をどう数式（行列）に落とすか」にあることが多い。この点が曖昧だと、ライブラリを利用しても実装すること自体が難しい上、実装できたとしても「それっぽく動くが解釈できない／どの成分が何を指すかが一意に決まらない」などの問題に直面してしまうことがある。しかし、ドメイン知識の数式への落とし込みは、状態空間モデルをモデリングするにあたっての非常に面白いところでもあると筆者は考える。

そこで前半では、実装の前提となる「数式への落とし込み」に集中し、最終的に Python で実装できるための準備を行う。具体的には、まず レベルと成長率を状態に入れ、次に週次成分(曜日効果)を組み込むところまでを、数式として明示する。後半では、この仕様をPythonで実装し、日経POSデータを題材に、推定された潜在状態(レベル・成長率・曜日効果など)がどのように解釈できるかを確認する。

■ ドメイン知識を数式に落とし込む

以下では、ある「状態」の$t$時点と$t-1$時点での推移式を示した後、そこを部分的にモデリングした際の状態推移行列$\mathbf{G}_t$と状態ノイズの共分散行列$\mathbf{W}_t$を示す、という形で少しずつモデルを組み立てていく。状態推移行列については、この推移式の右辺の係数を表していることを覚えておくと、理解がしやすいだろう。

まず、「全体の水準が少しずつ動く」という言い方は、数式に落とすと以下のように表せる。

ここで$\mu_t$は$t$時点での水準を表し、$\eta^{(\mu)}_t$は「水準がどれくらい揺れるか」を表すノイズである。これは、「今日の水準は昨日の水準を、揺れを伴いながら引き継ぐ」という最も素朴な仮定だ。

ただ、現場のデータを眺めると「水準の揺れ」だけでは表現しづらい局面がある。たとえば、増分がしばらく同じ方向に続いたり、成長率が加速・鈍化していくようなときだ。そういう場合は、傾き$\delta_t$を状態として持ち、それを水準$\mu_t$に反映したくなる。これは、ローカル線形トレンドモデル（local linear trend model）と呼ばれる。このモデルは次式で表せる。

ここで$\eta^{(\delta)}_t$は、「傾きがどれくらい揺れるか」を表すノイズである。レベル$\mu_t$の更新に前時点での傾き$\delta_{t-1}$を加えることで、「日にちが近い時の増分の値同士は比較的近くなる（増分が短期で激しく上下に振れにくい）」という解釈を数式で反映している。

この「水準」と「傾き」までを状態に入れると、$\mathbf{x}_t^{(\mathrm{trend})} = (\mu_t, \delta_t)^\top$となる。このときの状態遷移行列$\mathbf{G}_{\text{trend}}$と状態ノイズの共分散行列$\mathbf{W}_{\text{trend}}$は、以下のようになる。状態遷移行列の1行目がこのモデルの$\mu_{t}$の式の右辺の係数、2行目が$\delta_{t}$の式の右辺の係数に対応する。

さらに、次に「曜日で売上が決まっている」というドメイン知識を盛り込みたい。直感的には「月曜は弱い／週末は強い」といった、週次（7日周期）の繰り返し構造だ。しかし、曜日による影響も固定ではなく、流通オペレーションの変化などによって時によって変わりうることを反映させたい。そこで、曜日パターンが完全に固定という仮定と、日々自由に変わるという仮定の中間として、「週次性は保ちつつ形は緩やかに変化し得る」モデルを採用する。

まず直感のため、各曜日の効果を $\gamma_{t,1},\ldots,\gamma_{t,7}$と書く（$\gamma_{t,1},\ldots,\gamma_{t,7}$は時点$t$における「$j$番目の曜日」の効果を表す）。観測は「レベル＋その日の曜日効果」という形で入るため、概念的には

（$d(t)$は時点 $t$の曜日を表す添字）と書ける。

ただし、曜日効果を7個すべて自由に持つと、レベル $\mu_t$と曜日効果の「役割分担」が曖昧になる。実際、すべての曜日効果に同じ定数 $c$を足し、同じ $c$をレベルから引いても

となって、同じ観測を生成できてしまう。つまり、観測 $y_t$だけからは「それはレベルなのか、曜日効果なのか」を一意に区別できない（同定不能）状態になる。
そこで、週次成分には

という制約を課し、平均（全体水準）はレベル$\mu_t$が表し、曜日効果はそこからの偏差だけを表すようにする。この制約の下では自由度が1つ減るため、週次成分は7−1=6次元で表現できる。

次に、この「自由度 6」の週次成分を状態として持つための具体的な表現を導入する。以降、週次成分（観測に入るその日の曜日効果）を$s_t$と書き、実装上は $s_t=\gamma_{t,1}$となるように状態を定義する。つまり、1日進むごとに状態をシフトさせることで曜日が巡回する。そのために状態ベクトルとして

を導入し、7個目は和ゼロ制約により

と一意に定まる（これが「7−1=6 個で表す」の意味である）。

1日進むと状態はシフトし（昨日の効果が2番目に回る）、さらに「暗黙の7個目」を新しい先頭として計算する。こうすれば、7日間で(ノイズを除いて)状態が一周して曜日効果である$s_t$に入る。具体的には、次のように書ける。

ここで入れた $\eta^{(s)}_t$は「曜日パターンがどれくらいゆっくり形を変えるか」を制御するノイズである。もし$\eta^{(s)}_t$の分散である$\sigma_s^2$を0にすると、曜日パターンは固定（完全に同じ形が繰り返される）になり、分散を大きくすると、曜日の山谷が時間とともにより揺れやすくなる。この推移は、状態遷移行列のブロックとして次のように書ける。1行目が$\gamma_{t,1}$の式、2行目が$\gamma_{t,2}$の式…にそれぞれ対応する。

さらに状態ノイズの共分散行列については、ノイズは状態$\gamma_{t,1}$のみに加わるため以下のようになる。

加えて最後に、以前のコラムで示したような検索需要などの外生変数$z_t$が売上に影響し、その効き方が季節や市場環境によって刻々と変わることもモデリングに入れたい場合を考える。これは、時間で変化する回帰係数$\beta_{t}$そのものを状態として持つ場合としてモデリングできる。

この$\beta_{t}$の推移を推移行列$\mathbf{G}_t$の要素に加えるためには、新たなブロック要素としてスカラー値1を加えれば良い。$\mathbf{W}_t$については、新たなブロック要素としてスカラ値$\sigma_\beta^2$を右下に加えれば良い。

ここまでの「レベル＋傾き」の状態$(\mu_t,\delta_t)$と、この週次状態$\boldsymbol{\gamma}_t$、さらには$\beta_{t}$を結合して、状態ベクトル$\mathbf{x}_t$を

と置けば、式全体を組み合わせた状態遷移行列$\mathbf{G}_t$は以下のブロック対角の形になる。

このように、状態遷移行列の対角上に状態要素がどのように推移するかのドメイン知識を組み合わせることで、様々なドメイン知識を組み入れてモデリングできる。

さらに状態ノイズの共分散行列$\mathbf{W}_t$についても、以下のようなブロック対角の形になる。

そして売上の観測値は、レベル＋曜日効果＋外生変数の影響＋観測ノイズとして

と書ける。

観測行列$\mathbf{F}_t$は状態$\mathbf{x}_t$の中で観測値$y_t$に足される部分に対応する次元以外の成分を0にする。すなわち、

となる。($\beta_t$に対応する部分は外生変数との積をとる必要があるため、$z_t$を置いている)

ここまでで、数式への落とし込み方の例について説明した。後半ではこれをもとに、どのように実装するかと、実際の分析例について説明したい。