# 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)