2023の18日目の記事です。1日遅れですが、枠が空いていたので飛び入り参加してみました。Stanで実装してみました。Stanの推定結果のプロットにはbayesplotとggplot2を用いています。モデルをStanで実装することで、TOPIXのボラティリティを推定しました。推定されたボラティリティは、2008年のリーマンショックと、2011年の東日本大震災、2020年の新型コロナウイルスによる市場の急落局面で非常に高まっていることを確認できました。ボラティリティが一度高まるとしばらくボラティリティが高い日が続く現象であるボラティリティ・クラスタリングも確認できました。ボラティリティの定式化ボラティリティとは、金融商品の値動きの変動の大きさを指すパラメータです。Risk)に使われるなど、金融工学の根幹をなす値です。以下は株式に限らず、例えば為替など、金融商品であれば何でもよいですが、株式とします。S_{t-1}$となります。これを対数収益率といいます。このとき、ボラティリティとは以下の$\sigma_t$、あるいは$\sigma_t^2$を指します1。$$$\sigma_t$と$\sigma_t^2$のどちらをボラティリティと呼ぶかは文献によりますが、以降$\sigma_t$をボラティリティと呼びます。この$\sigma_t$を推定する方法は大きく分けて三通りあります。一つ目の方法は、$\sigma_t$を過去一定期間の$r_t$の標準偏差とする方法です。r_{t}\}$の標準偏差とします。直近のボラティリティを重視したいなら$j$をより小さい値に設定します。この$t$を1日ずつずらしてローリング計算することでボラティリティの時系列を得ることができます。この方法はシンプルですが、標準偏差を計算している$t-j+1,\dots,t$の間は$\sigma_t$が一定であることを暗に仮定しています。実際はそうではありませんので、精緻に求めるなら残りの二つの方法のどちらかを用いることになります。二つ目の方法は、$\sigma_t$を統計的なモデルで定式化するものです。モデルと呼ばれるモデルのうち、シンプルな以下のモデルを推定することにします。$$ここで$y_t$は対数収益率$r_t$であり、$\exp(x_t/2)$がボラティリティ$\sigma_t$です。$x_t$は定常な過程と仮定し、$|\phi|<1$です。\sigma_t^2$がAR(1)モデルに従うことを意味します。このモデルは1本目の式が観測方程式、2本目の式が状態方程式の非線形な状態空間モデルですので、パーティクルフィルタのような非線形な状態空間モデルでも推定できるタイプのフィルタ系の手法か、Stanなどを用いてMCMCで推定することになります。最後に三つ目の方法としては、以上二つのように収益率$r_t$の時系列から$\sigma_t$を推定するのではなく、分単位のような細かい収益率データを用いて直接$\sigma_t$を求めるアプローチがあります。Volatility2と呼びます)になることが知られています。Stanコードを書いてrstanというRのパッケージから呼び出します。Rを用いているのはベイズモデリングは文献もパッケージもRが豊富なためです。あとこの記事のコードは自分で昔書いたコードを流用しているのですが、それがRだったからというのもあります。(RStan)で実装し、MCMCによってパラメータを推定してみます。TOPIXの日次の終値が必要なので用意します。2008/5/7~2023/12/18(3824営業日)のTOPIXの終値をJ-Quantsから取得しました。J-QuantsはJPXが個人向けにリリースしている株価などの金融データのAPIです。この記事のモデリングをするにはTOPIXの終値だけあればいいので、J-Quantsではなくとも構いません。証券会社に口座を開いていれば、証券会社が提供しているトレーディングツールからCSVでエクスポートできたりもします。データの取得コードJ-Quantsは、認証エンドポイントを2回POSTしてトークンを取得したら、そのトークンをBearer認証のヘッダに入れてTOPIXのエンドポイントをGETするだけのシンプルな作りです。"hoge"}のようなJSONをbodyに入れます。auto_unbox=TRUE)とします。["hoge"]}となってしまいます。Rではmail@example.comのような文字列は長さ1のcharacter型のベクトルであるため、そのまま長さ1のリストができてしまうからです。これを防ぐためにauto_unbox=TRUEを指定します。TOPIXの終値から対前日の対数収益率を計算します。100倍して%表記にします。最初の1日の収益率はNAになるので最初の1日を除いておきます。slice(-1)こんな感じのデータです。上に挙げたSVモデルをStanコードで書きます。}}以下、実装のポイントです。ある以下の分布を用いました。IGは逆ガンマ分布です。無情報事前分布ではなく、事前分布を書いてあげると収束しやすくなります。quantitiesのvolが今回求めたいボラティリティ$\sigma_t$です。MCMCの結果から$x$の事後分布を取り出してR上で$\exp(x/2)$を計算することで求めることもできますが、Stanで求めておきます。quantitiesのy_predは、各パラメータの事後分布から乱数で生成した$y_t$の値です。これはパラメータの推定自体には不要ですが、あとでモデルの事後診断に使用するために求めています。あとは高速化のためにできるだけベクトル化をします。このStanコードを”svmodel.stan”で保存して以下のRコードでキックすることでMCMCの推定を行います。i9-9900Kで実行しました。thin=50で8時間程度かかりました。iterをかなり大きくしてthinで間引かないと$\phi$などのパラメータに自己相関が残っていました。seed=1234)MCMCのチェックMCMCがうまくいっているかチェックします。MCMCが終わったら早速パラメータの値を見たいところですが、収束していないということがよくあるのでチェックしましょう5。収束していなければパラメータが正しく推定されていないということです。詳細に知りたい方は、例えばこちらの分かりやすい記事をご覧ください。diaryこのあたりの事後診断の可視化はbayesplotが超便利です。ggplot2ベースなので出力のグラフにggplot2の関数で軸などのレイアウトを微調整できるのもナイスです。MCMCが収束しているかの目安であり、全てのパラメータで1.1を下回ることが必要とされます。パラメータが多くて個別にプロットするのは難しいので、ヒストグラムで描きます。bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit))次に有効サンプルサイズ(n_eff)です。Nで割ったものが0.1以上であることが望ましいとされます。0.1を下回るような小さい値のパラメータは、MCMCのサンプリングで自己相関が残っていることを示唆します。自己相関が大きいとパラメータの事後分布の分散を正しく推定できません。こちらもよさそうですね。bayesplot::mcmc_neff_hist(bayesplot::neff_ratio(fit))次は自己相関です。サンプルに自己相関が残っている場合、サンプルに定常性がないということになります。全部描けないので以下の3つのパラメータに絞ります。phiとsigma_etaにちょっと自己相関が残っていそうですが、おおむねいい感じです。iterとthinを増やせばphiとsigma_etaの自己相関も消える気がします。"sigma_eta"))ちなみにbayesplot::mcmc_acf_bar()をはじめ、bayesplotの描画関数はparsでパラメータを指定し忘れると全パラメータ描画します。巨大なモデルだとこれでRStudioがクラッシュすることがありますので気を付けましょう。時系列の状態空間モデルはパラメータの数が多いのでクラッシュしがちです(今回はxがx[1]からx[3823]まである)。モデルをRDSファイルに保存し忘れた状態だと推定結果が消えて悲しいことになります。トレースプロットを見ます。こちらもパラメータを絞っています。これはchainごとのサンプルの推移で、線が混ざり合っていると初期値によらず同じ値に収束している=局所解に落ちていないことを示します。Rhatが大きいときはトレースプロットが混ざり合っていないので、Rhatと合わせてチェックします。"sigma_eta"))最後に、モデルが現実のデータをよく説明するように定式化されているなら、パラメータの事後分布から乱数を振って得られる目的変数(ここでは$y_t$)は、実際の$y_t$と同じような分布で得られるはずです。これを図示してみます。濃い線は実際のy_t(対数収益率)の分布、薄い線はモデルから生成した8000系列の対数収益率y_predのうち最初の10系列の分布です。薄い線はそれぞれの系列で10本あります。8000系列は描けないので10本に絞っています。濃い線と薄い線が大体重なっているので悪くなさそうです。rstan::extract(fit)$y_pred[1:10,])自分でモデルを一から組む場合は、モデルのパラメータはMCMCで正しく推定されていても、そもそもモデル自体が元のデータを全然説明できていない誤ったモデルであることがあります(SVモデルは幅広く用いられているモデルなので、今回ちゃんと当てはまるのはある意味当然なわけですが)。モデルが元のデータに当てはまっているかどうかを示してくれるのが上のプロットです。bayesplot::ppc_*系の関数はこのような事後予測の確認に役立つものが色々あって便利です。以上がMCMCの結果の基本的なチェック方法です。以下がパラメータの推定結果です。printが長くなるので一部のパラメータに絞っています。Rhat=1).推定したボラティリティを見てみましょう。rstan::extract()を用いてパラメータの事後分布のサンプルからvolの時系列を作ります。点推定値として事後中央値を用います。あわせて95%ベイズ信用区間も示したいので、中央値と2.5%タイル点と97.5%タイル点を取り出します。8000個のサンプルを小さい順に並び変えて50%と2.5%と97.5%のタイル点を取り出すことで得られます。rows以下のグラフの上はボラティリティ$\sigma_t$(赤線は$\sigma_{t}$の事後分布の中央値、青いバンドは95%ベイズ信用区間)、下はTOPIXの終値です。a$であれば、TOPIXの収益率の標準偏差がa[%]であることを意味します。ncol=1)2008年のリーマンショック、2011年の東日本大震災、2020年の新型コロナウイルスの市場急落局面でボラティリティが高まっていることが分かります。グラフはpatchworkで並べました。複数のggplotオブジェクトを綺麗に並べられて重宝するパッケージです。最後に$\phi$のパラメータの事後分布を見てみます。0.951-0.978)でした。pars="phi")このパラメータは過去のボラティリティがどの程度後を引くかを示すパラメータであり、1に近いということは持続性がかなりあるということを示します。ボラティリティが一度上昇するとしばらくボラティリティが高い日が続くということであり、この現象をボラティリティ・クラスタリングといいます。SVモデルを推定した論文をサーベイすると$\phi$の推定値は0.8から0.995までの値となっているという論文がありますが6、この先行研究と整合的です。おわりにStanを使うと柔軟にモデルを組んで解釈ができて楽しいですね。SVモデルの一番シンプルなものを推定してみましたが、SVモデルには収益率の裾の厚さを表現するために誤差項をt分布としたり、収益率がマイナスの日の方がボラティリティが高まる(ボラティリティの非対称性といいます)ことを表現したりするモデルなど、色々な発展形があります。これらをStanで組んでみてもよいでしょう。最後に株価データで時系列モデリングをやってみたいという方におすすめの書籍をご紹介します。現代金融工学)統計モデル的アプローチによるボラティリティモデルです。ボラティリティモデルをしっかり学ぶならこれです。MCMCによる推定も少し触れられています。コード例はなく数式展開の理論の本です。(統計ライブラリー)前半はカルマンフィルタの解説、後半は経済・金融データへのカルマンフィルタの適用事例の紹介です。前半はカルマンフィルタの導出を丁寧に書いていて理解しやすかったです。後半はカルマンフィルタでベータ値を推定したりペアトレーディングの銘柄を発掘したりと面白い事例が豊富です。コードはないですが読みやすいです。memo式の通り、$z_t$は正規分布とは指定していないので、$r_t$も正規分布以外の分布を仮定することができます。実証的には$r_t$は裾が厚く正規分布ではないとされます。この辺りの論文をご覧ください。Anubis.suzuna。