ill-identified diary

5 million yen salary expected.

[R] [bsts, dlm, KFAS] マーケティングの状態空間モデリング

概要

  1. 岩波 DS Vol. 6 での佐藤忠彦 (2017, 状態空間モデルのマーケティングへの応用)の記事でなされた小売業の売上量のモデリングを R で再現してみる.
  2. dlmKFAS, そしてbstsパッケージを利用して, それぞれでプログラムを書いてみる.
  3. 最近はゆるめの読み物ばかり書いてきたので, 硬派にやっていく. よってそれなりにハードな内容である. 量はだいたい7ページ分くらい.

初めに

最近は技術的なことをあまり書かずに抽象的な話ばかりやってきたので, 久々に戻そうと思う. そこで, 状態空間モデルを用いてマーケティングにおける売上数予測モデルの構築の例をやってみる.

今回は読者が状態空間モデルの話を知っている前提で書いている. 状態空間モデルを解説した記事には, 例えば自分が昔書いたものがある.

もちろん他にも探せばいくらでもある.

では具体的になにをやるかというと, 佐藤忠彦 (2017, 状態空間モデルのマーケティングへの応用) などで紹介されている動的市場反応モデルを使う. これは名前の通り, 消費者の意思決定のルールが時間とともに変化するモデルである. このモデルは状態空間モデルとして扱えるので, dlm, KFAS, それから最近登場した bsts パッケージを使って実際のデータに当てはめることで, パラメータの推定や予測が可能となる.

この記事では, まず動的市場反応モデルがどのような理論に基づくものかについて説明する. 次に, R の3つのパッケージ, dlm, KFAS, bsts をそれぞれ用いた場合の計算プログラムを, 簡単な解説とともに掲載する.

動的市場反応モデル

モデルは 佐藤忠彦 (2017, 状態空間モデルのマーケティングへの応用), 佐藤 & 樋口 (2014, ビッグデータ時代のマーケティングベイジアンモデリングの活用) でその詳細が述べられているので, ここでは最低限の設定だけを書いておく. 小売店での売上点数の時系列予測モデルについて言及されている. 消費者が商品を買うかどうかの意思決定として最も単純なのは価格であるが, 価格と売上数の回帰分析だけでは有意義なモデリングとはならない. 価格以外の, 消費者の潜在的な購買決定の構造をモデルに取り入れるのがこのモデリングの目的となる. まず, より基本的な「市場反応モデル」について考えてみる. これは以下のように与えられる*1.

\begin{align}
\ln y_{t}= & c+(\alpha+\beta e_{t})\ln p_{t}+\varepsilon_{t}\label{eq:model1}\end{align}

y_{t} は売上点数で, e_{t} は商品の山積み陳列の有無を表すダミー変数, p_{t} は商品の販売価格である. そして \varepsilon_{t} はシステムノイズである. つまり, このモデルは販売価格と売上量の関係が, 山積み陳列を行うことで意思決定に変化が生じるかを見ている. これは最も基本的なモデルであるが,

実際のデータを観察した結果, パラメータが観察期間を通して一定であるという仮定が現実とそぐわないから, さらに, 定数項と係数が時間によって変化するバージョン, いわゆる時変 (time-varying) 係数モデル*2,


\begin{align}
\ln y_{t}= & c_{t}+\mathit{ar}_{t}+(\alpha_{t}+\beta_{t}e_{t})\ln p_{t}+\varepsilon_{t}\label{eq:model2}\end{align}

もモデル候補として考慮する. c_{t} はローカルトレンド項で, \mathit{ar}_{t} は2階の自己回帰項である. よって, 遷移方程式 (状態方程式とも呼ぶ) は以下のようになる.

\begin{align}
c_{t+1}= & c_{t}+v_{1,t},\\
\mathit{ar}_{t+1}= & \phi_{1}\mathit{ar}_{t}+\phi_{2}\mathit{ar}_{t-1}+v_{2,t}\\
\alpha_{t+1}= & \alpha_{t}+v_{3,t}\\
\beta_{t+1}= & \beta_{t}+v_{4,t}\end{align}

