2018-05-14

雑記：カルマンフィルタとか何とかカルマンフィルタとか

雑記

参考文献：データ同化入門 (予測と発見の科学) | 樋口知之 |本 | 通販 | Amazon
文字の置き方と表式をだいたい上の本に準拠していますが、違うこともあります。導出の仕方は上の本とは異なります。間違っていたらご指摘ください。今回はカルマンフィルタの話しかないです。キャラクターに元ネタはないです。

次回：まだ

f:id:cookie-box:20180513082851p:plain:w60

カルマンフィルタとその仲間たちみたいな何とかフィルタって色々ありますよね。

f:id:cookie-box:20180513082908p:plain:w60

逐次ベイズ推定のアルゴリズムのこと？確かに色々あるね。拡張カルマンフィルタとか、粒子フィルタとか。

その辺ってなんかわかった気になってしまいませんか？線形でベイジアンだったらちゃんと解けばよくて、そうじゃなかったら現在の $x_t$ の周りで1次形近似したり、 $x_t$ の分布を粒子のアンサンブルで近似すればいいんだ、って感じで…。

なんで！？いうほどわかった気にならないよ？

というか、そもそも何をしたかったのかもよくわからなくなってしまって…。

目的は状態の推定だよね！？

よく ARIMA のような伝統的な時系列モデルから状態空間モデルへ、って説明されますけど、話が変わってる気がするんですよね。ARIMA では AR や MA の各次数の係数を求めて、状態空間モデルは状態の分布を推定していくんですよね。なんか話が逆になってませんか？モデル $f(x_t)$ に対して時系列データ $x_t$ を所与として $f$ の係数を求めたいのか、 $f$ の形は所与として $x_t$ を推定したいのか…。

うーん…状態って言葉が曖昧なのかな？観測値の真の値であるかのように説明されることも多いけど、というよりは状態って「トラッキングしたいもののベクトル」なんだよね。所与なのはあくまで観測値 $y_t$ で、ARIMA ではむしろ AR や MA の各次数の係数たちが $x_t$ で（ふつう時不変だから $x$ だけど）。状態空間モデルではモデルの形 $f$ は所与とするけどもしそのパラメータが未知ならそれは $x_t$ 側に含めると思うよ。なんていうか結局「確信しているものを所与として不確かなものを推定したい」だね。トートロジーだけど。

なるほど…。じゃあ、このブログの2018年4月23日の記事にある、「（状態空間モデルで）もし『状態』が観測される変数そのもので、『観測モデル』が恒等写像だったら（伝統的な時系列モデルと）一緒の構図になる」というのは対応がおかしかったんですね。より正しくは、「もし『状態』が ARIMA モデルの係数（及び必要なステップだけの過去から現在までのノイズや観測値）で、『観測モデル』がそれらの線形和（ARIMA モデルの形）だったら一緒の構図になる」ですね。観測モデルは観測といっても日常的な意味のように「観測時に入り混じる誤差を考慮したモデル」ではなく、「不確かな変数を知覚できる数値に変換するモデル」なのに、筆者はそこを踏み外して、思考停止でシステムモデルに ARIMA を代入していたんですね。

う、うん。たぶんそうだね…。あくまでそういう状態空間モデルにすれば ARIMA と「係数を求めたい」という目的が大雑把にそろうというだけで、実際には ARIMA モデルの方では逐次的に係数の分布を更新してはいかないと思うけどね。あと、もしそのように修正する前の元の文章の立場に立つなら、「この時系列データはこの係数のARIMAモデルで時間発展する」ということの方に絶対の確信があって、知覚された時系列データがそれに従っていないならそっちがあやしい、観測ノイズによって真の値とずれているんだろう、真の値が知りたい、という姿勢にはなるね。

不確かなものを推定したいという目的はわかったと思います。では早速各種フィルタについてちゃんと追って行きたいと思うんです。まず一般的に、一期先予測とフィルタリングは以下の式で表せます。
　　 $\displaystyle p(x_t | y_{1:t-1}) = \int p(x_t | x_{t-1}) p(x_{t-1} | y_{1:t-1}) dx_{t-1}$ 　（一期先予測分布）
　　 $\displaystyle p(x_t | y_{1:t}) = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{p(y_t | y_{1:t-1})} = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{\displaystyle \int p(y_t | x_t) p(x_t | y_{1:t-1}) dx_t}$ 　（フィルタ分布）

