# 8章
# 8.4.1項
library(evd)
portpirie
## [1] 4.03 3.83 3.65 3.88 4.01 4.08 4.18 3.80 4.36 3.96 3.98 4.69 3.85 3.96 3.85
## [16] 3.93 3.75 3.63 3.57 4.25 3.97 4.05 4.24 4.22 3.73 4.37 4.06 3.71 3.96 4.06
## [31] 4.55 3.79 3.89 4.11 3.85 3.86 3.86 4.21 4.01 4.11 4.24 3.96 4.21 3.74 3.85
## [46] 3.88 3.66 4.11 3.71 4.18 3.90 3.78 3.91 3.72 4.00 3.66 3.62 4.33 4.55 3.75
## [61] 4.08 3.90 3.88 3.94 4.33
plot(portpirie, type="n") # プロットの枠だけ
length(portpirie)
## [1] 65
text(1:65, portpirie, paste(1:65)) # 番号(時刻)でプロットする
abline(h=4.2)

clusters(portpirie, u=4.2, r=3) # uとrのみ設定
## $cluster1
## 9 10 11 12
## 4.36 3.96 3.98 4.69
##
## $cluster2
## 20 21 22 23 24 25 26
## 4.25 3.97 4.05 4.24 4.22 3.73 4.37
##
## $cluster3
## 31
## 4.55
##
## $cluster4
## 38 39 40 41 42 43
## 4.21 4.01 4.11 4.24 3.96 4.21
##
## $cluster5
## 58 59
## 4.33 4.55
##
## $cluster6
## 65
## 4.33
##
## attr(,"acs")
## [1] 2.166667
clusters(portpirie, u=4.2, r=3, ulow=4, rlow=2) # ulowとrlowも設定
## $cluster1
## 9
## 4.36
##
## $cluster2
## 12
## 4.69
##
## $cluster3
## 20 21 22 23 24 25 26
## 4.25 3.97 4.05 4.24 4.22 3.73 4.37
##
## $cluster4
## 31
## 4.55
##
## $cluster5
## 38 39 40 41 42 43
## 4.21 4.01 4.11 4.24 3.96 4.21
##
## $cluster6
## 58 59
## 4.33 4.55
##
## $cluster7
## 65
## 4.33
##
## attr(,"acs")
## [1] 1.857143
clusters(portpirie, u=4.2, r=3, cmax = TRUE) # 群れの最大値を取る
## 12 26 31 41 59 65
## 4.69 4.37 4.55 4.24 4.55 4.33
## attr(,"acs")
## [1] 2.166667
clusters(portpirie, u=4.2, r=3, ulow=3.8, plot = TRUE) # 群れのプロットをする

tvu <- c(rep(4.2, 20), rep(4.1, 25), rep(4.2, 20))
clusters(portpirie, tvu, 3, plot = TRUE) # 閾値を時刻によって変える

exi(portpirie, u=4.2, r = 3, ulow = 3.8)
## [1] 0.5384615
tvu <- c(rep(4.2, 20), rep(4.1, 25), rep(4.2, 20))
exi(portpirie, u=tvu, r = 1)
## [1] 0.8
exi(portpirie, u=tvu, r = 0)
## [1] 0.9402985
exiplot(portpirie, tlim=c(2,4.3))
exiplot(portpirie, tlim=c(2,4.3), r=3, add=T)

