Stanで確率的ボラティリティ変動モデルの推定

Stan入門として確率的ボラティリティ変動モデルの推定をやってみた。前々回の投稿SVモデルでのオプションのプライシングでは潜在変数の推定結果を論文から引用していたが、今回はその推定を行う。

金融機関では、たとえばVaRを用いて保有資産がある額以下になる確率（分布の裾確率=tail riskにあたる）を一定以内に抑えるようにリスク管理しているが、その確率は資産価格のモデルの仮定に依存する。もっとも単純なランダムウォークモデルでは将来の価格の分布は正規分布となり、裾確率が過小評価されることで知られている。
デリバティブの価格はボラティリティの状態に大きく依存するので市場より正しく予測できてミスプライシングを発見できれば儲けられる。

正攻法(MCMC)で推定しようとすると死ぬほどめんどくさいが、Stanでやってみるとかなり簡単にできることがわかる。

【Stanとは】

モデルを記述してデータを与えるだけで勝手にベイズ推定してくれるソフト。推定にはHMC(ハイブリッドモンテカルロ)などを使っている。
似たソフトにWinBUGSやJAGSなどがある。時系列特化だとLiBiというソフトも。
ギブスサンプラーなどを自分で書かなくてよいのが利点。
収束が遅かったり、おかしな値にはまったりする可能性があるのが欠点。

【確率的ボラティリティ変動モデルとは】

ファイナンスで用いられるモデル。ボラティリティが常に一定であるランダムウォークとは異なり、ボラティリティ自体もランダムウォークするというモデル。
たとえば株式市場が長期にわたって不安定になる(=ボラティリティが高い状態は一定期間持続する)様子などがうまく記述できる。
確率的ボラティリティ変動モデルは非線形状態空間モデルというクラスに位置し推定が難しく、実務ではより推定が容易なGARCHなどが用いられてきたが、こちらはデータへの当てはまりが悪いことで知られている。
頻度主義的な推定法も存在するが、推定効率が悪いのでベイズ推定が主流。その場合はMCMCを用いて、「潜在パラメーターを既知として状態パラメーターをサンプリング→状態パラメーターを既知として潜在パラメーターをサンプリング」を繰り返す（MH in Gibbs Sampler）か、PMCMCやSMC^2といった逐次モンテカルロベースの手法を用いるなどして求めるが、実装はかなり面倒。
Rなどで手軽に利用できるライブラリはいまのところない。論文用等で作られたコードはいくつか公開されている。
式にすると

$y_t=log(x_t/x_{t-1}) \\y_t =exp(h_t/2)\epsilon_t \\h_{t+1}=\mu+\phi(h_t-\mu)+\sigma_{\eta}\eta_t \\\eta_t,\epsilon_t \sim_{i.i.d}N(0,1)$

...なので時間がない人はStanで求めてしまいましょう。

【インストール】

最初はwindowsでやろうとしたが、Rcppまわりでつまずいてギブアップ。RがProgram Files下にインストールされている（デフォルト設定）と、フォルダ名に空白がはいっているためうまく動かない。
LinuxではチュートリアルどおりにRのコマンドに打ち込めばすぐに使えるようになる。

【実装】

StanはモデルをStanのドメイン特化言語(テキスト形式)で記述する。
data,parameters,modelの三つのブロックで定義する
- data:外部から与えるデータ（データの次元の情報なども整数として与える）
- parameters:推定する潜在変数
- model:変数同士の分布を介した関係
今回はファイルを別にし、sv.stanにモデルを記述し、sv.Rでそれをテキストファイルとして読み込んでいる。
R側からStanを呼ぶには、stan関数を用いる。
- model_code: sv.stanに記述されているテキストデータ
- data: リストでデータを渡す
- iter: イテレーション回数
- chain: チェーンの本数

【結果】
イテレーション10000回（うちバーンイン5000回）でやって30min~1hくらいで終わりました。そんなに多くないですが、パスを見る限り収束は問題なさそうです。

パラメーターの推定値

パス

ボラティリティの推移

2008年のリーマンショックのときにボラティリティが急上昇しているのがわかる。2014年の7月現在はレンジの狭いボックス相場となっており、ボラティリティが低下している。
※中央の青の線が期待値、緑と赤の線がそれぞれ10%,90%の線。

【コード】
下の二つのファイルですべて。事前分布がソフトまかせなのはご愛嬌ということで。

#sv.R

library(rstan)

library(reshape)

library(ggplot2)
model.text=scan("sv.rstan",what = character(),sep="@" , blank.lines.skip = F)

model.text=paste(model.text,collapse="\n")
ind=read.csv("./topix_all.csv")[[1]]

dat=read.csv("./topix_all_return.csv")[[2]]
fitchan <- stan (model_code = model.text,

                 data = list(y=dat,N=length(dat)), 

                 iter  = 10000, 

                 chain = 1)

traceplot (fitchan,par=c("mu","sigeta","phi"))
fit=extract(fitchan,permuted=T)

fit.posterior=get_posterior_mean(fitchan)
h=fit$h

upper=apply(h,2,function(x){quantile(x,.9)})

lower=apply(h,2,function(x){quantile(x,.1)})

expected=apply(h,2,mean)
h_range=melt(data.frame(list(upper=exp(upper/2),lower=exp(lower/2),expected=exp(expected/2),date=ind[2:length(ind)])))

graph=ggplot(h_range,aes(x=as.Date(date),y=value,group=variable,color=variable))+geom_line()
ggdat=melt(fit[c("mu","phi","sigeta")])
graph  <- ggplot (ggdat, aes(x = value))

graph  <- graph + geom_density () + facet_grid(. ~ L1, scales = "free") + theme_bw()

#sv.rstan

data {

        int<lower=1> N;

        real y[N];

    }
parameters {

        real h[N];

        real mu;

        real<lower=-0.999,upper=0.999> phi;

        real<lower=0.0001> sigeta;

}
model {
        h[1] ~ normal(mu,sigeta/sqrt(1-phi*phi));

        y[1] ~ normal(0,exp(h[1]/2));
        for (i in 2:N){

           h[i] ~ normal((1-phi)*mu+phi*h[i-1],sigeta);

           y[i] ~ normal(0,exp(h[i]/2));

        }

}