その１はこちら（極値理論の簡単な解説、Block maximaの実装を行っています）
https://qiita.com/hrkz_szk/items/43debffda9697d9dd7a9
r-block maximaモデル
Block maximaでは最大値のみを抽出しましたが、問題となるのはサンプル数。特にファイナンス系のデータの場合、数十年単位で最大値を確保できないことが多いです。そこで、r-block maximaモデルとThresholdモデルがその対処法になります。r-block maximaは、上からr個の値を抽出します。例えば、2-block maximaモデルでは、最大値とその次に大きい値を抽出します。Thresholdモデルは次回紹介したいと思いますが、これらを統合したのが点過程というモデルです。
データ
ここでは、BMWの株価の損失について、2-block maximaモデルで推定していきたいと思います。具体的に、1年ごとの最大値とその次に大きい値を抽出し、ismevのrlarg.fitを使ってGEV分布への推定を行います。
データに関しては、evirというパッケージにあるdata(bmw)を使用します。こちらは1973年1月2日から1996年7月23日におけるBMW社の株価のログリータンがvectorで入っています。
library(evir)
data("bmw")
head(bmw)
[1] 0.047704097 0.007127223 0.008883307 -0.012440569 -0.003569961 0.000000000
plot(bmw, t="l")
データの抽出
ログリータンをプラスの損失に変換し、さらにパーセント表示に変換。その後、一年間の最大値とその次に大きい値を抽出。
BmwLoss <- -1.0 * bmw * 100
Years <- format(attr(BmwLoss, "time"), "%Y")
attr(BmwLoss, "years") <- Years
Yearu <- unique(Years)
idx <- 1:length(Yearu)
r <- 2
BmwOrder <- t(sapply(idx, function(x)
head(sort(BmwLoss[attr(BmwLoss, "years") == Yearu[x]], decreasing = TRUE), r)))
rownames(BmwOrder) <- Yearu
colnames(BmwOrder) <- paste("r", 1:r, sep = "")
サンプルの分布
plot(Yearu, BmwOrder[, 1], col = "black", ylim = range(BmwOrder),
ylab = "Losses BMW (percentages)", xlab = "",
pch = 20, bg = "black")
points(Yearu, BmwOrder[, 2], col = "blue", pch = 20, bg = "blue")
推定
library(ismev)
BmwOrderFit <- rlarg.fit(BmwOrder)
$conv
[1] 0
$nllh
[1] 77.76214
$mle
[1] 4.4798274 1.7524702 0.3034724
$se
[1] 0.3411369 0.2891531 0.1587266
convのゼロは、うまく推定されたことを示しています。nllhはマイナスのログ尤度の最小値で、mleが尤度関数を用いて推定されたパラメータ、seはパラメータの標準誤差になります。$\xi$が0.3034724とプラスなので、今回の分布はFrechet分布に従うことが分かります。このため、この分布はヘビーテイルであり、損失が有限ではないということになります。
パラメータの検証
信頼区間
私の知る限りismevにはパラメータの信頼区間を自動的に計算してくれるコードを知らないので、ここでは原始的な方法で95%の信頼区間を求めてみます。
lower <- BmwOrderFit$mle - 1.96*BmwOrderFit$se
upper <- BmwOrderFit$mle + 1.96*BmwOrderFit$se
data.frame(lower=lower, upper=upper, row.names = c("mu", "sigma", "xi"))
| lower | upper | |
|---|---|---|
| mu | 3.811199061 | 5.1484558 |
| sigma | 1.185730115 | 2.3192103 |
| xi | -0.007631649 | 0.6145765 |
muとsigmaはゼロをまたいでいませんが、最も大切なxiがゼロの可能性が。しかし、ゼロをまたいでいると行ってもとても小さな幅なので、ここでは深追いはせずに…。
モデルの検証
rlarg.diag(BmwOrderFit)
Block maximaモデルの事例と同じように、Probability PlotとQuantilte Plotとともに右上の方がラインに沿っておらず、推定したモデルでは説明しきれていないことを示しています。しかし、再現レベルはサンプルデータが95%の内側に入っているので、この点は良いのかなと思います。
続編
その３はこちら。Thresholdモデルを使った分析を行っています。
https://qiita.com/hrkz_szk/items/8b01552e0e09583d4c35
参考文献
