« 統計学者のみなさん、記法は統一してもらえないかしら (状態空間モデルの巻) | メイン | 読了: Petris & Petrone (2011), Petris (2010) dlmパッケージとそのライバルたち »
2014年10月17日 (金)
以下は純粋に個人的なメモであります。
ここにn変量時系列 $y_t$ がございます。
これを以下のように分解する。まずは
$y_t = B + Lf_t + u_t, \ \ \ u_t \sim N(0, \Sigma_u)$
$B$は切片(サイズ n x 1), $L$は係数行列(n x p), $f_t$はp変量の潜在動的因子(p x 1), $u_t$はn変量の不規則要素(n x 1)である。
潜在動的因子をトレンドと季節に分解する。季節要素を観測値のモデルにではなく因子の時系列モデルのほうに入れている点に注意。
$f_t = \alpha_t + \gamma_t$
トレンドのほうは
$\alpha_t = \alpha_{t-1} + \beta_{t-1} + \epsilon_t, \ \ \ \epsilon_t \sim N(0, \Sigma_\epsilon)$
トレンドの傾きも時間変動させて
$\beta_t = \beta_{t-1} + \delta_{t-1} + \eta_t, \ \ \ \eta_t \sim N(0, \Sigma_\eta)$
傾きの変化率をランダムウォークさせる。
$\delta_t = \delta_{t-1} + \zeta_t, \ \ \ \zeta_t \sim N(0, \Sigma_\zeta)$
季節のほうも確率的に変動させる。以下、添え字が煩雑になって混乱するので、季節周期を4に決め打ちする。
$\gamma_t = -\sum_{j=1}^{3} \gamma_{t-j} + \xi_t, \ \ \ \xi_t \sim N(0, \Sigma_\xi)$
以上、すごく複雑に見えるけど、因子の時系列構造を工夫しているだけで、要するに因子分析である。びびってはいけない。
さて、状態空間表現ではどうなるんだよ、という話である。さあ深呼吸！
観察方程式は
$y_t = A z_t + B + u_t, \ \ \ u_t \sim N(0, \Sigma_u)$
ここで状態変数ベクトルは
$z_t \equiv [\alpha_t \ \beta_t \ \delta_t \ \gamma_t \ \gamma_{t-1} \ \gamma_{t-2}]'$
実際には$\alpha_t, \beta_t, \ldots$は縦ベクトル(p x 1)なので、こういう書き方でいいのかどうかよくわかんないんだけど(自慢じゃないけど高校で数学の時間はずっと寝ていた)、原文でもこうなってるし、まあいいことにしよう。なお、原文で最後の項は$\gamma_{t-s-2}$となっているが、$\gamma_{t-(s-2)}$の誤りであろう。
係数行列は
$A \equiv [L_{n \times p} \ 0_{n \times p} \ 0_{n \times p} \ L_{n \times p} \ 0_{n \times p} \ 0_{n \times p}]$
順に、$\alpha_t$の係数, $\beta_t$の係数、$\delta_t$の係数, 季節ダミー1の係数($\alpha_t$の係数と同じ。なぜなら季節要素が状態方程式に突っ込んであるから), 季節ダミー2の係数、季節ダミー3の係数、である。結局、$z_t$の長さは$m = 6p$。季節周期を$s$として、$m=(2+s)p$である。
さあ状態方程式だ！
$z_t = C z_{t-1} + v_t, \ \ \ v_t \sim N(0, \Sigma_v)$
状態撹乱要素は
$v_t \equiv [\epsilon_t \ \eta_t \ \zeta_t \ \xi_t \ 0 \ 0]'$
問題は遷移行列$C$。原文では次のようになっている。
$C \equiv \left( \begin{array}{llllll} I_p & I_p & 0 & 0 & 0 & 0 \\ 0 & I_p & I_p & 0 & 0 & 0 \\ 0 & 0 & I_p & 0 & 0 & 0 \\ 0 & 0 & 0 & -I_p & \cdots & -I_p \\ 0 & 0 & 0 & I_p & 0 & 0 \\ 0 & 0 & 0 & 0 & I_{(s-3)p} & 0 \end{array} \right)$
本当かなあ？ ここにどうも納得がいかなかったのである。
いちから考えてみたい。$C$の右側にあるのは、たとえば$[\alpha_4 \ \beta_4 \ \delta_4 \ \gamma_4 \ \gamma_3 \ \gamma_2]'$である。これに$C$を掛け、つくりたいのは$[\alpha_5 \ \beta_5 \ \delta_5 \ \gamma_5 \ \gamma_4 \ \gamma_3]'$である。以下、大きな行列を書くのが面倒なので、$C$をp行p列づつ区切り、上の段から順に考える。
最初のp行は$\alpha_5$をつくる奴だから、
$(I \ I_ \ 0 \ 0 \ 0 \ 0)$
次のp行は$\beta_5$をつくる奴だから
$(0 \ I \ I \ 0 \ 0 \ 0)$
次のp行は$\delta_5$をつくる奴だから
$(0 \ 0 \ I \ 0 \ 0 \ 0)$
次のp行は$\gamma_5$をつくる奴だから
$(0 \ 0 \ 0 \ -I \ -I \ -I)$
次のp行は$\gamma_4$をつくる奴だから、
$(0 \ 0 \ 0 \ I \ 0 \ 0)$
最後のp行は$\gamma_3$をつくる奴だから、
$(0 \ 0 \ 0 \ 0 \ I \ 0)$
ほんとだ...！ 合っている！
やれやれ。分散行列 $\Sigma_v$はブロック単位行列で、左上から順に$\Sigma_\epsilon, \Sigma_\eta, \Sigma_\zeta, \Sigma_\xi, 0_p, 0_p$である。
なお、実際にはすべての分散行列を対角行列に制約した由。$\Sigma_v$のみならず、観察方程式の分散行列$\Sigma_u$も対角行列にしちゃうのである。不規則要素に共分散はないと考えたわけだ。ま、いいか。
パラメータ推定はEMアルゴリズムで行った由。詳細は略するが、ステップのなかでバリマクス回転している。いやーん。そんなんようせんわ。
なお、データはあらかじめ対数変換した由(歪度が大きかったから)。さらに平均0, 分散1に標準化した由(分析の主旨は動的因子にあるので、別にかまわない)。
以上、Du & Kamakura (2012, JMR)のWeb Appendix Bより。とても歯が立たないだろうと敬遠していたが、実際に読んでみたら、もちろん難しいんだけど、思ったほどではなかった。
本文を読んだときも思ったことだが、季節要素をこうして動的因子系列に入れるべきかどうかはデータの実質的な性質によるなあ、と思った。著者らは自動車ブランド名の検索量の時系列を分析していて、因子として「欧州車」というようなものを想定しているので、それなりに筋が通っている(欧州車に共通の季節要素があると考えているわけだ)。でもこれをマルチカテゴリに拡張して、著者らが主張するようにtrendspottingという観点から捉えたときには、観察方程式に入れたほうが自然ではないかと思う。
推定のところにオリジナルな工夫があるところがちょっと嫌だ。このモデルって、あらかじめ予備的分析で因子数と負荷行列にあたりをつけ(変量ごとに季節を除去しEFAにかけるとか)、$L$を何箇所か適当に固定すれば、最尤推定できたりしないのかなあ。
雑記：データ解析 - Du & Kamakura (2012) メモ