また端折ったね…。まず前提として、ここでは時刻 $t-1$ までの観測値 $y_{1:t-1}$ が得られたもとでの状態の分布 $p(x_{t-1} | y_{1:t-1})$ を所与とするんだね。かつ、時刻 $t-1$ の状態を $x_{t-1}$ と固定したもとでの次の時刻の状態の分布 $p(x_t | x_{t-1})$ も計算できるとする。これはシステムモデルを決めておきなさいってことだね。これらがわかれば、一期先予測の式は $x_{t-1}$ の分布にわたって次の時刻の $x_t$ がどうなるかの分布を積分するだけ。それで次にフィルタの式は、現時点までの観測値 $y_{1:t}$ が出そろったもとでの状態の分布 $p(x_t | y_{1:t})$ はどうなっているか、ということだけど…これは観測モデル $p(y_t | x_t)$ とさっきの一期先予測の式 $p(x_t | y_{1:t-1})$ の積からベイズの定理により求まるね。ただ、単にこれらの積だと確率分布が正規化されないから、 $p(y_t | y_{1:t-1})$ で割る必要がある。これは、 $p(y_t | y_{1:t-1})$ は $p(y_t | x_t) p(x_t | y_{1:t-1})$ を $x_t$ の分布にわたって積分すれば求まる。あと、ここまで暗に、下図のグラフィカルモデルで表される従属構造が仮定されていることに注意が必要だね。つまり、 $p(x_t | x_{1:t-1}, y_{1:t-1}) = p(x_t | x_{t-1})$ 及び $p(y_t | x_{1:t}, y_{1:t-1}) = p(y_t | x_t)$ が成り立つ。マルコフ性ともいうね。そうじゃなきゃ上の式にならない。

※ ここに状態空間モデルのグラフィカルモデルを貼る。

あ、はい、色々飛ばしてました…補足ありがとうございます。それでまずカルマンフィルタの導出です。ここからは、システムモデルと観測モデルがそれぞれ以下のようにかけると仮定しています。
　　 $\begin{cases} x_t = F_t x_{t-1} + G_t v_t \; \; \; \; & v_t \sim N_m(0, Q_t) \\ y_t = H_t x_t + w_t & w_t \sim N_l(0, R_t)\end{cases}$
上式に出てくる以下の文字は以下のような次元のベクトルや行列です。 $x_t \in \mathbb{R}^k, \; y_t \in \mathbb{R}^l, \; F_t \in \mathbb{R}^{k \times k}, \; G_t \in \mathbb{R}^{k \times m}, \; H_t \in \mathbb{R}^{l \times k}, \; Q_t \in \mathbb{R}^{m \times m}, \; R_t \in \mathbb{R}^{l \times l}$
このうち $Q_t$ と $R_t$ は正規ノイズの分散共分散行列ですね。それで、 $x_t, \; v_t, \; y_t, \; w_t$ が最初から決まっていない、確率変数です。目的は、 $x_t$ の分布をトラッキングすることです。

うん、観測値より状態の次元の方が大きいなら $H_t$ は横長の、逆なら縦長の行列になるね。

そうですね。それで、先に結論を書きます。
　　 $x_{t-1} \, | \, y_{1:t-1} \sim N_k (\mu_{t-1|t-1} , \, V_{t-1|t-1})$
を仮定すると、
　　 $x_{t} \, | \, y_{1:t-1} \sim N_k (\mu_{t|t-1} , \, V_{t|t-1}) = N_k (F_t \mu_{t-1|t-1} , \, F_t V_{t-1|t-1} {F_t}^{\top} + G_t Q_t {G_t}^{\top})$ 　（一期先予測分布）
　　 $x_{t} \, | \, y_{1:t} \sim N_k (\mu_{t|t} , \, V_{t|t}) = N_k (\mu_{t|t-1} + K_t e_t , \, V_{t|t-1} - K_t H_t V_{t|t-1})$ 　（フィルタ分布）
但し、
　　 $e_t = y_t - H_t \mu_{t|t-1} \; \in \mathbb{R}^l$
　　 $K_t = V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} \; \in \mathbb{R}^{k \times l}$
が成り立ちます。時刻 $t-1$ までの観測値 $y_{1:t-1}$ が得られたもとで推定されている状態の分布がガウシアンならば、一期先予測もガウシアンであり、さらに時刻 $t$ のフィルタ分布もガウシアンであるということです。かつ、その一期先予測とフィルタ分布の平均ベクトルと分散共分散行列は上式のように1つ手前の平均ベクトルや分散共分散行列の式として陽に書き下せます。

