Cross Entropy Minimisation

概要

  • エントロピー差が最小になるような最も近い分布を求める案件
  • 平均2.2の正規分布から平均2.3の分布を求める

Rコード

参考サイト1にMatlabコードがあるのでそれをR用に書き写した

結果

期待値2.2と2.3がそれぞれ求められているのであっているっぽい。

PDFをプロットすると以下の様な感じ。参考サイト1と同じに見えるっぽい?
150614

参考

Tutorial Demo on Cross Entropy minimization (I) [Matlab Octave]
http://www.quantcode.com/modules/mydownloads/singlefile.php?lid=478

Forecasting Using Relative Entropy
https://www.frbatlanta.org/-/media/Documents/filelegacydocs/wp0222.pdf
(元論文。あまり読んでいない。)

未知パラメータがある場合のパーティクルフィルタ

概要

  • 未知パラメータ(ノイズの分散など)を含む場合のパーティクルフィルタの実装
  • システムノイズが正規分布の場合とコーシー分布の場合の比較

システムノイズが正規分布の場合

詳しい理論に関しては参考文献のP.225辺りからを御覧ください。

対象データに関しては前回の記事と同様のものを使っております。前回との違いとしては、システムノイズと観測ノイズの分散が同時に推定されているという点ですね。
また、今回は重点遷移密度を得るためにAnxiliary Particle Filter(補助粒子フィルタ)を使っています。これは、未知パラメータを含んでいるからという理由と、他の分布への一般化も視野に入れているからという理由です。

また参考文献中では実装されていませんが、多項リサンプリングも行なっています。

状態の推定結果

粒子数は1万個です。以下のように、若干ずれる箇所はありますが、カルマンフィルタとほぼ同程度の推定を行なっていることが分かります。
150607_1

システム分散の推定結果

システム分散の真値は1ですが、まあまあという感じですかね。
赤の破線はそれぞれ25%タイル値、75%タイル値を表しています。
150607_2

観測分散の推定結果

こちらの真値は2ですが、まあこんな感じという感じなんですかね(― ―)
150607_3

システムノイズがコーシー分布の場合

ほとんど正規分布の場合と変わりません。正規分布のコードをコピーして再利用しているので、変数名などが適切ではないですがそこは目をつむってください(― ―)

今回推定対象としたデータ列は以下。コーシー分布に従う状態に、正規分布を上乗せしています。
状態がコーシー分布なので、正規分布よりもドラスティックなジャンプが時々発生します。
150607_0

状態の推定結果

粒子数は100万個。青線は状態の真値で、赤線がパーティクルフィルタで推定された状態です。
概ね推定されているような気がしなくもないですがやはり精度はかなり悪くなりますね。
100万個でこれですからね。MBPで走らせるとかなり時間がかかります(― ―)
150607_1

システムノイズのscaleの推定結果

真値は0.05ですがこんなもんなんですかね。時々ゼロに張り付いたりしていますが。
150607_2

観測分散の推定結果

こちらはまあいい感じですかね。真値は1.5です。
150607_3

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

パーティクルフィルタ・平滑化の実装

概要

  • パーティクルフィルタ・平滑化を実装し、DLMパッケージによる結果との比較を行なう
  • 平滑化の計算が合っているのか微妙。教えて偉い人!

Rコード

サンプル時系列の生成

参考文献[1]のP.219を参考に、dlmパッケージを使って適当な時系列を作成します。コレに対してパーティクルフィルタ・平滑化を適用させます。
今回は単純な線形ガウシアンモデル。パラメータも特定させておきます。

答え合わせように、フィルタ、平滑化も先に計算しておきます(各々modFilt, modSmoothに格納)。

パーティクルフィルタ本体

重点遷移密度の選び方について、DLMであれば最適重点核として\( y_t \)が与えられたもとでの事後分布を選ぶのが良いのですが、今回は事前分布を使っています。
非ガウス分布の実装も視野に入れているのでこういう仕様にしています。パッと見た雰囲気では、この事後分布を一般的に書くのはできないんですよね?(たぶん)

さらに、参考文献[1]のP.216に紹介されている多項リサンプリングを入れています。
これは少数の粒子が大きい重みをしめてしまうことによるモンテカルロ近似の劣化を防ぐための措置です。

パーティクルフィルタのメインの実装に関しては参考文献[2]のP.214のアルゴリズムをなぞっているつもりです。

プロットなど

y.pfに上述のパーティクルフィルタ関数の実行結果を格納しています。

結果

上記のプロットを行なったものが以下

フィルタ分布

フィルタに関してはカルマンフィルタと同じような結果が出ていそう。
赤破線はパーティクルフィルタからもとめた1標準偏差です。
150530_1

フィルタ分布の標準偏差に関しては、パーティクルフィルタは精度が落ちていますね。参考文献[1]よりも明らかに悪そうです。これは重点遷移密度に事後分布を選んでいなからなのでしょうか…
150530_2

平滑化

平滑化ラグは2を使用。なんとなく近い分布が得られていそうな気がしますが、微妙に違う感じにも。。
150530_3

標準偏差は以下、フィルタ分布よりは精度が上がっているような気がしますがそれでもかなり暴れん坊将軍ですな。
150530_4

平滑化(ラグが20の場合)

参考文献[2]によると、ラグは20くらいを選ぶと良いと書いてあるのでそれを選んだ場合の結果。
ラグを含む分布の精度が極端に悪い。なぜでしょう。
150530_5

150530_6

今後の課題

  • そもそも平滑化あってるのか?
  • 非ガウス分布
  • 未知パラメータを含むパーティクルフィルタ

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