# 8.4.2項
fpot(portpirie, 4.2) # 前章と同様
##
## Call: fpot(x = portpirie, threshold = 4.2)
## Deviance: -23.3281
##
## Threshold: 4.2
## Number Above: 13
## Proportion Above: 0.2
##
## Estimates
## scale shape
## 0.15380 -0.02513
##
## Standard Errors
## scale shape
## 0.08204 0.45577
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 23
## Gradient Evaluations: 5
fpot(portpirie, 4.2, r=3) # cmax=Fなので上と同じ
##
## Call: fpot(x = portpirie, threshold = 4.2, r = 3)
## Deviance: -23.3281
##
## Threshold: 4.2
## Number Above: 13
## Proportion Above: 0.2
##
## Estimates
## scale shape
## 0.15380 -0.02513
##
## Standard Errors
## scale shape
## 0.08204 0.45577
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 23
## Gradient Evaluations: 5
# fpot(portpirie, 4.2, r=3, cmax=T) # 計算できない
fpot(portpirie, 4.2, r=3, cmax=T, method="Nelder-Mead")
##
## Call: fpot(x = portpirie, threshold = 4.2, cmax = T, r = 3, method = "Nelder-Mead")
## Deviance: -28.88447
##
## Threshold: 4.2
## Number Above: 13
## Proportion Above: 0.2
##
## Clustering Interval: 3
## Number of Clusters: 6
## Extremal Index: 0.4615
##
## Estimates
## scale shape
## 0.7091 -1.4472
##
## Standard Errors
## scale shape
## 0.0000020 0.0006548
##
## Optimization Information
## Convergence: 10
## Function Evaluations: 293
# 群れの最大値のみで推定
# 8.4.3項
set.seed(100)
marma(5, p=1, q=2, psi=0.2, theta=c(0.3,0.2))
## [1] 9.5561305 2.8668391 1.9112261 0.8514774 10.7391557
mar(5, psi=0.2, init=1)
## [1] 4.000113 5.145979 1.904391 2.958200 0.591640
mma(5, theta=0.5, rand.gen=rexp)
## [1] 0.5655238 0.1902905 0.4216077 0.2108038 0.4979445
# 8.5節
library(evd)
set.seed(100)
temp <- marma(500, p=1, q=2, psi=0.2, theta=c(0.3,0.2)) # MARMA時系列を発生
plot(temp) # 時系列の実現値をプロット

plot(temp, type="l") # 表示方法を変える

plot(temp, ylim=c(0,100)) # 範囲を変える

plot(temp, ylim=c(0,50)) # 範囲を変える

mrlplot(temp) # 標本平均超過プロット

tcplot(temp, c(10,30)) # 推定プロット