有名な結論だね。

でもどうしてこういうことになるのか。

まず一期先予測には以下をつかおうか。つまり、多変量正規分布にしたがう確率ベクトルは線形変換しても多変量正規分布にしたがうし、多変量正規分布にしたがう確率ベクトルどうしを足しても多変量正規分布にしたがう。証明は、前者は $M_{Ax}(t) \equiv E(\exp(t^{\top} Ax))=M_x(A^{\top} t)$ を、後者は、互いに独立な確率ベクトル $x_1, \, x_2$ について $M_{x_1 + x_2}(t) = M_{x_1}(t) M_{x_2}(t)$ を利用すれば導出できるかな。あ、 $M_x(t)$ はモーメント母関数で、これは確率分布に対して一意に定まる関数で、多変量正規分布のモーメント母関数は $M_x(t) = \exp(t ^{\top} \mu + \frac{1}{2} t ^{\top} \Sigma t)$ だね。

$x \sim N_k (\mu_x, \, \Sigma_x)$ で $y = Ax$ のとき、 $y \sim N_l (A \mu_x, \, A \Sigma_x A^{\top})$

$x_1$ と $x_2$ が独立で、 $x_1 \sim N_k (\mu_{x_1}, \, \Sigma_{x_1}), \; x_2 \sim N_k (\mu_{x_2}, \, \Sigma_{x_2})$ のとき、 $x_1 + x_2 \sim N_k (\mu_{x_1} + \mu_{x_2}, \, \Sigma_{x_1} + \Sigma_{x_1})$

それらをつかう……？えっと、いま $x_{t-1} \, | \, y_{1:t-1}$ が多変量正規分布にしたがうと仮定しているので、その線形変換である $F_t x_{t-1} \, | \, y_{1:t-1}$ も多変量正規分布にしたがう。それとは独立にノイズ $G_t v_t$ も多変量正規分布にしたがうから… $x_{t} \, | \, y_{1:t-1} = F_t x_{t-1} \, | \, y_{1:t-1} + G_t v_t$ は多変量正規分布にしたがいますね。それで、その平均ベクトルと分散共分散行列は、
　　 $\mu_{t|t-1} \equiv {\rm E}(x_{t} \, | \, y_{1:t-1}) = {\rm E}(F_t x_{t-1} \, | \, y_{1:t-1}) + E(G_t v_t) = F_t \mu_{t-1|t-1}$
　　 $V_{t|t-1} \equiv {\rm Var}(x_{t} \, | \, y_{1:t-1}) = {\rm Var}(F_t x_{t-1} \, | \, y_{1:t-1}) + {\rm Var}(G_t v_t) = F_t V_{t-1|t-1} {F_t}^{\top} + G_t Q_t {G_t}^{\top}$

一期先予測分布が導出できたね。フィルタ分布はベイズの定理を使うとどうなる？さしあたり正規化因子は考えなくていいよ。

ベイズの定理をつかうと、
$\begin{split} p(x_{t} \, | \, y_{1:t}) &\propto p(y_t | x_t) p(x_t | y_{1:t-1}) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) \right) \exp \left( -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \end{split}$
ここからどうすれば…。

$\begin{split} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) + (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \end{split}$
この式を $x_t$ があらわれない項は無視して展開して整理すると、

$\begin{split} & - {y_t}^{\top} {R_t}^{-1} H_t x_t - {(H_t x_t)}^{\top} {R_t}^{-1} y_t + {(H_t x_t)}^{\top} {R_t}^{-1} H_t x_t + {x_t}^{\top} {V_{t|t-1}}^{-1} x_t - {x_t}^{\top} {V_{t|t-1}}^{-1} \mu_{t|t-1} - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1} x_t \\ &= {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) x_t - {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} y_t + {V_{t|t-1}}^{-1} \mu_{t|t-1}) - ({y_t}^{\top} {R_t}^{-1} H_t - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1}) x_t \end{split}$

