金融市場が混乱する間隔を日経平均ボラティリティー・インデックス（日経平均VI）から推定したい

日経ボラティリティー・インデックス（日経平均VI）祭りはまだ続く。
やりたいことは簡単で

日経VIが跳ねるタイミング＝金融市場が混乱

だと考えて、"日経VIが跳ねるタイミング"の確率分布がどのような形になるのかを推定するというお話。
個人的には指数分布なんじゃね？と思っていたが、違った。

データの準備

ほぼ

日経平均ボラティリティー・インデックス（日経平均VI）の半減期が知りたい - My Life as a Mock Quant

でやった内容とおなじで、日経VIが跳ねたポイント（下図赤点）を手でやるのがめんどうなので、Rに同定させるためのコード。
↑のLINKとはrollapp()関数がちょっと違ってるので注意

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")
#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    index_s <- max(i-w, 1)
    index_e <- min(i+w, length(x))
    middle <- index_s + i-w
    fun(1+i-index_s, x[index_s:index_e])
  })
}
x <- vi$Close
max_point <- rollapp(x, 74, function(i, x){x[i]==max(x)})
max_point <- replace(max_point, is.na(max_point), FALSE)
plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)

このデータの赤点同士の時間間隔がどういう確率分布に従っているのかを評価したい。

推定(指数分布)

指数分布で確率密度関数を推定してみる

# 適当にdata.frame化
df <- dplyr::group_by(data.frame(index=cumsum(max_point)), index) %>% summarize(n=n())
lambda <- 1/mean(df$n)
# 下記でも同じ
# library("MASS")
# exponential <- coef(fitdistr(df$n, "exponential"))
hist(df$n, freq=FALSE)
curve(dexp(x, rate=lambda),add=TRUE,lty=2)

う、うーん、これではない感。

推定(ワイブル分布)

ワイブル分布 - Wikipedia

によると、ワイブル分布は「時間に対する劣化現象や寿命を統計的に記述するためにも利用される」とあるので、こいつでやってみたい。
金融市場が破裂する＝寿命が尽きると考えるわけだ。

library("MASS")
weibull <- coef(fitdistr(df$n, "weibull"))
hist(df$n, freq=FALSE)
curve(dweibull(x, weibull["shape"], weibull["scale"]), add=TRUE, lty=2)