MCMCサンプリング練習

概要

・前回の最尤推定の練習からの続き
・二項分布のパラメータ最尤推定
・参考文献(後述)の二項分布のパラメータ「ふらふら最尤推定」
・同じく参考文献の二項分布のパラメータをMCMCサンプリング(メトロポリス法)

まずは普通に最尤推定

ということで前回は正規分布のパラメータ推定をoptimメソッドを用いて行いましたが,
正規分布はパラメータが2つなので,もっと単純にということで二項分布のパラメータ推定

今回は最適化アルゴリズムとして,準ニュートン法を使用しています.
実行すると

推定結果は0.44ということで,真値(0.45)が推定されてそうですというのが最尤推定

ふらふら最尤推定

参考文献P.176からの抜粋ですが,ふらふら最尤推定とは

1. パラメータ\(q\)の初期値を選ぶ
2. \(q\)を増やすか減らすかをランダムに決める(新しく選んだ\(q\)の値を\(q^新\)としましょう)
3. \(q^新\)において尤度が大きくなる(あてはまりが良くなる)なら\(q\)の値を\(q^新\)に変更する

この手順を\(q\)の値が変化しなくなるまで実行するものとなります.

Rでの再帰の練習もかねてfor文は使わずに書いています(後々面倒なことになるのですが…)
また,最尤推定の時に使った対数尤度を計算する関数Lも使いまわしています.
これを実行すると,

131229_1

というような感じに,\(q\)が徐々に収束していきます.
最初の条件分岐のところでランダムにどちらに進むか(尤度が改善される方向か,されない方向か)を決めるので,
ステップを進むにつれて段々と\(q\)の値が更新されていきます.

MCMCサンプリング

上記のふらふら最尤推定のアルゴリズムを少し改変したのがメトロポリス法によるMCMCサンプリングとなります(らしい)
上述したアルゴリズムに4番目が追加されて以下のようになります.

1. パラメータ\(q\)の初期値を選ぶ
2. \(q\)を増やすか減らすかをランダムに決める(新しく選んだ\(q\)の値を\(q^新\)としましょう)
3. \(q^新\)において尤度が大きくなる(あてはまりが良くなる)なら\(q\)の値を\(q^新\)に変更する
4. \(q^新\)で尤度が小さくなる(あてはまりが悪くなる)場合であっても,確率\(r\)で\(q\)の値を変更する

ここで確率\(r\)は尤度比
$$r = \frac{L\left( q^新 \right)}{L\left( q \right)}$$
で与えられます.

これをRで実装すると

となります.尤度が改善されない方向に進もうとした場合,尤度比(likelihood ratio)を計算して,
尤度比の確率で,尤度が改善されない方向にも\(q\)を更新していきます.

131229_2

初期値0.05から始めた場合も何やら0.45周辺でサンプリングされてるっぽいですね!!
試しに初期値を色々動かした場合にも同じようにサンプリングされるかというと,

131229_3

ということで大体100ステップ以降くらいからはqの定常分布になってるっぽいですね!何かすごい!

余談

で,

調子に乗って10万ステップくらいで計算しようとしたら

どうやら再帰で計算しているため,関数のネストが深くなりすぎて処理できなくなっている様子.
教えてジェネラル!

参考

Amazon.co.jp: データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学): 久保 拓弥: 本
http://www.amazon.co.jp/dp/400006973X

最適化 最尤推定の練習

概要

・Rの最適化メソッドoptimの練習
・最尤法の復習も兼ねて正規分布のパラメータ推定

Rコード

結果

結果としては,$parに推定結果,$valueに目的関数の値などが入るようです.

最適化アルゴリズムとしてはデフォルトではNelder-Mead法が使われる模様.
面倒な最適化も簡単に,そうoptimならね.

参考

関数の最大・最小化 – RjpWiki
http://www.okada.jp.org/RWiki/?%B4%D8%BF%F4%A4%CE%BA%C7%C2%E7%A1%A6%BA%C7%BE%AE%B2%BD#s8619974
ここの「負の二項分布のパラメータの最尤推定」など

セントルイス連銀のデータをpython経由で取得する

概要

・セントルイス連銀が提供しているFRED Economic Dataはデータが豊富なので欲しい
・PythonでFRED APIを叩いてデータを取得する

ユーザー登録

まずはユーザー登録を行い,apikeyを取得する必要がある様子

Register User Account – St. Louis Fed
http://research.stlouisfed.org/useraccount/register/step1

データ取得に関するドキュメント
St. Louis Fed Web Services: fred/series/observations
http://api.stlouisfed.org/docs/fred/series_observations.html

データ取得方法

APIが公開されているので例えば
http://api.stlouisfed.org/fred/series/observations?series_id=PAYEMS&api_key=[YOUR_APIKEY]&file_type=json&realtime_start=2000-01-01