こうなるけど、 $x_t$ に挟まれている ${H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}$ は、逆行列補題より逆行列 $V_{t|t-1} - ( V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1}) H_t V_{t|t-1}$ をもつ。加えて、 ${x_t}^{\top}$ に右側からかかっている行列と $x_t$ に左側からかかっている行列が対称だから、上の式は結局 $\begin{split} & ({x_t} - \alpha)^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) (x_t -\alpha) \end{split}$ の形に整理することができる。これを展開して係数比較すれば $\alpha$ が求まるね。

※ ここに逆行列補題をかく。

そして、 $\begin{split} & ({x_t} - \alpha)^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) (x_t -\alpha) \end{split}$ の形に整理できるということは、フィルタ分布 $p(x_{t} \, | \, y_{1:t})$ が多変量正規分布にしたがうということに他ならないよね。平均 $\alpha$ で、分散共分散行列 $({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1})^{-1}$ のね。

…ごめん寝てた。

なんで寝るの！？部長のために式展開したのに！

次のフィルタの話に移りたいです。ざっとモンテカルロ手法まで復習したく…というのも、そろそろ部活動の予算申請の時期なので、我らがベイズ統計部も膨大な計算資源の必要性を申請書にしたためなければ。

部活動予算でどんなマシン買う気なの！？運動部じゃないんだからそんな予算下りないよ？

（次回があれば）つづく

2018-04-24

論文読みメモ： Unsupervised Anomaly Detection with GAN（その1）

論文読み

以下の論文を読みます。

Thomas Schlegl, Philipp Seeböck, Sebastian M. Waldstein, Ursula Schmidt-Erfurth, Georg Langs. Unsupervised Anomaly Detection with Generative Adversarial Networks to Guide Marker Discovery. arXiv: 1703.05921, 2017. https://arxiv.org/abs/1703.05921

※ 以下、キャラクターが会話します。原作とは関係ありません。

次回：まだ

f:id:cookie-box:20180405221013p:plain:w60

この論文は、何が異常かって正解を与えることなしに、画像から異常を検知するって話みたいです！「この画像はこういう異常のデータ」って正解ラベルを付けて学習するのって、大量にラベルを付けるのが大変だし、学習データになかった異常に対応できないから、正解なしで学習したいって。GANってモデルと、画像から隠れ変数空間へのマッピング手法？を組み合わせたらしいです。あ、この論文で想定している画像というのは、OCT（光干渉断層撮影）スキャンといって、病院での検査などに用いる、皮膚や眼球の断面が撮影できる技術みたいです。網膜のOCTスキャンデータに対する検証で、病気状態の画像を異常として正しく特定できたそうです。…でも絵理さん！このモデルは、画像を入力すると、異常な画像かどうかを出力してくれるんですよね？何が異常かっていう正解を与えてないのに、どうやってこの画像は異常だろうってわかるんですか？おかしくないですか！？

f:id:cookie-box:20180405222314p:plain:w60

確かに、画像から隠れ空間へのマッピングによって異常を判定するってことは、隠れ空間内のどの座標にマッピングされたら異常っていうのがあるんだよね。じゃあGANっていうモデルは、正常な画像は正常な画像に対応する領域に、異常な画像は異常な画像に対応する領域にうつるように画像から隠れ変数空間への写像を学習するはずだけど、画像に正常か異常かのラベルが貼ってないのにそんな学習はできないんじゃないかな。

f:id:cookie-box:20180405220800p:plain:w60

…この論文では、網膜のOCTスキャン画像のうち、網膜下液という症状のない箇所の画像のみを学習データにしている。だから、モデルは正常な画像と異常な画像をごちゃ混ぜに正解なしで与えられてるわけじゃない。正常な画像のみを学習データとして、正常な画像の特徴を学んでおいて、未知の画像の異常度を判定する。

うーん…それって、教師なし学習っていうんですか？正常な画像か異常な画像かを判定するのに、正常な画像のみで学んでおきましたって…。犬の画像か猫の画像しか含まれていないデータを分類するのに、犬の画像をバッチリ学習しておきました、みたいな話ですよね。それってある意味正解ラベル付きみたいなものだし、予め学習したデータにあてはまりそうかそうでないかを判断するだけじゃないですか。

そう、予め知っているデータ側なのか、そうでない側なのかという話。でもそれを判定する方法は、自明じゃない。例えば、何かデータ間の距離を定義して、予め知っているデータに近いか判定する？未知の画像が犬か猫かを判定したいとして、予め学習した犬の画像のどれにも近くなくても、犬の画像かもしれない。予め学習したある犬の画像に近くても、猫の画像かもしれない。教師なし学習で代表的なクラスタリングは、データの中にクラスタの分布たちを見出す。異常検知もまた、与えられた正常データを表現する正常な分布を見出す。異質なものとの比較によって自分の特徴を把握できないから、後者はきっと前者よりもっと難しい。教師なし学習か教師あり学習かという区別は、いま本質的じゃない。