sum(temp>13)
## [1] 42
clusters(temp, u=13, r=3) # 群れを調べる
## $cluster1
## 15 16 17
## 13.962130 4.188639 13.038252
##
## $cluster2
## 49
## 38.16393
##
## $cluster3
## 53
## 28.6373
##
## $cluster4
## 61 62 63
## 124.95730 37.48719 24.99146
##
## $cluster5
## 80 81 82 83 84
## 73.60953 22.08286 98.30403 29.49121 19.66081
##
## $cluster6
## 158
## 16.74253
##
## $cluster7
## 194 195 196 197 198 199 200
## 46.642619 21.243652 9.328524 31.425522 186.014291 55.804287 37.202858
## 201
## 13.035300
##
## $cluster8
## 236
## 17.86497
##
## $cluster9
## 241 242 243 244 245
## 43.419027 13.025708 8.683805 44.494197 13.348259
##
## $cluster10
## 259 260 261
## 15.723594 4.717078 39.181266
##
## $cluster11
## 265
## 22.35944
##
## $cluster12
## 280
## 26.35825
##
## $cluster13
## 294
## 13.81951
##
## $cluster14
## 316 317 318
## 65.09524 19.52857 13.01905
##
## $cluster15
## 337
## 34.20305
##
## $cluster16
## 411
## 26.26472
##
## $cluster17
## 419
## 16.62869
##
## $cluster18
## 427
## 17.12813
##
## $cluster19
## 463
## 24.97662
##
## $cluster20
## 472 473 474 475 476
## 32.451590 9.735477 65.551418 19.665425 13.110284
##
## attr(,"acs")
## [1] 2.1
exi(temp, u=13, r=3) # 極値収縮率の推定
## [1] 0.4761905
fpot(temp, 13) # 超過データをすべて使う
##
## Call: fpot(x = temp, threshold = 13)
## Deviance: 342.6551
##
## Threshold: 13
## Number Above: 42
## Proportion Above: 0.084
##
## Estimates
## scale shape
## 13.8281 0.4522
##
## Standard Errors
## scale shape
## 4.0802 0.2611
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 34
## Gradient Evaluations: 17
fpot(temp, 13, r=3, cmax=T) # 群れの最大のみ使う
##
## Call: fpot(x = temp, threshold = 13, cmax = T, r = 3)
## Deviance: 177.5535
##
## Threshold: 13
## Number Above: 42
## Proportion Above: 0.084
##
## Clustering Interval: 3
## Number of Clusters: 20
## Extremal Index: 0.4762
##
## Estimates
## scale shape
## 20.3315 0.4267
##
## Standard Errors
## scale shape
## 8.3108 0.3552
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 25
## Gradient Evaluations: 15
fpot(clusters(temp, u=13, r=3, cmax=T),13) # 群れの最大ベクトルをすべて使う
##
## Call: fpot(x = clusters(temp, u = 13, r = 3, cmax = T), threshold = 13)
## Deviance: 177.5535
##
## Threshold: 13
## Number Above: 20
## Proportion Above: 1
##
## Estimates
## scale shape
## 20.3315 0.4267
##
## Standard Errors
## scale shape
## 8.3108 0.3552
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 25
## Gradient Evaluations: 15
fpot(temp, 13, r=3, model="pp",cmax=T) # 群れの最大のみでppモデル
##
## Call: fpot(x = temp, threshold = 13, model = "pp", cmax = T, r = 3)
## Deviance: 97.82808
##
## Threshold: 13
## Number Above: 42
## Proportion Above: 0.084
##
## Clustering Interval: 3
## Number of Clusters: 20
## Extremal Index: 0.4762
##
## Estimates
## loc scale shape
## 119.8724 57.1187 0.3394
##
## Standard Errors
## loc scale shape
## 46.5858 45.3085 0.3449
##
## Optimization Information
## Convergence: iteration limit reached
## Function Evaluations: 158
## Gradient Evaluations: 100
dim(temp) <- c(25,20) # ブロック最大で調べる準備 20ブロックにする
temp2 <- apply(temp, 2, max) # 各ブロックの最大を取る
fgev(temp2) # ブロック最大のGEVモデル
##
## Call: fgev(x = temp2)
## Deviance: 185.4294
##
## Estimates
## loc scale shape
## 16.7558 13.2754 0.8758
##
## Standard Errors
## loc scale shape
## 3.6359 4.3373 0.3371
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 54
## Gradient Evaluations: 38
# 8.6節
library(evd)
library(ismev)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
data(dowjones) # データをロードする
# dowjones
plot(dowjones, type="l") # 時系列としてデータをグラフ化する

dj <- dowjones[,2]
dldj <- diff(log(dj)) # 対数差分を取る
plot(dowjones[2:1304,1], dldj, type="l") # 対数差分のグラフ

