TOPIXのボラティリティをStochastic Volatilityモデル + R + Stanで推定する
はじめに
この記事はマケデコ Advent Calendar 2023の18日目の記事です。1日遅れですが、枠が空いていたので飛び入り参加してみました。
状態空間モデルによるボラティリティモデルのStochastic Volatilityの論文をR + Stanで実装してみました。Stanの推定結果のプロットにはbayesplotとggplot2を用いています。
- 非線形な状態空間モデルによるボラティリティモデルであるStochastic Volatility (SV) モデルをStanで実装することで、TOPIXのボラティリティを推定しました。
- 推定されたボラティリティは、2008年のリーマンショックと、2011年の東日本大震災、2020年の新型コロナウイルスによる市場の急落局面で非常に高まっていることを確認できました。
- ボラティリティが一度高まるとしばらくボラティリティが高い日が続く現象であるボラティリティ・クラスタリングも確認できました。
ボラティリティの定式化
ボラティリティとは、金融商品の値動きの変動の大きさを指すパラメータです。どのくらいの損失がどのくらいの確率で発生するかを示すVaR（Value at Risk）に使われるなど、金融工学の根幹をなす値です。
以下は株式に限らず、例えば為替など、金融商品であれば何でもよいですが、株式とします。
$S_t$を$t (1, \dots,\ T)$日における株式の価格とするとき、$t$日における対前日の収益率$r_t$は$r_t=\log S_t - \log S_{t-1}$となります。これを対数収益率といいます。
このとき、ボラティリティとは以下の$\sigma_t$、あるいは$\sigma_t^2$を指します1。
$$ \begin{aligned} r_t &= E_{t-1}[r_t] + \epsilon_t \\\ \epsilon_t &= \sigma_t z_t, \quad \sigma_t > 0, \quad z_t \sim i.i.d., \quad E[z_t] = 0, \quad Var[z_t] = 1 \end{aligned} $$
$\sigma_t$と$\sigma_t^2$のどちらをボラティリティと呼ぶかは文献によりますが、以降$\sigma_t$をボラティリティと呼びます。
この$\sigma_t$を推定する方法は大きく分けて三通りあります。
一つ目の方法は、$\sigma_t$を過去一定期間の$r_t$の標準偏差とする方法です。過去1年～3年＝250～750営業日ということで$j=250$や$j=750$などとして、$\{r_{t-j+1},r_{t-j+2},\dots r_{t}\}$の標準偏差とします。もちろん直近のボラティリティを重視したいなら$j$をより小さい値に設定します。この$t$を1日ずつずらしてローリング計算することでボラティリティの時系列を得ることができます。
この方法はシンプルですが、標準偏差を計算している$t-j+1,\dots,t$の間は$\sigma_t$が一定であることを暗に仮定しています。実際はそうではありませんので、精緻に求めるなら残りの二つの方法のどちらかを用いることになります。
二つ目の方法は、$\sigma_t$を統計的なモデルで定式化するものです。
このタイプのモデルは有名なGARCHモデルをはじめ色々ありますが、本記事ではStochastic Volatility (SV) モデルと呼ばれるモデルのうち、シンプルな以下のモデルを推定することにします。
$$ \begin{aligned} y_t &= \exp(x_t/2) \epsilon_t, \quad \epsilon_t \sim N(0,1) \\\ x_{t+1} &= \mu + \phi(x_t - \mu) + \eta_t, \quad \eta_t \sim N(0,\sigma_{\eta}^2) \\\ x_1 & \sim N(0,\sigma_{\eta}^2/(1-\phi^2)) \end{aligned} $$
ここで$y_t$は対数収益率$r_t$であり、$\exp(x_t/2)$がボラティリティ$\sigma_t$です。$x_t$は定常な過程と仮定し、$|\phi|<1$です。
このSVモデルは、最初の式で$E_{t-1}[r_t]=0$とした上で、$x_t = \log \sigma_t^2$がAR(1)モデルに従うことを意味します。
このモデルは1本目の式が観測方程式、2本目の式が状態方程式の非線形な状態空間モデルですので、パーティクルフィルタのような非線形な状態空間モデルでも推定できるタイプのフィルタ系の手法か、Stanなどを用いてMCMCで推定することになります。
最後に三つ目の方法としては、以上二つのように収益率$r_t$の時系列から$\sigma_t$を推定するのではなく、分単位のような細かい収益率データを用いて直接$\sigma_t$を求めるアプローチがあります。本記事からは外れるので詳細は触れませんが、$t$日における1分間隔や5分間隔程度の細かい間隔の収益率の2乗を1日分足し合わせたものが$\sigma_t$の精度のよい推定量（Realized Volatility2と呼びます）になることが知られています。
環境
R + RStanです。Stanコードを書いてrstanというRのパッケージから呼び出します。
Rを用いているのはベイズモデリングは文献もパッケージもRが豊富なためです。あとこの記事のコードは自分で昔書いたコードを流用しているのですが、それがRだったからというのもあります。
- Windows 10
- R 4.3.1
- httr 1.4.7
- rstan 2.32.3
- bayesplot 1.10.0
- patchwork 1.1.3
library(tidyverse)
library(httr)
library(rstan)
library(bayesplot)
library(patchwork)
データの取得
今回は上で挙げたSVモデルをRとStan (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するだけのシンプルな作りです。
# J-Quantsに登録したメールアドレスとパスワード
mail_address <- "MAIL_ADDRESS"
password <- "PASSWORD"
resp <- httr::POST(
"https://api.jquants.com/v1/token/auth_user",
body=jsonlite::toJSON(
list(mailaddress=mail_address, password=password),
auto_unbox=TRUE
)
)
refresh_token <- httr::content(resp)$refreshToken
resp <- httr::POST(
"https://api.jquants.com/v1/token/auth_refresh",
query=list(refreshtoken=refresh_token)
)
id_token <- httr::content(resp)$idToken
resp <- httr::GET(
"https://api.jquants.com/v1/indices/topix",
query=list(from="2008-05-07", to="2023-12-18"),
httr::add_headers(Authorization=glue::glue("Bearer {id_token}"))
)
topix <- httr::content(resp)$topix |>
dplyr::bind_rows() |>
tibble::as_tibble()
最初のPOSTでは{"mailaddress": "mail@example.com", "password": "hoge"}のようなJSONをbodyに入れます。RでこのJSONを作るには、jsonlite::toJSON(list(mailaddress=mail_address, password=password), auto_unbox=TRUE)とします。
jsonlite::toJSON()の引数auto_unboxはデフォルトではFALSEなのですが、jsonlite::toJSON(list(mailaddress=mail_address, password=password), auto_unbox=FALSE)は{"mailaddress": ["mail@example.com"], "password": ["hoge"]}となってしまいます。Rではmail@example.comのような文字列は長さ1のcharacter型のベクトルであるため、そのまま長さ1のリストができてしまうからです。これを防ぐためにauto_unbox=TRUEを指定します。
TOPIXの終値から対前日の対数収益率を計算します。100倍して%表記にします。また、最初の1日の収益率はNAになるので最初の1日を除いておきます。
df <- topix |>
mutate(Date=as.Date(Date, "%Y-%m-%d")) |>
mutate(ret=(log(Close) - log(lag(Close, 1)))*100) |>
slice(-1)
こんな感じのデータです。
df
#> # A tibble: 3,823 × 6
#> Date Open High Low Close ret
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2008-05-08 1384. 1386. 1373. 1373. -1.47
#> 2 2008-05-09 1372. 1374. 1341. 1342. -2.30
#> 3 2008-05-12 1331. 1345. 1327. 1343. 0.0767
#> 4 2008-05-13 1351. 1364. 1344. 1360. 1.28
#> 5 2008-05-14 1360. 1376. 1351. 1373. 0.951
#> 6 2008-05-15 1382. 1404. 1382. 1393. 1.43
#> 7 2008-05-16 1405. 1412. 1391. 1396. 0.215
#> 8 2008-05-19 1400. 1410. 1397. 1404. 0.599
#> 9 2008-05-20 1402. 1410. 1394. 1400. -0.315
#> 10 2008-05-21 1385. 1386. 1361. 1370. -2.15
#> # ℹ 3,813 more rows
Stanで実装
上に挙げたSVモデルをStanコードで書きます。
data {
int N;
vector[N] y;
}
parameters {
vector[N] x;
real mu;
real<lower=-1,upper=1> phi;
real<lower=0> sigma_eta;
}
transformed parameters {
real phi_beta;
phi_beta = (phi + 1)/2;
real sigma_eta_square;
sigma_eta_square = sigma_eta^2;
}
model {
mu ~ normal(0, 1);
phi_beta ~ beta(20, 1.5);
sigma_eta_square ~ inv_gamma(5.0/2, 0.05/2);
x[1] ~ normal(0, sigma_eta/sqrt(1 - phi^2));
x[2:N] ~ normal(mu + phi * (x[1:(N-1)] - mu), sigma_eta);
y ~ normal(0, exp(x/2));
}
generated quantities {
vector[N] vol;
vol = exp(x/2);
vector[N] y_pred;
for (i in 1:N) {
y_pred[i] = normal_rng(0, exp(x[i]/2));
}
}
以下、実装のポイントです。
$\mu, \phi, \sigma_{\eta}$の事前分布は、SVモデルの元の論文であるKim, Shephard and Chib (1998) 3や、それを日本株に適用した大森, 渡部 (2007)4 にある以下の分布を用いました。IGは逆ガンマ分布です。無情報事前分布ではなく、事前分布を書いてあげると収束しやすくなります。
$$ \begin{aligned} \mu & \sim N(0,1) \\\ \frac{\phi+1}{2} & \sim Beta(20,1.5) \\\ \sigma_{\eta}^2 & \sim IG(5/2, 0.05/2) \end{aligned} $$
generated quantitiesのvolが今回求めたいボラティリティ$\sigma_t$です。MCMCの結果から$x$の事後分布を取り出してR上で$\exp(x/2)$を計算することで求めることもできますが、Stanで求めておきます。
generated quantitiesのy_predは、各パラメータの事後分布から乱数で生成した$y_t$の値です。これはパラメータの推定自体には不要ですが、あとでモデルの事後診断に使用するために求めています。
あとは高速化のためにできるだけベクトル化をします。
このStanコードを”svmodel.stan”で保存して以下のRコードでキックすることでMCMCの推定を行います。
i9-9900Kで実行しました。chains=4, iter=110000, warmup=10000, thin=50で8時間程度かかりました。iterをかなり大きくしてthinで間引かないと$\phi$などのパラメータに自己相関が残っていました。
# Stanのおまじない（上は並列化、下はstanコードが変わらない限り再コンパイルしない）
options(mc.cores=parallel::detectCores())
rstan_options(auto_write=TRUE)
# C++コードへのコンパイル
mod <- rstan::stan_model("svmodel.stan")
# MCMCサンプリング
fit <- rstan::sampling(
mod,
data=list(N=nrow(df), y=df$ret),
chains=4, iter=110000, warmup=10000, thin=50, seed=1234
)
MCMCのチェック
MCMCがうまくいっているかチェックします。
MCMCが終わったら早速パラメータの値を見たいところですが、収束していないということがよくあるのでチェックしましょう5。収束していなければパラメータが正しく推定されていないということです。
詳細に知りたい方は、例えばこちらの分かりやすい記事をご覧ください。
[R] [stan] bayesplot を使ったモンテカルロ法の実践ガイド - ill-identified diary
このあたりの事後診断の可視化はbayesplotが超便利です。ggplot2ベースなので出力のグラフにggplot2の関数で軸などのレイアウトを微調整できるのもナイスです。
まずはRhatです。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の自己相関も消える気がします。
bayesplot::mcmc_acf_bar(fit, pars=c("mu", "phi", "sigma_eta"))
ちなみにbayesplot::mcmc_acf_bar()をはじめ、bayesplotの描画関数はparsでパラメータを指定し忘れると全パラメータ描画します。巨大なモデルだとこれでRStudioがクラッシュすることがありますので気を付けましょう。特に時系列の状態空間モデルはパラメータの数が多いのでクラッシュしがちです（今回はxがx[1]からx[3823]まである）。モデルをRDSファイルに保存し忘れた状態だと推定結果が消えて悲しいことになります。
そしてトレースプロットを見ます。こちらもパラメータを絞っています。
これはchainごとのサンプルの推移で、線が混ざり合っていると初期値によらず同じ値に収束している＝局所解に落ちていないことを示します。Rhatが大きいときはトレースプロットが混ざり合っていないので、Rhatと合わせてチェックします。いい感じですね！
bayesplot::mcmc_trace(fit, pars=c("mu", "phi", "sigma_eta"))
最後に、モデルが現実のデータをよく説明するように定式化されているなら、パラメータの事後分布から乱数を振って得られる目的変数（ここでは$y_t$）は、実際の$y_t$と同じような分布で得られるはずです。これを図示してみます。
濃い線は実際のy_t（対数収益率）の分布、薄い線はモデルから生成した8000系列の対数収益率y_predのうち最初の10系列の分布です。薄い線はそれぞれの系列で10本あります。8000系列は描けないので10本に絞っています。
濃い線と薄い線が大体重なっているので悪くなさそうです。
bayesplot::ppc_dens_overlay(df$ret, rstan::extract(fit)$y_pred[1:10,])
自分でモデルを一から組む場合は、モデルのパラメータはMCMCで正しく推定されていても、そもそもモデル自体が元のデータを全然説明できていない誤ったモデルであることがあります（SVモデルは幅広く用いられているモデルなので、今回ちゃんと当てはまるのはある意味当然なわけですが）。モデルが元のデータに当てはまっているかどうかを示してくれるのが上のプロットです。
bayesplot::ppc_*系の関数はこのような事後予測の確認に役立つものが色々あって便利です。
以上がMCMCの結果の基本的なチェック方法です。
パラメータの推定結果
以下がパラメータの推定結果です。ただしprintが長くなるので一部のパラメータに絞っています。
print(fit, pars=c("mu", "phi", "sigma_eta"), digits_summary=3)
#> Inference for Stan model: anon_model.
#> 4 chains, each with iter=110000; warmup=10000; thin=50;
#> post-warmup draws per chain=2000, total post-warmup draws=8000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> mu 0.127 0.001 0.106 -0.079 0.058 0.127 0.197 0.336 7883 1
#> phi 0.966 0.000 0.007 0.951 0.961 0.966 0.971 0.978 4028 1
#> sigma_eta 0.213 0.000 0.020 0.176 0.198 0.212 0.226 0.255 3167 1
#>
#> Samples were drawn using NUTS(diag_e) at Wed Dec 20 00:41:18 2023.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
推定したボラティリティを見てみましょう。rstan::extract()を用いてパラメータの事後分布のサンプルからvolの時系列を作ります。
点推定値として事後中央値を用います。あわせて95%ベイズ信用区間も示したいので、中央値と2.5%タイル点と97.5%タイル点を取り出します。volの事後分布の(110000 (iter) - 10000 (warmup))/50 (thin) * 4 (chain) = 8000個のサンプルを小さい順に並び変えて50%と2.5%と97.5%のタイル点を取り出すことで得られます。
# 8000 (サンプル) x 3823 (vol[1] - vol[3823])のmatrix
mat <- rstan::extract(fit, "vol")$vol
vol_stat <- tibble::tibble(
vol_median=apply(mat, 2, \(x) quantile(x, 0.5)),
vol_lower=apply(mat, 2, \(x) quantile(x, 0.025)),
vol_upper=apply(mat, 2, \(x) quantile(x, 0.975))
) |>
mutate(Date=df$Date) |>
relocate(Date)
res <- left_join(df, vol_stat, by="Date")
res
#> # A tibble: 3,823 × 9
#> Date Open High Low Close ret vol_median vol_lower vol_upper
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2008-05-08 1384. 1386. 1373. 1373. -1.47 1.36 0.892 2.20
#> 2 2008-05-09 1372. 1374. 1341. 1342. -2.30 1.37 0.908 2.15
#> 3 2008-05-12 1331. 1345. 1327. 1343. 0.0767 1.34 0.894 2.10
#> 4 2008-05-13 1351. 1364. 1344. 1360. 1.28 1.33 0.896 2.10
#> 5 2008-05-14 1360. 1376. 1351. 1373. 0.951 1.33 0.894 2.06
#> 6 2008-05-15 1382. 1404. 1382. 1393. 1.43 1.33 0.892 2.05
#> 7 2008-05-16 1405. 1412. 1391. 1396. 0.215 1.32 0.887 2.07
#> 8 2008-05-19 1400. 1410. 1397. 1404. 0.599 1.33 0.903 2.05
#> 9 2008-05-20 1402. 1410. 1394. 1400. -0.315 1.36 0.918 2.06
#> 10 2008-05-21 1385. 1386. 1361. 1370. -2.15 1.39 0.963 2.09
#> # ℹ 3,813 more rows
以下のグラフの上はボラティリティ$\sigma_t$（赤線は$\sigma_{t}$の事後分布の中央値、青いバンドは95%ベイズ信用区間）、下はTOPIXの終値です。$\sigma_t = a$であれば、TOPIXの収益率の標準偏差がa[%]であることを意味します。
code
p_vol <- res |>
ggplot(aes(Date))+
theme_light()+
geom_ribbon(aes(ymin=vol_lower, ymax=vol_upper), fill="lightsteelblue1", alpha=0.5)+
geom_line(aes(y=vol_upper), color="lightsteelblue1")+
geom_line(aes(y=vol_lower), color="lightsteelblue1")+
geom_line(aes(y=vol_median), color="firebrick")+
scale_x_date(breaks=scales::date_breaks("1 year"), date_labels="%y")+
labs(
x="date (year)",
y="volatility (sigma_t)",
subtitle="red: estimated (median), light blue: 95% CI"
)
p_topix <- res |>
ggplot(aes(x=Date, y=Close))+
theme_light()+
geom_line()+
scale_x_date(breaks=scales::date_breaks("1 year"), date_labels="%y")+
labs(x="date (year)", y="TOPIX close")
patchwork::wrap_plots(p_vol, p_topix, ncol=1)
2008年のリーマンショック、2011年の東日本大震災、2020年の新型コロナウイルスの市場急落局面でボラティリティが高まっていることが分かります。
なお、グラフはpatchworkで並べました。複数のggplotオブジェクトを綺麗に並べられて重宝するパッケージです。
最後に$\phi$のパラメータの事後分布を見てみます。先程示した通り、$\phi$ = 0.966 (95%CI: 0.951-0.978)でした。
bayesplot::mcmc_hist(fit, pars="phi")
このパラメータは過去のボラティリティがどの程度後を引くかを示すパラメータであり、1に近いということは持続性がかなりあるということを示します。ボラティリティが一度上昇するとしばらくボラティリティが高い日が続くということであり、この現象をボラティリティ・クラスタリングといいます。
SVモデルを推定した論文をサーベイすると$\phi$の推定値は0.8から0.995までの値となっているという論文がありますが6、この先行研究と整合的です。
おわりに
Stanを使うと柔軟にモデルを組んで解釈ができて楽しいですね。
SVモデルの一番シンプルなものを推定してみましたが、SVモデルには収益率の裾の厚さを表現するために誤差項をt分布としたり、収益率がマイナスの日の方がボラティリティが高まる（ボラティリティの非対称性といいます）ことを表現したりするモデルなど、色々な発展形があります。これらをStanで組んでみてもよいでしょう。
最後に株価データで時系列モデリングをやってみたいという方におすすめの書籍をご紹介します。
- ボラティリティ変動モデル (シリーズ
現代金融工学)
- 統計モデル的アプローチによるボラティリティモデルです。ボラティリティモデルをしっかり学ぶならこれです。MCMCによる推定も少し触れられています。コード例はなく数式展開の理論の本です。
- 経済・ファイナンスのための カルマンフィルター入門
(統計ライブラリー)
- 前半はカルマンフィルタの解説、後半は経済・金融データへのカルマンフィルタの適用事例の紹介です。前半はカルマンフィルタの導出を丁寧に書いていて理解しやすかったです。後半はカルマンフィルタでベータ値を推定したりペアトレーディングの銘柄を発掘したりと面白い事例が豊富です。コードはないですが読みやすいです。
関連記事
式の通り、$z_t$は正規分布とは指定していないので、$r_t$も正規分布以外の分布を仮定することができます。実証的には$r_t$は裾が厚く正規分布ではないとされます。 ↩︎
例えばこの辺りの論文をご覧ください。渡部敏明, 佐々木浩二 (2006), 「ARCH型モデルと”Realized Volatility”によるボラティリティ予測とバリュー・アット・リスク」, 金融研究, 25 別冊(2), 39-74. ↩︎
Kim, S., N. Shephard, and S. Chib (1998), “Stochastic Volatility: Likelihood Inference and Comparison with ARCH Models”, Review of Economic Studies, 65, 361-393. ↩︎
大森裕浩, 渡部敏明 (2007), 「MCMC法とその確率的ボラティリティモデルへの応用」CIRJEディスカッションペーパー, J-173, 1-39. ↩︎
何日もかけて回したMCMCの結果をチェックしたら全然収束していない悲しみもあるある ↩︎
Jacquier, E., N. Polson, and P. Rossi (2004), “Bayesian Analysis of Stochastic Volatility Models (with Discussion)”, Journal of Business & Economic Statistics, 12, 371-417. ↩︎