パーティクルフィルタのテスト

概要

・パーティクルフィルタ(aka 粒子フィルタ or 逐次モンテカルロ法)のテスト
・とりあえず特定化されたモデルを対象
・勉強したことから成敗するためという名目で写経状態。ちなみにこの記事の全内容は参考文献に書いてあります。
・パーティクルフィルタスゲー

モチベーション

そもそもなんでパーティクルフィルタをやる必要があるかというところですが、一言で言うと「オンライン推定ができるから」ということになります。

・カルマンフィルタ:完全に特定化されたモデルの場合はカルマンフィルタでオンライン推定が可能。ただし通常は未知パラメータがモデルに含まれるため推定が必要
・MCMC:無いならパラメータを弾けば良いじゃないとなりますが計算が重い。しかも新たな観測値を得られたら、最新の状態を推定するために1からモンテカルロをぶん回す必要がある。そんな悠長なことは言っていられない。

未知パラメータを含みつつも、カルマンフィルタのように状態を逐次的に更新していくアルゴリズムがパーティクルフィルタというクラスに属する手法と。

モデルの超概要(という名の自分の中での理解)

結局求めたいのは状態\( \theta_t \)の期待値(と標準偏差くらい)なわけでして、

$$
E\left( \theta_t \right) = \sum_{i=1}^N w^{\left( i \right)} \delta_{\theta_t^{\left( i \right)}}
$$

という感じで、確率x値 という期待値の計算のように分解して状態を推定していきます。
ここで、\( w^{\left( i \right)} \)は重み、\( \delta_{\theta_t^{\left( i \right)}} \)が密度という感じで、イメージ的には各粒子にこの重みと密度を持たせてるという理解ですな。

で、あとは各粒子の重みと密度を新たな観測値が得られる度に更新していくわけですが、今回の記事では最適重点核(optimal importance kernel)という更新方法を使っています。
また、粒子を更新していくと次第に外れていく不良粒子が現れてきます。そのときはリサンプリングというプロセスを経由して、粒子を新たに生成させます。
最適重点核やリサンプリングの詳しい方法は参考文献参照ということで^^;;

Rコード

これはひどい写経だ!!

変数名まで一緒という手抜き使用^^;;;
まず下図が\( \theta_t \)の推定結果。実際の観測値、カルマンフィルタ、パーティクルフィルタの結果です。
カルマンフィルタとほぼ近い値が推定されています。

140627_1

そして下図が粒子数1000の場合のフィルタリング分布の標準偏差。パーティクルフィルタでの推定は若干ブレがあります。
ただしこれは粒子数を増やせばある程度精度が改善されます。

140627_2

粒子数10000の場合は以下の様な感じ
140627_3

参考文献

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

ドル円のボラティリティをGARCHで推定

概要

・GARCHでドル円のボラティリティを推定
・tseriesのgarch関数を使えるようになるまで

モデル

GARCHを何となく把握する必要があるような気がしてきたので概観とRの使い方だけメモ

リターン系列\( r_t \)をAR(m)モデルに適用した場合の残差系列\( \epsilon_t \)としたとき、GARCH(p,q)は以下。
$$
\epsilon_t \sim N \left( 0, \sigma_t^2 \right) \\
\epsilon_t = \sigma_t z_t \\
z_t \sim N \left(0, 1 \right) \\
\sigma_t^2 = w + \sum_{i=0}^p \beta_i \sigma_{t-i}^2 + \sum_{j=1}^q \alpha_j \epsilon_{t-j}^2
$$

2-3行目が解釈しづらかったのですが、\( \epsilon_t \)が期待リターンからの実現した乖離。一方で推定したいボラティリティ\( \sigma_t \)は実際にどの程度平均からバラつき得るか( \( \epsilon_t \)の分布の標準偏差)という感じですかねえ。

Rでの使い方

大きくtseriesパッケージとfSeriesパッケージが主流のようですが、fSeriesは安定性に欠けているとのコメントが散見されたのでここは大人しくtseriesで。

データはFRED先生からドル円レートを落としてきて使います。

Japan / U.S. Foreign Exchange Rate – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/DEXJPUS

適当にfxという変数に格納し、とりあえず値を概観。

140622_1

FREDから取れるだけデータを取ると1971年から取得できるのですが、今回は1973年4月からのデータを使います。
(1973年2-3月に日本も変動相場制に移行したため by Wikipedia)

garch関数を使うことで、系列に対してGARCHを適用してくれます。デフォルトはGARCH(1,1)。
モデル選択は面倒なので今回は割愛。
summaryの出力結果は以下