max(dldj)
## [1] 0.04860535
-min(dldj)
## [1] 0.07454905
length(dldj)
## [1] 1303
dldjm <- matrix(dldj[1:1300], 20, 65) # 後ろ3つを捨て20*65型行列にする
(dldjp <- apply(dldjm, 2, max)) # 各列の最大値を取る
## [1] 0.008944056 0.008753856 0.011532093 0.007997048 0.011720876 0.016628651
## [7] 0.020007092 0.010845452 0.011613771 0.007262826 0.013237362 0.016105784
## [13] 0.012986306 0.011469094 0.008024566 0.015748897 0.019792525 0.015647241
## [19] 0.014981749 0.014655216 0.020747639 0.026048668 0.017745915 0.020023686
## [25] 0.019405219 0.033206080 0.022384149 0.046008413 0.023994220 0.014620484
## [31] 0.015732055 0.025138889 0.014701341 0.013253654 0.012417665 0.018667895
## [37] 0.018768467 0.016283430 0.037535399 0.048605349 0.040646996 0.023171651
## [43] 0.017271696 0.024797411 0.022520844 0.027984254 0.018613661 0.020699032
## [49] 0.016108149 0.017770396 0.011090735 0.021462702 0.012404013 0.021662878
## [55] 0.016289916 0.022138190 0.023648100 0.018669970 0.048095195 0.027097229
## [61] 0.019856400 0.021886599 0.014633657 0.013722624 0.013361791
(dldjp.gev.e <- fgev(dldjp)) # evdでブロック最大GEVモデル
##
## Call: fgev(x = dldjp)
## Deviance: -455.7098
##
## Estimates
## loc scale shape
## 0.015274 0.005933 0.130391
##
## Standard Errors
## loc scale shape
## 0.0008154 0.0005410 0.0989643
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 88
## Gradient Evaluations: 9
oldpar <- par()
par(mfrow=c(2,2))
plot(dldjp.gev.e) # そのEDA

dldjp.gev.i <-gev.fit(dldjp) # ismevでブロック最大GEVモデル
## $conv
## [1] 0
##
## $nllh
## [1] -227.9024
##
## $mle
## [1] 0.015195615 0.005744125 0.131994411
##
## $se
## [1] 0.0007981760 0.0004965542 0.0949590451
gev.diag(dldjp.gev.i) # そのEDA

(dldjn <- -apply(dldjm, 2, min)) # 各列の最小値を取り正負を入れ替える
## [1] 0.009054944 0.010544110 0.006104115 0.019805679 0.019126694 0.011681518
## [7] 0.030822786 0.015697412 0.013898281 0.009272651 0.029661283 0.008268274
## [13] 0.011396220 0.003519488 0.007312398 0.004462789 0.015553314 0.015556984
## [19] 0.013300061 0.023061376 0.023583977 0.019519912 0.005827825 0.024967581
## [25] 0.019333221 0.031642328 0.017035719 0.074549045 0.021044824 0.016445389
## [31] 0.028890968 0.008971548 0.011176801 0.009397597 0.016347586 0.023709699
## [37] 0.011427796 0.034671660 0.065781933 0.032234183 0.009519996 0.023473635
## [43] 0.020604578 0.024757511 0.017160517 0.015282108 0.022358175 0.017608585
## [49] 0.022219221 0.012454829 0.017269780 0.016000850 0.021046077 0.026289894
## [55] 0.009517253 0.010602088 0.032173916 0.028461346 0.037514172 0.058217077
## [61] 0.023667780 0.020271787 0.025091906 0.017297387 0.010044842
(dldjn.gev.e <- fgev(dldjn)) # evdでブロック最大GEVモデル
##
## Call: fgev(x = dldjn)
## Deviance: -415.2338
##
## Estimates
## loc scale shape
## 0.014389 0.007894 0.155708
##
## Standard Errors
## loc scale shape
## 0.001097 0.000787 0.098860
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 53
## Gradient Evaluations: 8
plot(dldjn.gev.e) # そのEDA

dldjn.gev.i <- gev.fit(dldjn) # ismevでブロック最大GEVモデル
## $conv
## [1] 0
##
## $nllh
## [1] -207.6322
##
## $mle
## [1] 0.014328890 0.007749394 0.156441392
##
## $se
## [1] 0.001085028 0.000753832 0.096728304
gev.diag(dldjn.gev.i) # そのEDA

par(oldpar)
mrlplot(dldj) # 標本平均超過プロット

mrlplot(dldj, tlim=c(0,0.04))

mrlplot(dldj, tlim=c(0.005,0.03))

sum(dldj>0.008)
## [1] 286
# tcplot(dldj, tlim=c(0.005,0.03)) # 推定プロット
tcplot(dldj, tlim=c(0.005,0.03), method="Nelder-Mead")


