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

R 金融工学

引き続き日経VIネタ。
日経VIの挙動は、株価のような対数幾何分布に従うような確率過程（そうは見えないけどそうだと思っておく！）ではなく、平均回帰的であり、上昇するも下落するも急なものであるというのが(私的に)一般的に知られた事実である。

ここでは日経VIがどのくらいの速さで平均回帰するのか？、すなわち急上昇したのちにどれくらいの日数をかけて落ちるのかについて考えたい。
物理学では元の値の半分の大きさ（すなわち、金融でいうところの半値に相当）になるまでにかかる時間(日数)を半減期（とか、特徴的な長さスケール）だなんて呼ぶのでそれを計算してやる。
これがわかるとおおよそ日経VIがピークになってからその半分の値になるまでの時間（日数）がわかるのである。

データの取得

データは従前の記事同様、Quandlパッケージを活用して取得する。

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")

ピーク位置の同定

これが今回の一番悩んだところ。
ここでは以下のアイディアに基づいて抽出しているが、これが一般的なのかどうかはさっぱりわからない*1。
後から調べた感じ、

この辺を参考にパッケージ使えばよかった…。


ここではまず、時系列データの各点を中心として半径(単系列での半径なので、ただの一次元円の半径、すなわちただの距離)k/2以内にあるデータに対してfunという処理をかませるrollapp()関数を定義する。
これは、時系列の一番後ろ（時間的に昔）からデータを舐めて処理するzooパッケージにあるrollapplyとは違い、各時点の未来・過去k/2単位(今回の場合日次データなので、単位は日)にまたがったデータに対して処理を適用する関数。

#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    if((i-w) <= 0 || (i+w) > length(x)){
      NA
    } else{
      fun(x[(i-w):(i+w)])
    }
  })
}

こいつを使って

  • 上値ピークポイント：直近150日(未来・過去方向それぞれ75日)分のデータの中で最も大きく、かつ、日経VIの値自体は30以上
  • 下値ピークポイント：直近125日(未来・過去方向それぞれ62日)分のデータの中で最も小さい

点を同定し、それぞれmax_point, min_pointとする。
NAはよしなに削除で。

x <- vi$Close
max_point <- rollapp(x, 150, function(x){x_middle <- x[1+length(x)/2]; (x_middle==max(x)) && (x_middle > 30)})
max_point <- replace(max_point, is.na(max_point), FALSE)
min_point <- rollapp(x, 125, function(x){x_middle <- x[1+length(x)/2]; (x_middle==min(x))})
min_point <- replace(min_point, is.na(min_point), FALSE)

ここで計算したピークポイントがどの辺にあるのかというと、下図のとおり、まぁ見た感じよさそうだぞと。
（上(下)値ピークが赤(青)点、日経VIの値は黒線で表示）
（逆に、こういう見た目になるように↑の日数とか日経VIの水準を足きりしている点に注意）

plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)
points(seq_len(length(x))[min_point], x[min_point], col = "blue", pch=16)


上昇→下落フェーズのデータの抽出

ここで、「定義した上値・下値の間にあるデータ」のみを抽出する。
つまり上の図でいうところの、各赤点・青点の間にある日経VIのデータを抽出するということだ。
そのためには汚いが、データを全舐めして、下記のようなコードでデータを抽出する。
index変数が0以外の値で、かつ、共通の値を持つ箇所がおなじ上昇→下落局面に対応している。

index <- numeric(length(x))
counter <- 0
mode <- "max"
for(i in seq_len(length(x))){
  if(mode == "max" && max_point[i]){
    mode <- "min"
    counter <- counter + 1
    index[i] <- counter
  } else if(mode == "min" && min_point[i]){
    mode <- "max"
  } else if(mode == "min"){
    index[i] <- counter
  }
}
xs <- lapply(seq_len(max(index)), function(i){x[which(index == i)]})

結果は、各赤点・青点の間にあるデータごとのリストになっている。

半減期の同定

上で作ったリストの各要素が上昇→下落局面に対応する
NIKKEI-VI = a*\exp\{-t/\tau\}

という回帰式に対してτを推定することがメインであり、このτが半減期、すなわち日経VIがどの程度の速度で下落するのかの指標になっている*2。
推定する部分を関数にし、元データ・τ・推定結果(データ自体をlog化し、そのデータの誤差が正規分布になってほしいなと仮定して、線形回帰に直して計算させる)を返却させるようにしよう。

# x = a*exp(-t/tau)
# log(x) = log(a) - 1/tau * t
decay <- function(x) {
  df <- data.frame(t=seq_len(length(x)), logx=log(x))
  regression <- lm(logx~t, data=df)
  tau <- -1/coef(regression)["t"]
  list(x=x, tau=tau, regression=regression)
}

これを、上で作った全データに適用してやると

result <- lapply(xs, decay)

となる。
得られたτは

> tau <- sapply(result, function(x){x$tau})
> tau
        t         t         t         t         t         t         t         t         t 
 56.60291  58.95572 165.28740 395.43998 201.99652 109.37092  40.24398  33.13756  27.06593

となり、30くらいから400くらいまでとかなり幅広だ。
各々の結果を出したデータと回帰の当てはまり具合を下記のように目視で確認すると、長い半減期の系はだいたいフィッティングがおかしいので、無視してよさそうだ。
これは、どちらかというと、データの中に小さなピークが立っていることからもわかるように、ピークの同定がうまくないからなんだろうなと。

従って、短い方を意識して大体、半減期は10-30日程度と意識して日経VI先物に投資すればよいだろう(と期待したい)。
τ=56.60291

τ=58.95572

τ=165.28740

τ=395.43998

τ=201.99652

τ=109.37092

τ=40.24398

τ=33.13756

τ=27.06593

*1:まじめに調べてもいない

*2:より正確にはこの値の0.7倍程度か https://ja.wikipedia.org/wiki/%E5%8D%8A%E6%B8%9B%E6%9C%9F