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

概要

  • パーティクルフィルタ・平滑化を実装し、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

ExcelとPython(とR)を連携させる – ExcelPython編

概要

  • ExcelPythonを使うとエクセルのビルドイン関数のようにpythonの関数を呼び出すことができる
  • 当然numpyにも対応。ただしpythonのクラスには対応していない(?)
  • VBAを1行も書く必要がない。xlwingsではVBAに、呼び出したいpython関数を書かなければならなかった。
  • VBAを駆逐可能という代物

ExcelPythonのインストール

まずは下記サイトダウンロードし、インストールします。今回はv2.0.8を使用しています。

ExcelPython – Write Excel UDFs in Python!
http://ericremoreynolds.github.io/excelpython/

その後エクセルを起動させるとExcelPythonのリボンが追加されています。
150516_1

Pythonコード

ここでは前回の記事で使用した関数を流用します。
ポイントとしてはハイライトしていますが、xlpythonをロードするのと、エクセルから呼び出し可能にしたい関数の前に@xlfuncを書くことです。
あとはこれを使用したいエクセルシートと同階層に保存します。

簡単に説明するとgetDataFromQuandlは指定したティッカーのデータをQuandl経由でロードします。
garchはGARCH(1,1)のボラティリティを返します。

Excelでpythonコードをリンクさせる

といってもそこまで難しいことはなく

150516_2

  1. エクセルシートをマクロ使用可な形式(xlsm)で保存する。その際に先ほど保存した.pyファイルと同じ名前にする必要があります。
  2. ExcelPythonリボンにある「Setup ExcelPython」をクリック
  3. 「Import Python UDFs」をクリック

通常何も手を加えていないエクセルの場合以下の様なエラーが出ると思います。

150516_3

その際は
「Excelのオプション」→「セキュリティセンター」→「セキュリティセンターの設定」→「マクロの設定」→「VBAプロジェクト オブジェクト モデルへのアクセスを信頼する」のチェックをオンにします。

「Import Python UDFs」のクリックが完了すると、先ほど@xlfuncを指定した関数がエクセルから呼び出すことが可能です。

ExcelPythonを使ってみる

まずはaddNumを使ってみましょう。関数を使用しようとするとちゃんとサジェストにも表示されるようになります。
150516_4

Ctrl+Shift+Aの引数表示も対応しています。
150516_5

ということでめでたく計算することが出来ました。
150516_6

ExcelPythonを使ってみる – 返り値が行列編

どうやらExcelPythonはPandasのDataFrameには対応していないので行列に変換する必要があります。
その変換はto_records()で行っています。こうすることでPandasのindexも含めて返すことが可能になります。

先ほど定義したgetDataFromQuandlを使うと以下。1列目の日付は見やすいようにエクセル上でフォーマットを変更しています。

150516_7

ExcelPythonを使ってみる – インプットが行列編

エクセルからの入力が行列の場合、@xlargを使ってその引数のサイズを指定する必要があります。

縦ベクトルで返したいのでPython上で返すときにサイズを変換しています。

150516_8

ということでGARCHボラティリティの計算もすることができました。

感想

xlwingsとの比較という点ですと、
長所

  • マクロを書かなくて良い
  • ビルドイン関数のように呼び出すことができる

短所

  • Pandasに対応していない
  • 行列を返す場合の処理が面倒(予めエクセル上の範囲を指定する必要がある)

という感じでしょうか。

参考

ExcelPython – Write Excel UDFs in Python!
http://ericremoreynolds.github.io/excelpython/

excelpython/docs at master · ericremoreynolds/excelpython · GitHub
https://github.com/ericremoreynolds/excelpython/tree/master/docs

xlwingsを使ってExcelとPython(とR)を連携させる

概要

  • xlwingsというpythonパッケージを使うことでExcelからPython関数を呼ぶ
  • さらにRのPypeRパッケージを使うことでPythonからR関数を呼ぶ
  • Quandlから日経平均を取得し、そこからGarchモデルを計算してみる

xlwingsを用いたPythonとExcelの連携

公式HPとドキュメントが非常にわかりやすいのでそこに従うだけで連携は簡単にできます。(URLは記事一番下にあります)
xlwingsのダウンロードはpipがある環境では

で行うことができます。(この記事では0.3.5を使用しています。)
今使っているエクセルシートとPythonをリンクさせるには、1つモジュールをロードさせる必要があります。

エクセルのVBエディター(alt + F11)から「ファイル」→「ファイルのインポート」→「xlwings.bas」を選択します。
xlwings.basの保存場所はpython上で

とすることで確認できます。xlwings.basをインポートすることで標準モジュール内にxlwingsというモジュールが追加されます。

使ってみる

Pythonコード

コードを使う際にはQuandlのトークンを自分のものに置き換えてくださいね。
pypeRの使い方に関しては手前味噌ですがhttp://nekopuni.holy.jp/?p=1345に書いています。

VBA

VBAに関しては適当なモジュールを追加し、そこに以下のように記述しました。

Excel

適当に以下の様な感じのシートを作ります。
Generateボタンには上述したGarchマクロが割り当てられています。
このボタンをクリックするとGarchマクロからPython上のxl_getDataFromQuandl()関数が実行され、
日経平均とGarchボラティリティが返されます。
excel_ann

xlwingsはnumpyやpandasに対応しているので日付の処理などもスムーズに行なうことができそうです。
A列に表示されている日付はQuandlから返されたPandasをそのままExcelに貼り付けていますが、ちゃんとExcel形式の日付に変換されていますね。

これを使えば例えばRの機械学習モジュールなどもPython経由で呼び出すことが出来ますね。
夢が広がりングです。

参考

xlwings – Replace Excel VBA with Python!
http://xlwings.org/

xlwings – Make Excel Fly!
https://media.readthedocs.org/pdf/fzumstein/latest/fzumstein.pdf