sum(dldj>0.016)
## [1] 76
(dldjp.pot <- fpot(dldj, 0.016, r=3, cmax=T)) # 母数の推定
##
## Call: fpot(x = dldj, threshold = 0.016, cmax = T, r = 3)
## Deviance: -525.9453
##
## Threshold: 0.016
## Number Above: 76
## Proportion Above: 0.0583
##
## Clustering Interval: 3
## Number of Clusters: 65
## Extremal Index: 0.8553
##
## Estimates
## scale shape
## 6.437e-03 2.344e-14
##
## Standard Errors
## scale shape
## 0.001112 0.104221
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 23
## Gradient Evaluations: 1
par(mfrow=c(2,2))
plot(dldjp.pot) # EDA

(dldjp.pot <- fpot(dldj, 0.016, r=3, cmax=T, method="Nelder-Mead"))
##
## Call: fpot(x = dldj, threshold = 0.016, cmax = T, r = 3, method = "Nelder-Mead")
## Deviance: -527.6451
##
## Threshold: 0.016
## Number Above: 76
## Proportion Above: 0.0583
##
## Clustering Interval: 3
## Number of Clusters: 65
## Extremal Index: 0.8553
##
## Estimates
## scale shape
## 0.005348 0.172258
##
## Standard Errors
## scale shape
## 0.0009312 0.1474533
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 69
plot(dldjp.pot)

par(oldpar)
mrlplot(-dldj) # 標本平均超過プロット

mrlplot(-dldj, tlim=c(0,0.04))

mrlplot(-dldj, tlim=c(0.015,0.03))

tcplot(-dldj, tlim=c(0.015,0.03)) # 母数の推定プロット


sum(-dldj>0.018)
## [1] 53
(dldjn.pot <- fpot(-dldj, 0.018, r=3, cmax=T, method="Nelder-Mead"))
##
## Call: fpot(x = -dldj, threshold = 0.018, cmax = T, r = 3, method = "Nelder-Mead")
## Deviance: -307.0608
##
## Threshold: 0.018
## Number Above: 53
## Proportion Above: 0.0407
##
## Clustering Interval: 3
## Number of Clusters: 42
## Extremal Index: 0.7925
##
## Estimates
## scale shape
## 0.008101 0.160500
##
## Standard Errors
## scale shape
## 0.001704 0.157299
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 61
length(grep("1999", dowjones[,1])) # 1年は何日か?
## [1] 261
fpot(-dldj, 0.016, r=3, cmax=T, npp=261, mper=10, method="Nelder-Mead") # 母数の推定
##
## Call: fpot(x = -dldj, threshold = 0.016, npp = 261, cmax = T, r = 3, mper = 10, method = "Nelder-Mead")
## Deviance: -419.2387
##
## Threshold: 0.016
## Number Above: 73
## Proportion Above: 0.056
##
## Clustering Interval: 3
## Number of Clusters: 56
## Extremal Index: 0.7671
##
## Estimates
## rlevel shape
## 0.07259 0.20603
##
## Standard Errors
## rlevel shape
## 0.0191 0.1554
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 57
max(dldjn)
## [1] 0.07454905
# fgev(dldjn, prob=1/3650, method="Nelder-Mead") # 分位点の推定
# 計算できなかった
(temp <- fgev(dldjn, method="Nelder-Mead")$estimate) # 入力の省略のため
## loc scale shape
## 0.014328079 0.007748444 0.156375607
qgev(1-1/3650, loc=temp[1], scale=temp[2], shape=temp[3]) # 分位点の推定
## loc
## 0.1434644
# 8.7節
nh <- read.csv("NikkeiHeikin.csv")
dim(nh)
## [1] 848 6
head(nh)
## Date Open High Low Close Volume
## 1 1949-05-31 121.14 176.52 121.14 176.52 NA
## 2 1949-06-30 175.97 175.98 146.14 146.92 NA
## 3 1949-07-31 147.37 150.00 134.80 144.86 NA
## 4 1949-08-31 109.91 176.52 109.91 175.18 NA
## 5 1949-09-30 176.89 176.89 160.41 163.29 NA
## 6 1949-10-31 162.95 162.95 131.09 141.54 NA
plot(nh[,1], nh[,5], type="l") # 事前にデータを眺める