じゃあ、GAN ってモデルはどうやって正常データの分布を学習するんだろう。

生成モデルと識別モデルを同時に学習するって書いてありますね。ある1つのコスト関数を最小化するのではなく、生成モデルの表現力（識別器に本物と判定される確率）と識別モデルの真贋の識別力のナッシュ均衡を目指すみたいです。

（その2があれば）つづく

2018-04-23

論文読みメモ： Forecasting Stock Price with EKF

論文読み

以下の論文を読みます。

H. Haleh, B. Akbari Moghaddam, and S. Ebrahimijam. A New Approach to Forecasting Stock Price with EKF Data Fusion. International Journal of Trade, Economics and Finance, 2(2):109-114, 2011. http://www.ijtef.org/papers/87-F00046.pdf

※ 以下、キャラクターが会話します。それぞれの原作とは関係ありません。内容の誤りは本ブログ筆者に帰属します。

この論文の内容は、拡張カルマンフィルタで株価を予測する？

株価が予測できるんですか？夢がある話ですね！

何そのテンプレみたいな反応。

株式の真の価値と市場価格に乖離があるために、従来の分析では十分なパフォーマンスが出ないって問題提起してる。

市場価格とは別に株式の真の価値ってのがあるんですか？市場で出回ってる価格が株式の価値そのものなんじゃないんですか？

株式の本来の価値を測る方法っていうと、企業が将来獲得するキャッシュフローを現在価値に割り引く方法とか、配当方針から逆算する方法とかあるけど、そういうのが絶対的に株式の真の価値ってわけでもないし、「株式の真の価値と市場価格の乖離があるために上手くいかない」っていうくだりは「何か真の価値というものがあると考えた方が上手くいく」と読み換えた方がいいんじゃないかな。というかそれでちゃんと上手くいったんだよね？

テクニカル分析とファンダメンタル分析の情報を拡張カルマンフィルタで融合することで、イランのある工業関連の会社の株価について、回帰やニューラルネットよりもよい短期予測ができたとアブストに書いてある。

テクニカル分析とファンダメンタル分析？何かつかみどころがない言葉です。

この時点では株価を過去時系列に基づいて予測する手法をテクニカル分析といい、株式会社の持分であるところの株式の価値を業績予測を通して見積もる手法をファンダメンタル分析というイメージでいいと思う。大丈夫、何もつかみどころがないのは株価だけじゃない。トレーダーでなくたって、あらゆる人は明日どうなるか読めない人生を生きている。だから心配ない。

心配だよ！？

イントロは、株式投資の重要性？こうくるのはちょっと面白いかも。のっけから株価予測は投資家の関心事である、みたいなのが多いから。株価に影響するたくさんの情報を上手く扱う方法がない、だから上手く扱えるようにしたい、と。そして効率的市場仮説の紹介とそれへの反論。ここでいわれている市場の非効率性は元論文までみていないので詳しくわからない。そして既存の手法の欠点を非線形性への未対応や訓練データ不足に帰している。

訓練データが足りないって理由になるんですか？訓練データ足せばいいじゃないですか。

ここでの訓練データ不足は、訓練データが足りないために未知の状況に遭遇するといった意味みたいだから、確かに違和感があるかも。あと予測が困難な理由として、市場価格が理論モデルにしたがわないとも書いてる。理論モデルとは何だったのか。それでこの論文では拡張カルマンフィルタを用いて、非線形性への対応と、テクニカル指標とファンダメンタル指標の融合を実現したって。

そうそう絵理ちゃん、その拡張カルマンフィルタって？

この論文でも説明されているけど…ちゃんとした本を参照した方がよさそうだからとばす。

とばすの！？

111ページの Figure 2. が提案モデルの模式図だけど、本文を読み進めると、ファンダメンタル指標から配当割引モデルのゴードンモデルによって直接観測できない株価 $x_k$ を計算している。でも、テクニカル指標がどう入り込んでいるかがよくわからない。これだと、観測モデルをかませることそのものがテクニカル指標の融合だと言わんばかりのような。あ、違う。 $h(x_k, v_k)$ というのが (28) 式で、これはこれでこういう風に株価を推測する方法？なんで2次関数？