Residualsが残差系列の要約、
CoefficientsがGARCHの推定されたパラメータ、定数項を含めた3つのパラメータはすべて有意となっています。
Diagnostic Testsでは正規性のテストとしてJarque Bera検定、系列相関のテストとしてBox-Ljung検定の結果も表示されます。
Jarque Bera検定は帰無仮説「残差系列が正規分布に従う」となります。p値が非常に小さいことから帰無仮説は棄却され、残差系列は正規分布に従っていないと言えそうです。
Box-Ljung検定は帰無仮説「残差の2乗に1次の系列相関がある」となります。こちらは有意水準5%のもとで帰無仮説が棄却されず、系列相関が無いとは言えないという結果に。

ボラティリティ推定結果

\( \sigma_t \)は日次ボラティリティで、少し見づらいためfx.sigmaに格納するときに年率換算しています。
またこの値はfitted.valuesという変数に格納されています。

140622_2

平均では年率換算で10%をうろちょろし、時々ジャンプして時間とともに落ち着いてという感じですね。
直近(2014年5月末)はニュースなどになっていますように非常に落ち着いた環境であります。
まあこの推定結果を見ると1970年代とかはボラティリティが落ち着いたまま何年も過ごしている期間がありますし、足元でボラティリティが下がっているからといって近々ボラティリティが上昇するという見方はちょっと違うんじゃないんですかという感じですかね。まあいつかは上がるわけですが。

不明点

以下個人的に謎だったメモ
・残差系列の計算方法:教科書的には残差系列はAR(m)過程に適用した残差系列となっておりますが、これが今回使用したgarch関数ではどのように計算されているのかが謎(ソース追いかければ分かりますががが)
・というかそもそもリターン系列をAR過程にぶち込むという発想自体が違和感。特に為替レートなんてゼロサムなもんで期待リターンってあるんですかねえという疑問。普通にリターン2乗を残差系列としても良いような。

参考文献

Amazon.co.jp: ボラティリティ変動モデル (シリーズ 現代金融工学): 渡部 敏明, 木島 正明: 本
http://www.amazon.co.jp/dp/4254275048

Package “tseries”
http://cran.r-project.org/web/packages/tseries/tseries.pdf

Time series analysis using R: Estimation of stationary and non-stationary models(日本語)
http://www.geocities.jp/atsmatsumoto/

日次データから月末のみを抽出する

概要

・xts型の日次データから月末データのみを抽出するメモ
・もっと効率よい方法を募集中です → あった(補足参照) → まだあった(補足2参照)
・応用すれば月初や、月末1日前なども取得可能かと

データ

例えば以下の様な日次データがxts型として格納されているとする。
以下は10年スワップレート(FREDより)

方法

一行で書くと以下

順を追って書くと以下

1行目でまず月毎に分解します。tmp1には月ごとに分解されたxtsのリストが格納されます。
2行目で、各リストの要素に対して、tailを使って一番下の行(月末データ)を抽出します。
3行目でリストをマージします。

すると上のように月末データのみを保有するxtsが返ってきます。
月初を取得する場合はtailを使っているところをheadに。

補足(追記)

この記事を書いた数分後に某鎌倉方面から「まだそんな方法でやってるのプゲラその方法が許されるのは小学生までよねーゲラゲラ」というコメントを拝受しました

日次データから月末のみを抽出する – My Life as a Mock Quant
http://d.hatena.ne.jp/teramonagi/20140614/1402745254

to.monthlyで四本値が取得できるので

とすることでそれぞれ月初、月末を取得することができます。さんざん悩んでいたワタクスの時間はものの数分で虚無へと帰りました。