Amazon.co.jp| 時系列解析入門| 北川 源四郎| 本| 微積分・解析
http://www.amazon.co.jp/dp/4000054554

モンテカルロ・フィルタおよび平滑化について
http://www.ism.ac.jp/editsec/toukei/pdf/44-1-031.pdf

逐次モンテカルロ/(粒子|パーティクル|モンテカルロ)フィルタを実装してみた – My Life as a Mock Quant
http://d.hatena.ne.jp/teramonagi/20140525/1400996808

Particle Filter
http://daweb.ism.ac.jp/koza/koza2008/PF_Nakano20081030.pdf

主成分分析を使った簡単なRelative Value戦略 – 米スワップレート編

概要

  • 単に主成分分析を使って主なファクターで説明される部分と、説明されない部分を用いるという物
  • ワークするかワークしないかは保証しません

Rコード

図表

図1. まずはスワップレートの推移

特にコメントするまでもありませんが、日本のイールドカーブとは違ってダイナミズあふれるものとなっております。2000年代前半とリーマンショック直前まではイールドカーブはほぼフラット。
直近も利上げ期待に呼応する形で短い金利は上昇傾向です。
150524

図2. ファクターローディング

これもよく言われる話ですがイールドカーブの”変動”は3つのファクターで大方説明されます。
ただし今回はRelative value戦略を見たいので、金利”水準”を主成分分析対象としています。
結果は省いていますが、水準で主成分分析を行なった場合も、3つファクターで大方説明されるようです。
(というか1つめのファクターで96%が説明される模様。ほんまかいな。)
ということはイールドカーブ(の標準化したもの)は下図の線形和でだいたい説明されるということになりますね。
150524_2

図3. スコアの推移

第1ファクターローディングがマイナスなので見づらいですが、ここ15年間金利は基本的に低下していたため、第1スコアは上昇傾向となっています。で、直近は第2ファクターがマイナスの値をとっておりまして、要するにベアフラットニングみたいのを説明しているんですかね、よう知りませんが。
150524_3

図4. 5月20日時点でのイールドカーブ

out.Actualとなっているのが実際に観測されたイールドーカーブ、out.Estimatedとなっているのが3つのファクターにより説明される部分のみを抽出したイールドカーブとなります。
今回のRelative value戦略の思想としては、「3つのファクターで説明されない部分はノイズであり将来修正される」としていますので、この差が大きい部分をロングなりショートなりするという形になります。
150524_4

図5. イールドカーブの差分

図4の差分をプロットしたものが以下。Actual – Estimatedを図示していますので、プラスの値が割安(観測される金利が高い)、マイナスの値が割高(観測される金利が低い)となります。
今回はスワップなので3年固定を受けて、10年固定を払うという感じですかね。
150524_5

元本をデルタニュートラルになるように調整して、このイールドカーブの差分が解消されるタイミングでアンワインドするというのが理想的なシナリオでしょうか。

ということで簡単に主成分分析でRelative value戦略というのを考えることもできますが、この分析での問題点としては

  • ファクターで説明されない部分が修正されるとは言っていない
  • 仮に修正されるとしても、それにどのくらいの時間が必要かは分からない
  • 3年受け10年払いのポジションの場合、スティープニングのポジションを取ることになりますが、おそらく要らんファクターエクスポージャー(リスク)を持っている。

とかですかね。

今回は金利スワップなので個人は取引できないですが、例えば同じ業界内の株式などでも同様の分析をつかったRelative value戦略はできそうですね。ショートするのが大変かとは思いますが。
あとは個別銘柄リスクを負うことになりますね。。。
日経平均とTOPIXとかは分析していないので分かりませんが、差分が出てこなさそうなイメージ。コスト割れしそうですね。

為替などでも出来るんですかね。あんまり調べてないですが。
AUD対NZDとかSEK対NOKとかですかね。

コモディティだと今ホットな原油系ですかね。誰かお願いします。

メモ: ggplot2に使いやすさを求めるのは間違っているだろうか – 折れ線グラフ重ね書き編

概要

  • data.frameから折れ線グラフを1つのプロットに重ね書きする方法のメモ
  • reshape2パッケージのmelt関数で一回データフレームを変形させる必要がありそう
  • autoplot()を使えば超絶簡単に書けるとの有力情報を頂きました

データ例

対象としているデータはこのようなフォーマット

時系列データが縦に並んでいる。xts形式。
ちなみにデータの内容はFREDの米スワップレートをQuandl経由で取得したもの。まあそれはどうでもよい。

Rコード

ひと通りの流れを説明すると、ycはxts形式なので、yc.dfにデータフレーム形式として変換
ついでにその際に日付の系列をデータフレームに追加しておく。
次に、melt関数を使ってデータフレームを変形させ、yc.df.Lに格納。このyc.df.Lを描画対象とする。

150524

このやり方が比較的一般的に使えそうなのでとりあえずこの方法でやっていこうと思います。。。
1つのデータフレームから棒グラフと折れ線グラフを重ね書きしたいというケースも有ると思いますが、それをやろうとすると折れる骨がいくつあっても足りないと思うので次回に回したいと思います。

後日談

…というポストを書いたあと秒速で以下のようなご指摘を頂きました。

ということで今回のケースの場合、

と1行書くだけで所望のグラフを書くことができました。ありがとうございます。
…今までの苦労は…!!

参考

r – Plotting two variables as lines using ggplot2 on the same graph – Stack Overflow
http://stackoverflow.com/questions/3777174/plotting-two-variables-as-lines-using-ggplot2-on-the-same-graph