上記の v_{\cdot,t} は全て, 平均ゼロのホワイトノイズ項である. パラメータが時間によって変化していることが示唆されるが, 実際にこの変化がどういうメカニズムで発生しているかは説明できない. そこで, 平滑化事前情報アプローチの考え方に従い, 上記のような設定とした. この設定は, パラメータが時間の経過とともに少しづつ「滑らかに」, かつランダムに変化することを意味する.

さらに, 消費者行動理論にある内部参照価格 (reference price, 単に参照価格とも) の理論に基づき, モデルをもう1段階拡張する. 参照価格の理論とは, 消費者が絶対的な価格の大小に反応して購買を決めているのではなく, 本人が適正と考えた価格より安いか高いかで購買を決めているという理論である. この判断基準となる価格が参照価格である. そして, この参照価格もまた一定ではなく, 日々の買い物での経験を通して変化していく. いわば, 経験的に形成された「相場観」のようなものである.

参照価格を取り入れたモデルは動的市場反応モデルと呼ばれ, 以下のよう変更される.


\begin{align}
\ln y_{t}= & c_{t}+\mathit{ar}_{t}+(\alpha_{1,t}+\beta_{1,t}e_{t})z_{1,t}+(\alpha_{2,t}+\beta_{2,t}e_{t})z_{2,t}+\varepsilon_{t}\label{eq:model3}\end{align}

2017/9/6: ynakahashi さんの指摘により式の書き間違い修正

\begin{align}
Z_{1,t}:= & \begin{cases}
\frac{p_{t}-\mathit{RP}_{t}}{100} & \text{if }p_{t}>=\mathit{RP}_{t}\\
0 & \text{otherwise }
\end{cases},\\
Z_{2,t}:= & \begin{cases}
\frac{\mathit{RP}_{t}-p_{t}}{100} & \text{if } \mathit{RP}_{t} \gt p_{t} \\
0 & \text{otherwise}
\end{cases}\\
\end{align}

\mathit{RP}_{t} は参照価格を表すので, これは (販売価格) >= (参照価格) の 場合と大小逆の場合とでスイッチするモデルである. 時変パラメータが \alpha_{t},\beta_{t} から \alpha_{1,t},\beta_{1,t}\alpha_{2,t},\beta_{2,t} に増えているが, それぞれ前と同様に平滑化事前情報アプローチに従った変化をするものとして扱う. なお, なぜ大小関係が異なる場合でそれぞれ係数を別にするのかというと, (販売価格) > (参照価格), つまり消費者が損をすると感じた場合と, 逆に得をしたと感じる場合とで, 反応の大きさが異なるというプロスペクト理論を取り入れているからである.Z_{1,t} が損の大きさで, Z_{2,t} が得の大きさになる.

最後に, データとして得ることの出来ない 参照価格 \mathit{RP}_{t} を考える. 参照価格が日々の経験によって形成されるという仮説から, 0\leq a\leq1 のパラメータを用いて, 1期前の参照価格と実際の価格に応じて決定されるような遷移方程式 \begin{align}
\mathit{RP}_{t+1}= & a\mathit{RP}_{t}+(1-a)p_{t}\end{align} を仮定する. ただし, 参照価格とは本来, 個人レベルの意思決定において想定されるもので, 今回のモデルは集計されたデータであり, 理論的に厳密でない, 近似的な設定であることに注意する.

要約すると, 今回の動的市場反応モデルは, 以下のような仮説に基づいていることになる.

  1. 消費者の価格に対する反応は一定ではなく, 日々変化する.

  2. 消費者は経験から参照価格を形成し, 参照価格と販売価格の差から購買を決める.

  3. 消費者は, 参照価格と販売価格の差を見て損をしていると判断した場合と得をしたと判断した場合とで, 反応の大きさが同じではない (損を過大評価する, または得を過大評価する).

実践

今回は岩波データサイエンス Vol. 6 の 佐藤忠彦 (2017, 状態空間モデルのマーケティングへの応用) のデータを利用する.
データは以下のサポートサイトでダウンロードできる. ただし, 原典で使われたデータとは少し違うらしいので注意.