dlnh <- diff(log(nh[,5])) # 対数差分
plot(nh[2:848,1], dlnh, type="l") # 対数差分を眺める

max(dlnh)
## [1] 0.2300247
max(-dlnh)
## [1] 0.2721623
nh2 <- cbind(nh, c(0,dlnh))
head(nh2)
## Date Open High Low Close Volume c(0, dlnh)
## 1 1949-05-31 121.14 176.52 121.14 176.52 NA 0.00000000
## 2 1949-06-30 175.97 175.98 146.14 146.92 NA -0.18354596
## 3 1949-07-31 147.37 150.00 134.80 144.86 NA -0.01412046
## 4 1949-08-31 109.91 176.52 109.91 175.18 NA 0.19004626
## 5 1949-09-30 176.89 176.89 160.41 163.29 NA -0.07028626
## 6 1949-10-31 162.95 162.95 131.09 141.54 NA -0.14294540
nh3 <- nh2[9:848,] # 初めの年の5月から12月を捨てる
head(nh3)
## Date Open High Low Close Volume c(0, dlnh)
## 9 1950-01-31 108.56 108.56 92.54 92.54 NA -0.17202087
## 10 1950-02-28 92.56 113.40 92.56 107.44 NA 0.14929157
## 11 1950-03-31 108.19 108.19 97.61 97.61 NA -0.09595261
## 12 1950-04-30 109.91 114.55 88.70 99.31 NA 0.01726632
## 13 1950-05-31 103.78 103.78 93.61 94.22 NA -0.05261380
## 14 1950-06-30 94.22 94.68 86.17 86.17 NA -0.08931038
tail(nh3)
## Date Open High Low Close Volume c(0, dlnh)
## 843 2019-07-31 21566.27 21823.07 20993.44 21521.53 12048781700 0.01147791
## 844 2019-08-31 21361.58 21556.69 20110.76 20704.37 13904127400 -0.03870904
## 845 2019-09-30 20625.75 22253.97 20554.96 21755.84 14379994100 0.04953744
## 846 2019-10-31 21831.44 23008.43 21276.01 22927.04 13481080600 0.05243477
## 847 2019-11-30 22730.49 23595.87 22707.39 23293.91 13046528300 0.01587495
## 848 2019-12-31 23388.63 24076.00 23044.78 23821.11 9317493900 0.02238021
nh.mat <- matrix(nh3[,7], 12,70) # 1列12ヶ月70年分とする
(nh.pm <- apply(nh.mat, 2, max)) # 毎年の正の最大
## [1] 0.16023638 0.12457421 0.15259818 0.23002474 0.08568216 0.07918956
## [7] 0.09101456 0.08703709 0.09664070 0.07415778 0.09743517 0.09629403
## [13] 0.13517551 0.07299106 0.07986814 0.12726371 0.04973710 0.04669268
## [19] 0.07673670 0.06510909 0.05853917 0.10520903 0.08540286 0.07537376
## [25] 0.08857037 0.11317796 0.10205323 0.06812277 0.04934199 0.03453559
## [31] 0.04609770 0.04254525 0.08404979 0.05972421 0.08573951 0.03824395
## [37] 0.15072025 0.07621545 0.09525116 0.05907200 0.18284563 0.12555197
## [43] 0.12680708 0.11798024 0.14967048 0.13877011 0.06660659 0.06181652
## [49] 0.09281057 0.09735501 0.06960855 0.06942858 0.06443522 0.07843969
## [55] 0.05920442 0.08938565 0.05682098 0.02695128 0.10051475 0.12088807
## [61] 0.09092905 0.03702559 0.09948199 0.11154113 0.06178395 0.09301419
## [67] 0.06182364 0.07818058 0.05342911 0.05243477
(nh.nm <- -apply(nh.mat, 2, min)) # 毎年の負の最大
## [1] 0.172020866 0.048018156 0.040830815 0.244506309 0.087457628
## [6] 0.043124888 0.013457498 0.124922367 0.022180456 0.110328267
## [11] 0.098884158 0.137859139 0.094180942 0.124316256 0.053955862
## [16] 0.070430197 0.035896381 0.094891662 0.052741574 0.061845064
## [21] 0.177041914 0.147323522 -0.020116364 0.135558743 0.098817793
## [26] 0.064954024 0.044309379 0.035734661 0.017443733 0.022758963
## [31] 0.031336313 0.035716906 0.047890811 0.004882895 0.099599486
## [36] 0.049244389 0.054207264 0.108822304 0.030020611 0.039221880
## [41] 0.213489783 0.105924556 0.139936364 0.183062835 0.053006710
## [46] 0.089520005 0.085096660 0.109133082 0.149265102 0.035957719
## [51] 0.123548044 0.101730754 0.102106260 0.047797537 0.045994568
## [56] 0.058227431 0.088952135 0.065231372 0.272162319 0.102799587
## [61] 0.123916490 0.093512758 0.108407703 0.020657728 0.088295931
## [66] 0.085916223 0.101214908 0.014098287 0.110403246 0.077385450
max(nh.pm)
## [1] 0.2300247
max(nh.nm)
## [1] 0.2721623
(nhp.gev <- fgev(nh.pm)) # 正についてGEVモデル
##
## Call: fgev(x = nh.pm)
## Deviance: -278.0989
##
## Estimates
## loc scale shape
## 0.07170 0.02811 0.01469
##
## Standard Errors
## loc scale shape
## 0.003774 0.002732 0.086348
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 48
## Gradient Evaluations: 6
oldpar <- par()
par(mfrow=c(2,2))
plot(nhp.gev) # EDA