テクニカル指標に詳しくないのでよくわかりません。

ぐだぐだだけどおわり

2018-04-23

隼時系列本：ノート7

本読み

以下の本を読みます。

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装
馬場真哉

プレアデス出版 2018-02-14
売り上げランキング : 7742

Amazonで詳しく見る by G-Tools

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。
※ アイコンまだです。

前回：ノート6 ／次回：まだ

f:id:cookie-box:20180305231302p:plain:w60

後半からはナツキが本読み手伝ってくれるんだな。よろしく。

f:id:cookie-box:20180423175223p:plain:w60

うん。ジュンは学校も仕事も忙しいから、最後まで付き合うのは難しい。だから、後は俺に任せて…。

いや俺も学校も仕事も忙しいけどな！？ …でも、ナツキって前半読んでないんだよな？いきなり後半からで大丈夫なのか？

大丈夫、この本の前半と後半は比較的別々に読めるはず…ハヤト、前半の話をまとめられる？

この本の前半の話は…まず、時系列分析っていうのは、手元の時系列データの将来値の予測をしたりその時系列データについて何かしらの知識を得ることが目的。その目的のために何をするかっていうと、時系列データには各時点での値の背後に確率分布があると仮定して、尤もらしい確率分布列を与えるモデルを特定する。そのモデルの候補として、ARIMAやVAR、GARCHとかが紹介されてた。ただ、このモデルを同定するときの注意点として、予め対象時系列データを定常にしておくことや、あてはめたモデルで説明できない残差の正規性や自己相関の確認が要るんだったな。そのために用いる統計的検定も色々出てきた。

まとめありがとう。後半の「状態空間モデル」も、目的はだいたい一緒だよ（「状態」を明示的に扱うなりの目的も加わってくるかもしれないけど）。でも、目的のために何をするかがもっと一般的になる。時系列データには各時点での値の背後に、直接観測できない「状態」と、状態をその時点まで運んできた「システムモデル」、状態を目に見える観測値の分布に写す「観測モデル」があると考えるよ。もし「状態」が観測される変数そのもので、「観測モデル」が恒等写像だったら、前半までの話と一緒の構図になるけどね。

直接観測できない「状態」があると考えるって？

本の179ページにも魚の数の例があるけど、「状態」を導入する意義の1つは、時系列データって、「観測するときにだけ混ざってくるノイズ」があるかもしれない。例えば、俺たちが持ち回りで部室の気温を毎日記録するとして、ハヤトが記録する日だけはハヤトが温度計に近づくせいで温度計の指示値が少し高めに記録されるかもしれない。でも、そのハヤトのせいで記録値が少し大きくなった分は、翌日以降の気温に影響する成分とは思えないよね。

なんていうか、それはもう俺を記録係から外した方がいいな…。

この本の前半のモデルにも誤差項はあるけど、それはその場限りの誤差ではなくて、その影響はその後ずっと（MAモデルの場合はどこかで打ち切られるけど）続いていくものになっている。そうじゃなくて、現実に何かを測るときはその場限りの誤差ってあるはず。だったら、その場限りの誤差が付加される前の「状態」の推移を考えればいい。それが状態空間モデルの意義。それに加えて、第4部の冒頭には明示的に書いてないと思うんだけど、状態って何も「観測される変数にその場限りのノイズが付加される直前のもの」とは限らなくて、「観測される変数そのものではないけど、その推移を追うと有用な量」だったり「観測できないけど興味がある量」かもしれない。状態空間モデルならそういう設計もできる。例えば、ハルナが毎日食べるドーナツの個数は直接観測できて、ハルナの毎日の空腹度は直接観測できないけど、したがう方程式がわかっているのは後者という場合とか。

だから、ハルナの空腹度がしたがう方程式が研究で明らかにされてたらおかしいだろ！ …でもナツキ、状態空間モデルでは、その状態がしたがう方程式が予めわかってないといけないってこと？前半ではさ、ARIMA だったら ARIMA の次数と係数は AIC なり何なりを指標に計算機にゴリ押しで決めさせるって感じで、「対象時系列はこんなメカニズムをもつんじゃないか」って考えるプロセスがなかった気がするんだけど…。