https://sites.google.com/site/iwanamidatascience/vol6/time-series

計算方法

動的市場反応モデルの状態空間モデルとしての表現を明示しておくと, 以下のようになる. \boldsymbol{x}_{t} が状態変数ベクトル, \mathbf{Z}, \mathbf{Q} が遷移方程式の構造パラメータ, \mathbf{H}_{t} が観測方程式のデザイン行列である (見やすさのため, 自己回帰項の状態変数を最後に表記している).

\begin{align}
\boldsymbol{x}_{t}:=\begin{bmatrix}c_{t}\\
\alpha_{1,t}\\
\beta_{1,t}\\
\alpha_{2,t}\\
\beta_{2,t}\\
\mathit{ar}_{t}\\
\mathit{ar}_{t-1}
\end{bmatrix}, & \boldsymbol{Z}:=\begin{bmatrix}1 & 0 &  &  &  &  & \mathbf{O}\\
0 & 1 & 0\\
 & 0 & 1 & 0\\
 &  & 0 & 1 & 0\\
 &  &  & 0 & 1 & 0\\
 &  &  &  & 0 & \phi_{1} & \phi_{2}\\
\mathbf{O} &  &  &  & 0 & 1 & 0
\end{bmatrix},\\
\boldsymbol{Q}:= & \begin{bmatrix}1 &  &  &  &  &  & \mathbf{O}\\
 & 1\\
 &  & 1\\
 &  &  & 1\\
 &  &  &  & 1\\
 &  &  &  &  & 1\\
\mathbf{O} &  &  &  &  &  & 0
\end{bmatrix}\\
\mathbf{H}_{t}= & \begin{bmatrix}1 & Z_{1,t} & E_{t}Z_{1,t} & Z_{2,t} & E_{t}Z_{2,t} & 1 & 0\end{bmatrix}\end{align}

\begin{align}
\boldsymbol{x}_{t+1}=\mathbf{Z}\boldsymbol{x}_{t}+\boldsymbol{Q}\boldsymbol{v}_{t}, & \boldsymbol{v}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{V})\\
\ln y_{t}=\mathbf{H}_{t}\boldsymbol{x}_{t}+\varepsilon_{t}, & \varepsilon_{t}\sim\mathcal{N}(0,e)\end{align}

動的市場反応モデルは, 参照価格 \mathit{RP}_{t} という観測されていない変数があるが, a さえ決まれば p_{t} から一意に求めることができる. よって, まず, a を決めて \mathit{RP}_{t} を計算することで, 観測変数であるかのように扱う. 具体的には, 以下のような手順となる.

  1. モデルを dlm/KFAS/bsts の関数で記述する.
  2. a をグリッドサーチしつつモデルの尤度が最大になるように以下の操作を行う.
    1. 未知のパラメータである, 自己回帰項の構造パラメータ \phi_{1}, \phi_{2} と各ノイズの分散パラメータを最尤法で求める.
  3. \mathit{RP}_{t} を計算する.
  4. 全てのモデルについて, AIC を計算して最も値の良いものを採用する
  5. 状態空間モデルにカルマンフィルタを適用し, 各期間の状態変数を計算する.

状態空間モデルの予測計算を行う R パッケージには, dlm, KFAS, bsts などがある. このうち KFAS は日本語の情報が少ないが, 野村 (2016, カルマンフィルタ–Rを使った時系列予測と状態空間モデル–), 伊東 (2017, R による状態空間モデリング) などに言及がある. bsts パッケージは最近登場したため, 英語・日本語ともにいっそう情報が少ない. 日本語ではわずかに某有名ブログで紹介されている程度である.

英語ならば, 「Google 非公式データサイエンスブログ」では, 開発者の Scott 氏本人による わかりやすいチュートリアルが書かれている*3.

http://www.unofficialgoogledatascience.com/2017/07/fitting-bayesian-structural-time-series.html

また, Scott 氏のサイトにも解説スライドがある.

https://sites.google.com/site/stevethebayesian/googlepageforstevenlscott/course-and-seminar-materials/bsts-bayesian-structural-time-series

2017/9/8: bsts パッケージのより詳しい説明を自分で書いた.