(nhn.gev <- fgev(nh.nm)) # 負についてGEVモデル
##
## Call: fgev(x = nh.nm)
## Deviance: -218.0339
##
## Estimates
## loc scale shape
## 0.06070 0.04421 -0.02697
##
## Standard Errors
## loc scale shape
## 0.005895 0.004204 0.081234
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 67
## Gradient Evaluations: 8
plot(nhn.gev) # EDA

mrlplot(dlnh, tlim=c(0,0.2)) # 標本平均超過プロット
mrlplot(dlnh, tlim=c(0.03,0.15)) # 範囲を限定して詳しく見る
tcplot(dlnh, tlim=c(0.03,0.15)) # 推定プロット

sum(dlnh>0.05) # 標本サイズ
## [1] 174
(nhp.pot1 <- fpot(dlnh, 0.05)) # GPモデル
##
## Call: fpot(x = dlnh, threshold = 0.05)
## Deviance: -885.7451
##
## Threshold: 0.05
## Number Above: 174
## Proportion Above: 0.2054
##
## Estimates
## scale shape
## 2.886e-02 3.423e-11
##
## Standard Errors
## scale shape
## 0.003164 0.075245
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 18
## Gradient Evaluations: 1
(nhp.pot1 <- fpot(dlnh, 0.05, method="Nelder-Mead")) # 数値計算を変える
##
## Call: fpot(x = dlnh, threshold = 0.05, method = "Nelder-Mead")
## Deviance: -885.9101
##
## Threshold: 0.05
## Number Above: 174
## Proportion Above: 0.2054
##
## Estimates
## scale shape
## 0.02791 0.03292
##
## Standard Errors
## scale shape
## 0.003129 0.083215
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 41
plot(nhp.pot1)