ま、まあこれだと月末1日前とか取得できないから…(震え声

補足2(追記)

同じく某鎌倉方面から完膚なきまでに追撃を受けましたのでここに追記させていただきます。

とありますように、endpoints関数で取得したいindex値を取得し、日次xtsから該当の行だけ抽出すればいいんじゃねーのという方針

上から月末、月初、月末1日前のデータを抽出。
月初のデータならto.monthly使った方が楽なような気がしますが、endpointsを使えばほぼどのような条件でも抽出できそうですね。

米国金利を主成分分析+ファクターの推移

概要

・米国金利(今回はスワップレートを使用)の推移
・主成分分析結果 よく言われるように3つの主成分で大体の変動は説明される
・どの主成分がいつ効いているか(ファクター)の推移

米国金利(スワップレート)の推移

なんで米債の金利じゃないんですかと聞かれたらFRED先生がなぜかスワップレートしか提供してくれていないからです。
米債の金利に関してはFEDか財務省あたりが公表しているとかしていないかという話も聞きますがデータの取得がめんど(― ―)
まあ金利スワップも米債利回りも同じようなもんだからということで、データは以下からいつものように取得。

Interest Rate Swaps – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/categories/32299

xts型としてswpという変数に格納したとします。
まずはplotで動きを概観

140608_1

plotのデフォルトの色だと黄色とかが見づらかったので初めてRColorBrewerというカラーパレット作成パッケージ(?)を導入。
2000年7月とかは7%台でほぼフラットなイールドカーブ。
そこからITバブル崩壊とともに徐々に右肩上がりなイールドカーブに、その後連続利上げで短中期を中心に上昇し2006~2007年ではまた5%台あたりでフラットに。
その後リーマン・ショックと続く事実上のゼロ金利政策、QEなどでカーブは水準を切り下げぺっちゃんこ相場に。
2013年5月のバーナンキ前議長のtaper発言でまた金利が上昇するもそこまで上がらず今に至るという感じですな。

日本と異なり非常にダイナミズムあふれるイールドカーブで何よりです。

各スワップレートの水準(変動ではない)に対する標準偏差を見ると

140608_2

スワップレートの水準推移を見たとおり、順イールドからフラットまで短期はどったんばったんしますんで、それを反映する形で短期になればなるほど標準偏差は大きくなる傾向に。

主成分分析とファクターローディング

本当はスワップレートの変動に対して主成分分析を行うほうがベターなんでしょうが、後にみるようにファクターの推移も見たいと思いまして、その関係上、水準に対して主成分分析をします。
また、先ほど水準の標準偏差を見ると短期のスワップレートになるほど標準偏差が大きくなる傾向にありました。
こういったデータに対して主成分分析を行う場合、相関係数行列に対して主成分分析を行ったほうがベターだと思うのですが、それは正規化(平均0、分散1)のスワップレートに対して主成分分析を行っているということでして、それもそれでファクターの推移として見た時にイマイチ解釈しづらいというか結局見たいのってスワップレートの水準ですしお寿司ということで泣く泣く水準の分散共分散行列に対して主成分分析をします。

princompのcor=Tとすると、相関係数行列に対して主成分分析を行いますが、今回はFにして分散共分散行列を使います。
summaryを見ると以下

ということで累積寄与率(Cumulative proportion)を見ると、3つの主成分までで変動の99.95%が説明されるという結果になります。
もちろん8つ全ての主成分を使えば100%の変動を記述することもできます。

ファクターローディング(固有ベクトル)は以下。

140608_3

ということで、黒線=レベル、赤線=スロープ、緑線=カーベチャーの3成分が出てきましたとさ。

ファクター推移

ということで以上3つのファクターローディング(軸に相当)で、スワップレートというのは表現できるということがわかりました。
で、じゃあこの3つのファクターが過去どのように推移してきたのかよというのを見たのが以下。

coefには日次での、ファクターローディングに対するスワップレート水準のOLSの回帰係数をぶち込んでいます。
その後見やすさも兼ねてxts型に変換。プロットすると以下

140608_4

黒線はレベルファクターということで、一番最初に見ましたスワップレートの推移と同じような動きをしています。
(第一主成分はマイナス方向に出ていましたので、plotするときに-1をかけて直感に合うようにしています)
一方赤線がスロープファクターということで、イールドカーブが右肩上がりだったITバブル崩壊後には比較的高い値を出している一方、利上げを繰り返しイールドカーブがフラットになった2005-2006年あたりはファクターの値が小さくなっています。
最後に緑線がカーべチャーファクターで、えーこれはなんとも解釈しづらい。

ファクター推移おまけ

140608_5
2010年から2012年の3ファクターの推移を拡大したのが上図。
一方、QE2が2010年11月~2011年6月で、オペレーションツイストが2011年10月~2012年6月。
ということで、政策発表前および発動中の期間に関して、特にレベルとスロープに対して押し下げるように働いたと見ることもできそうですね。

おまけ2

なんでいきなり主成分分析なんかやり出したかと言いますと、このファクターの推移ってのをDLMにぶちこめないかなあと妄想するわけでして、主成分同士って直交するから独立で良いし。
ボラティリティ(VIXなど)の期間構造を主成分分析した場合も同じようなファクターローディングが現れるようですし色々やってみたら面白いんじゃないかと思います。

ま、結局予測力という観点からだと移動平均と大して変わらないっていう結果が出てきそうな気もしますけどね!!

参考文献

Amazon.co.jp: 主成分分析―講座 情報をよむ統計学〈8〉 (講座情報をよむ統計学 8): 上田 尚一: 本
http://www.amazon.co.jp/dp/4254127782

1961年から2011年、オペレーション・ツイストの差はなにか
http://www.mizuho-ri.co.jp/publication/research/pdf/market-insight/MI110922.pdf

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