ill-identified.hatenablog.com

残りのdlmKFAS だが, 上記のサポートページではこのうち dlm のサンプルプログラムが掲載されているが, 自己回帰項が省略されている, \mathit{RP}_{t} の値が既に用意されている, などの違いがある. そこで, dlm, KFAS, bsts の3つのプログラムを書いて比較した*4. 詳しい使い方まで書くのは大変なので, それぞれのパッケージの使い方を知りたい場合は, dlm ならば Petris (2010, An R Package for Dynamic Linear Models), KFAS ならば Helske (2017, KFAS : Exponential Family State Space Models in R)野村 (2016, カルマンフィルタ–Rを使った時系列予測と状態空間モデル–) をそれぞれ読んで欲しい*5.

dlm の場合

dlm は, ローカルトレンド成分, 季節調整成分, 回帰成分, といったように状態空間モデルの部品ごとに関数が与えられており, これを ggplot2 パッケージのように足していくことで状態空間モデルを記述できる. また,dlm() 関数を使えば構造パラメータ行列で記述できる.

パラメータを推定する場合は, 状態空間モデルを記述した dlmMod オブジェクトを返す関数を定義して dlmMLE() 関数に与えて最適化計算を行う.

ところで, 今回は自己回帰項を加える必要があるが, 定常自己回帰とするには係数に制約条件を与えなければならない. dlmMLE() 関数は optim() を用いてパラメータを推定するので, 厳密な意味で定常条件を満たしているとは言えないが, 自己回帰項の係数がいずれも \pm1 になるように制約を与えている.

また, dlm の弱点として, 時変パラメータのあるモデルは長期予測ができないということが挙げられる.

r - Gaussian state space forecasting with regression effects - Cross Validated では, 1期先予測を繰り返すことで長期予測をするというトリックが紹介されているが, 信頼区間を計算することはできない.

以下がそのプログラムである. ただし, パッケージやデータの読み込みなど, 共通の部分は common.R というプログラムにまとめてある. なお, 最尤法計算の際にいくつかエラーが出る箇所があるが, たぶん初期値が悪いからだと思われる (時間がかかるのでチェックしていないが, アルゴリズムか初期値を変更すればうまくいくと思う).

common.R

gist1fb70932d63725b29490edfd54726ff3

dlm.R

gist7640134a0330411e8db872ce084bf451

KFAS の場合

一方, KFASSSModel() 関数の中で formula オブジェクトのように, y~x1+x2 というふうに状態空間モデルを記述することができる. たとえばローカルトレンドモデルならば y ~ SStrend(1), 時変パラメータ回帰ならば y ~ SSMregression(~x1+x2) などというふうにできる. 状態空間モデルを記述した SSModel オブジェクトの要素は, dlm と同様に, 例えば SSModel$T などにアクセスすると, 構造パラメータを見ることができ, さらには, SSMCustom() 関数で構造パラメータ行列を直接記述することもできるから, 状態空間表現に慣れている人間には使いやすいだろう. 構造パラメータの名称と状態空間モデルは以下のように対応している.

\begin{align}
\boldsymbol{y}=\mathbf{Z}\boldsymbol{\alpha}_{t}+\boldsymbol{\varepsilon}_{t}, & \boldsymbol{\varepsilon}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{H})\\
\boldsymbol{\alpha}_{t+1}=\mathbf{T}\boldsymbol{\alpha}_{t}+\mathbf{R}\boldsymbol{\eta}_{t}, & \boldsymbol{\eta}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{Q})\end{align}

ただし, ARの係数など, 分散パラメータ以外に未知のパラメータがある場合の推定方法が少し面倒である. こちらも最尤推定の計算には optim() が利用されているが, このような場合にパラメータを更新するために自作する関数が, dlmより面倒になっている.

パラメータを推定したら, KFS() 関数によりフィルタリングと平滑化を両方行える.

kfas.R

gistc4036e56131437a35732a95acce6faf7

bsts の場合