といったURLにアクセスすればデータを取得することが出来ます.
(ユーザー登録後に取得したapikeyを[YOUR_APIKEY]のところに指定する必要あり)
ちなみに上記URLではPAYEMSという系列(アメリカの非農業部門雇用者数)を取得しています.
詳しくは上述したドキュメントを参照ということで(― ―)

ファイル・タイプとしては,
・XML
・JSON
・TXT
・CSV
といった方式が対応しているようです.今回はJSONを使用.

Pythonコード

まずはURLクローラということで,getPageText関数をパクじゃなくて参考にして作成します.
変数responseは通常str型で返ってきますので,ここではJSON形式に変換して返すようにしています.

これを実行すると

といった具合に値が返ってきます.
なお,これは非農業部門雇用者数となっておりますので,通常着目される非農業部門雇用者数の増加数自体は
136765と136562の差分をとって203となります.単位が1000人となっておりますので,
10月~11月に非農業部門雇用者数は203,000人増加したということですね.
この値は米国労働省のHPからも確認できます.

U.S. Bureau of Labor Statistics
http://www.bls.gov/

131221

サイトの右下のPayroll Employmentに +203,000とありそれっぽい値が取れてそうということが分かります.

参考

プログラマはWebページをデータとして扱う — ExSoft
http://www.exsoft.net/blog/entry/v8eik8

ウェブで数式を載せる方法 – MathJaxの導入

ウェブでの数式の表現について

ブログで数式を載せたい場合,何種類か方法があると思います;
・画像(texclip: http://maru.bonyari.jp/texclip/texclip.php)を作成
・x^2など代替表現
・画像生成サービス(よく知らないけどプラグインか何かを入れて表示させるやつ)

が,どれもしっくりこない
そんなときにとあるサイト様を見ていたらきれいに数式が表示されているではありませんか!
ということでやり方をパクではなくて参考にさせて頂きました.

MathJax

MathJax
http://www.mathjax.org/

どうやらJavaScriptベースで,一文字一文字ポジションやフォントサイズなどを駆使して表示させているっぽい(すごい人力!)
しかもLaTeX表現が使える!
かつクソIE6まで対応しているという親切っぷり!

いやはやディベロッパーの汗と涙がうかがい知れる代物なのでありますのでこれは敬意を表する意味でも是非導入しようという運びで使わせていただこうかなと思います.

導入

まずはプラグインの読み込み

とはいいつつも読み込んで終了なのであります.

あとは適当に数式を入力してみます(Texは数式を$で囲みましたが,MathJaxの場合は$$で囲むようです.ただこれも設定が変更できるようですが.)

$$ \def\rd{{\partial}} \frac{\rd V}{\rd t} + \frac{1}{2}\sigma^2 S^2 \frac{\rd^2 V}{\rd S^2} + rS \frac{\rd V}{\rd S} – r V = 0 $$

このようにBlack-Scholes微分方程式なんかも綺麗にかけると!
どうやら\defによる定義も使えるようでほぼTex文章と同じように使うことが出来そうな様子,という自己満記事でしたm(__)m

非定常過程同士の見せかけの回帰

概要

・非定常過程(eg. ランダムウォーク)同士を回帰させると有意なモデルができてしまう
・非定常過程同士の線形結合は基本的に非定常過程
・金融経済データのほとんどはランダムウォークなので気を付けようというお話

非定常データ同士の単回帰

適当な時系列データを作りまして

単回帰で回帰結果を表示させますと

Screen Shot 2013-12-07 at 2.19.36 PM

というような結果が帰ってきまして、切片項、説明変数共に0.1%の水準で有意という結果となりました。
決定係数も0.46ということで非常に約半分の変動が説明されているという結果に
(ランダムウォークによって回帰結果は異なるので上は例)

からくり

基本的に非定常データ同士を回帰させるとt値が発散して、決定係数も1に収束するとか。
この辺はハミルトン先生に聞くのが良さそうなのであるけどその辺は(´ `)y-~
(「基本的に」というのは、残差項が非定常過程になる場合。
残差項が定常だった場合は共和分という関係になるのですがそれはまた別のお話)

で、原因としては、単回帰における残差項が非定常となるためと。
定常かどうかを判定するためにadf検定を行ないます(要tseriesパッケージ)

Screen Shot 2013-12-07 at 2.22.53 PM

Screen Shot 2013-12-07 at 2.23.00 PM

まず元データであるxとyが非定常であるかどうかを判定。
p値が大きいので帰無仮説が採択されます。

同様に、単回帰における残差項も判定してみると
(単回帰の残差項のデータ自体は$residから取れるようです)

Screen Shot 2013-12-07 at 2.33.46 PM

こちらも非定常と判定されます。
ということで非定常過程同士の見せかけの回帰でしたというお話でした。
見せかけでない回帰結果も計算できるようですが、本日は以上m(__)m

参考

・経済・ファイナンスデータの計量時系列分析 (統計ライブラリー): 沖本 竜義: 本 – http://www.amazon.co.jp/dp/4254127928