(nhp.pot2 <- fpot(dlnh, 0.05, r=2, cmax=T,method="Nelder-Mead")) # 群れを入れる
##
## Call: fpot(x = dlnh, threshold = 0.05, cmax = T, r = 2, method = "Nelder-Mead")
## Deviance: -491.9695
##
## Threshold: 0.05
## Number Above: 174
## Proportion Above: 0.2054
##
## Clustering Interval: 2
## Number of Clusters: 105
## Extremal Index: 0.6034
##
## Estimates
## scale shape
## 0.03797 -0.07206
##
## Standard Errors
## scale shape
## 0.005105 0.093117
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 45
plot(nhp.pot2)

(nhp.pot3 <- fpot(dlnh, 0.05, r=3, cmax=T,method="Nelder-Mead")) # 別の群れ
##
## Call: fpot(x = dlnh, threshold = 0.05, cmax = T, r = 3, method = "Nelder-Mead")
## Deviance: -328.2445
##
## Threshold: 0.05
## Number Above: 174
## Proportion Above: 0.2054
##
## Clustering Interval: 3
## Number of Clusters: 74
## Extremal Index: 0.4253
##
## Estimates
## scale shape
## 0.04616 -0.14239
##
## Standard Errors
## scale shape
## 0.006988 0.098747
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 57
plot(nhp.pot3)

mrlplot(-dlnh, tlim=c(0,0.2)) # 負のデータの標本平均超過プロット
mrlplot(-dlnh, tlim=c(0.05,0.13)) # 範囲を限定して詳しく見る
tcplot(-dlnh, tlim=c(0.05,0.13)) # 推定プロット

sum(-dlnh>0.08) # 標本サイズ
## [1] 71
(nhn.pot1 <- fpot(-dlnh, 0.08)) # GPモデル
##
## Call: fpot(x = -dlnh, threshold = 0.08)
## Deviance: -329.0275
##
## Threshold: 0.08
## Number Above: 71
## Proportion Above: 0.0838
##
## Estimates
## scale shape
## 0.03427 0.05662
##
## Standard Errors
## scale shape
## 0.006135 0.134642
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 59
## Gradient Evaluations: 6
plot(nhn.pot1) # EDA

(nhn.pot2 <- fpot(-dlnh, 0.08, r=2, cmax=T)) # 群れを入れる
##
## Call: fpot(x = -dlnh, threshold = 0.08, cmax = T, r = 2)
## Deviance: -216.728
##
## Threshold: 0.08
## Number Above: 71
## Proportion Above: 0.0838
##
## Clustering Interval: 2
## Number of Clusters: 49
## Extremal Index: 0.6901
##
## Estimates
## scale shape
## 0.03884 0.03714
##
## Standard Errors
## scale shape
## 0.008317 0.160036
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 34
## Gradient Evaluations: 5
plot(nhn.pot2) # EDA

(nhn.pot3 <- fpot(-dlnh, 0.08, r=3, cmax=T)) # 別の群れ
##
## Call: fpot(x = -dlnh, threshold = 0.08, cmax = T, r = 3)
## Deviance: -187.4666
##
## Threshold: 0.08
## Number Above: 71
## Proportion Above: 0.0838
##
## Clustering Interval: 3
## Number of Clusters: 44
## Extremal Index: 0.6197
##
## Estimates
## scale shape
## 0.04514 -0.03221
##
## Standard Errors
## scale shape
## 0.009887 0.159311
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 46
## Gradient Evaluations: 4
plot(nhn.pot3) # EDA

fpot(-dlnh, 0.08, r=3, cmax=T, npp=12, mper=10) # 10年で最大下落
##
## Call: fpot(x = -dlnh, threshold = 0.08, npp = 12, cmax = T, r = 3, mper = 10)
## Deviance: -187.4667
##
## Threshold: 0.08
## Number Above: 71
## Proportion Above: 0.0838
##
## Clustering Interval: 3
## Number of Clusters: 44
## Extremal Index: 0.6197
##
## Estimates
## rlevel shape
## 0.16021 -0.03183
##
## Standard Errors
## rlevel shape
## 0.0118 0.1594
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 19
## Gradient Evaluations: 5
par(oldpar)