bsts パッケージの名称は, ベイズ構造時系列, Bayesian structural time series の略に由来する.bsts も, 状態空間モデルの記述が容易である. 状態空間モデルは list オブジェクトで記述され, 組み込みの rbind() 関数のような使い方をする, 様々な Add*() 関数を用いて要素を追加していく. たとえばローカルレベル要素なら AddLocalLevel(), 季節性トレンドならば AddSeasonal() がある. また, 非ガウシアンな状態空間モデルにもある程度対応している. bsts が他2つと大きく違うのは, ベイジアンの名を冠するゆえんでもあるが, パラメータの推定を, 最尤法ではなく, 事前分布を設定してモンテカルロ法で計算する点である. 結果を視覚化するのも, plot() 関数だけでできて大変便利である.

ただし, 使用中に気づいた不具合? が1点ある. 観測データに欠損値が存在すると, ylim が設定されずにエラーになる. 現時点 (0.7.1) では欠損値のあるデータを扱う場合, ylim を手動で指定する必要がある.

以下がプログラムである. なお, bsts パッケージの理念を知ればわかるのだが, 今回の例にこのパッケージを適用するのはあまり適切ではない.

bsts.R

gistb1637eebee64b17ec255bfcacc8062f1


結論

しかし, 3通りのどのプログラムでも, 自己回帰項のないサンプルプログラムよりかなり当てはまりが悪いという結果になった. 明確に AIC が悪化しており, 平滑化パラメータも以下のようになっている.


AR項なし
f:id:ill-identified:20170906000146p:plain


AR項あり
f:id:ill-identified:20170906000218p:plain

「データが本書掲載と異なる」と注意書きがあったが, 自己回帰成分を除去したデータなのだろうか……? あるいは定常ARの推定がうまくできてないような気もする. きれいな終わり方ではないが, 3種類のパッケージの使い方の具体例を掲載するのにけっこう時間をかけてしまったのでこのまま掲載することにする*6.

では, 3種類のパッケージの違いについて要約する.

  1. dlmは全般的に構文が簡単だが, 非ガウシアンモデルには対応しておらず, 時変パラメータのあるモデルでは長期予測の信頼区間を計算できない.

  2. KFAS は多くの構文が簡単だが, パラメータ推定のプログラムが少し煩雑になる. 非ガウシアンモデルにも対応している. そのため, 時変パラメータのあるモデルでも長期予測の信頼区間を求められる.

  3. bsts は最も操作がシンプルで, 特に, 計算結果のプロットが他2つと比較して非常に簡単である. モンテカルロ法を使用するため非ガウシアンモデルにも対応できるが, パラメータ推定の計算時間が不安定で, 事前分布の設定など知識が必要になってくる.

そして, どのパッケージも, ARMA の定常条件の制約を厳密には設定できないことに注意すべきである.

最後に, 動的市場反応モデルについて. 今回の事例では, 「小売店における山積み陳列」という施策が売上にどう影響するかを組み込んだモデルだった. 同様の小売業, あるいはEコマースなどの時系列データであれば, 山積みに限らず様々な施策のもたらす影響について分析できることだろう.

ただし, 自分はこのモデルに対し, 1つだけ気になる点がある. 今回の動的市場反応モデルは, 毎期の参照価格は過去の参照価格価格によって決まるという設定になっている. しかし, スーパーマーケットが定期的に値引きするということがよくある. 例えば毎月29日は「肉の日」で肉製品を値引きする, というふうに. もし消費者が, 過去だけでなく将来の価格についても予測して参照価格を形成できるとしたら, 過去の価格だけを参照するという今回のモデルは不適切になりうる. 思うに, このような数値の外にある情報をもうまくモデルに取り込むのが, 統計モデリングの面白さではないだろうか.

参考文献


*1:原典と記号の使い方が異なるので注意.

*2:dynamic regression とも呼ぶ.

*3:英語が読めるのならこちらを読むだけで機能のほとんどを理解できるだろう.

*4:ただし, Z の計算部分は本質的に変わらないので, dlm のプログラムでのみ書いた.

*5:もしかしたらそのうち使用法解説記事を書くかもしれない

*6:似たようなデータを扱うツテがあれば, もっと面白いことが出来たかもしれないが……