ARIMA を使ってみる時点である意味「ARIMA で表現できるんじゃないか」っていう考えがあるとも思うけど…状態空間モデルでは、対象時系列に対して「こんなメカニズムをもつんじゃないか」というのがもっと明示的にある場合が多いんだろうね。もちろん、データに対する強い信念なしに、色々な次数の ARIMA を試すような用法もあると思うよ？システムモデルの係数を未知パラメータにして。182ページの最初で「このあたりの手順に関しては」って別の文献を紹介しているね。でもそういうモデル探索よりは、何かこういうシステムモデルだろうという信念が予めある程度あって、直接観測できない「状態」の推移やその分散を捉えたいって向きが強いんじゃないかな。それが180ページの「過去の知見（中略）自由にモデル化できる」「（前略）モデルの解釈が容易」ってことだと思うから。俺も普段から状態空間モデルを使うわけじゃないから、よくわからないけど…。

普段から状態空間モデルを使う高校生はいないだろ！

まあそれで、この後の中心話題は未知パラメータをどうやって推定するかだね。この本では2つの異なるアプローチのそれぞれの代表として、カルマンフィルタとHMC法を大きく取り上げてる。この辺について、ちょっとプロデューサーさんが昔書いた記事を掘り出しておくね。

カルマンフィルタに関する記事は以下など。
- 状態空間モデルのまとめ絵
- R dlm をつかってみるだけ
1時点前までの観測値を得た下での現在の状態の分布は、一時点前の状態を用いて、

だけど、これに現在の観測値が手に入った下での現在の状態はベイズの定理より

になる。もし状態空間モデルのシステムと観測が線形で、システムノイズと観測ノイズがガウシアンなら（線形・ガウス状態空間モデル）、をガウス分布とすると、もガウス分布、さらにもガウス分布、…となって、1つ手前を含む形の式になって、どんどん逐次的に計算していけるんだよね。これがカルマンフィルタだね。
HMC法に関する記事。後者の記事、収束してなさすぎだけどスクリプト合ってるのかなあ…。

もっと具体的な計算手順や、これら以外の状態空間モデルの解き方とかはまた今度にする？

うーん、183～185ページに今ナツキが言ったような2つのアプローチのパラメタ推定があるってあるけどさ、なんでその2つなのかよくわからないんだよな。要するに、真面目に解くか、乱数シミュレーションで解くかって違い？でも、逐次的かどうかも違って…真面目に同時に解くとか、乱数的に逐次的に解くとかはないの？

ないとは言い切れない、けど…真面目に全期間同時に解けることなんてまずないんじゃない？逐次的に解けるような特殊なケースだから真面目に書き下せたってだけで…。あと185ページに、HMC法ではパラメタの推定を同時にやるから逐次的にはしにくいとも書いてある。

そっか、どちらかというと、特殊なケースでは真面目に解けて、副産物的に逐次解法になっている、って感じなのかな…。あとナツキ、209ページの散漫カルマンフィルタってのは？

あまり他では聞かない用語、だけど…状態の分散の初期値をものすごく大きくしておいて、状態の期待値の初期値の影響をすぐ消滅させるテクニックを指しているみたいだね。

（ノート8があれば）つづく

2018-04-03

隼時系列本：ノート6

本読み

以下の本を読みます。

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。

前回：ノート5 ／次回：まだ

144ページからは、VAR モデル？ …これは、複数の時系列の過去の値で複数の時系列を回帰するモデルで、複数時系列を用いることによってよりよい予測ができるようにしたり、時系列化間の相互作用を調べる目的で利用するのか…具体的に、どんな時系列に対して適用するんだろう？

f:id:cookie-box:20180305232608p:plain:w60

例えば春名さんの毎月の収入と、毎月のドーナツ消費量という2つの時系列データがあったとして、春名さんは収入が増えるほどドーナツを消費するかもしれませんが、好物のドーナツを消費することによって仕事に精が出て収入も増加するかもしれません。と考えると、春名さんの収入とドーナツ消費量の間には相互作用があるのではないかと疑われます。まあこの例は本に載っていないのでGranger因果性が存在するかはわかりませんが。

本に載ってたらおかしいだろ！ハルナのドーナツ消費行動が！

また沖本氏の本からの引用ですが、こちらの本の82ページに載っている例は日本の株式収益率、英国の株式収益率、米国の株式収益率でGranger因果性の存在が認められていますね。
経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) | 沖本竜義 |本 | 通販 | Amazon

ジュン、ごめん、そのグレンジャー何とかの前に、こっちの本の145ページの、 SUR モデルっていうのがよくわからなかったんだけど。 $x_t$ と $y_t$ が互いに説明変数になっていなくて、ノイズどうしが相関をもってると SUR モデルっていうのか？でも $x_t$ と $y_t$ の説明変数のセットは共通なのに「見かけ上無相関な回帰モデル」と呼ぶのもよくわかんないし、そう呼んだとしてそれが何っていう…。

調べてみましょう。ウィキペディア（Seemingly unrelated regressions - Wikipedia）によると、SUR モデルというのは時系列の VAR モデルに限らず、一般の線形回帰モデルを複数束ねたものが、被説明変数が互いには説明変数に入っていない状態であって、ノイズが相関をもつかもしれない場合を指しているようですね。例えば以下のようなイメージです。

$\begin{cases} x = a_0 + a_1 x_1 + a_2 x_2 + a_3 x_3 + \varepsilon_x \\ y = b_0 + b_1 y_1 + b_2 y_2 + \varepsilon_y \end{cases}$

この2つの回帰式はぱっと見個別に推定できそうで、「無相関」に見えます。しかし、 $\varepsilon_x$ と $\varepsilon_y$ が無相関である保証はありません。なので、あくまで「見かけ上無相関」であり、それぞれの回帰式を個別に推定するのではなく2つの回帰式を同時にFGLS（実行可能一般化最小二乗法）で推定する必要があります（個別に最小2乗推定すると、パラメータの推定値の有効性が劣ります）。しかし、以下の2つの特殊なケースでは、個別に最小2乗推定してもFGLSと同じ結果が得られるんです。

$\varepsilon_x$ と $\varepsilon_y$ が無相関である場合（この場合は、見かけ上ではなく実際に無相関なので当然ですね）。
全ての回帰式で、右辺に登場する説明変数のセットが共通の場合。なんでこの場合に個別に最小2乗推定できるかは面倒なんで以下の論文の351ページまで見といてください。

Zellner, A. (1962). An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias. Journal of the American Statistical Association, 57(298), 348-368.

そして、ハヤトの言う通り、VARモデルはこの特殊なケースの後者なので、一般的に「見かけ上無相関な回帰モデル」というよりは、「説明変数が共通なのであまり見かけ上無相関にも見えないし、普通に無相関でないが、個別に最小2乗推定しても同時にFGLSしたときと同じ結果が得られる特殊なケース」と言った方がわかりやすいかもしれません。

なっが！

次にGranger因果性ですが、来期の春名さんのドーナツ消費量を、春名さんのドーナツ消費量のみを利用して予測した場合と、春名さんのドーナツ消費量に加えて春名さんの収入も利用して予測した場合で、後者の方が2乗誤差が小さくなるなら、春名さんの収入から春名さんのドーナツ消費量へのGranger因果性が存在するといいます。この検証にはF検定っぽいことをやります。

144ページに書いてあった「普段私たちが使っている因果という言葉の意味とは（後略）」ってのは？

それは146～147ページに説明がありますが、Granger因果性はその定義から「その変数があった方がよい予測ができる」でしかないです。Granger因果性の存在が確かめられても、通常の意味の因果関係、つまり、「春名さんの収入を変化させると、春名さんのドーナツ消費量も変化する」を抑えたわけにはならないんです。因果関係がなくてもGranger因果性の存在が確認されてしまうかもしれませんので。なので、複数時系列間の相互作用を別の切り口で調べる方法として、147ページにインパルス応答が挙げられていますね。

そういうことか。…161ページからは、ARCH・GARCH…また新しいのが出てきた…。

ARCHはノイズの分散を過去のノイズの大きさで回帰するようなモデルですね。ここまでずっとノイズ項の分散は $\sigma^2$ で固定でしたが、もはやこれも時間変化するということです。GARCHはノイズの分散を過去のノイズの大きさ及び過去のノイズの分散でも回帰します。過去のノイズの正負を考慮した拡張モデルも紹介されていますね。

その GJR モデルって実質、過去のノイズの大きさへの回帰係数を、過去のノイズの正負によって別の値に切り替えてるんだよな。それだと 0 で不連続になるし、結構適当な対応だなって思うんだけど、これでも強力ってことなのかな。…ともかくやっと本の前半終わった！！

（ノート7があれば）つづく