原文連結：http://tecdat.cn/?p=17375
原文出處：拓端資料部落公衆号
為了幫助客戶使用POT模型，本指南包含有關使用此模型的實用示例。本文快速介紹了極值理論（EVT）、一些基本示例，最後則通過案例對河流的極值進行了具體的統計分析。
EVT的介紹
單變量情況
假設存在歸一化常數an> 0和bn使得：
根據極值類型定理（Fisher和Tippett，1928年），G必須是Fr'echet，Gumbel或負Weibull分布。Jenkinson（1955）指出，這三個分布可以合并為一個參數族：廣義極值（GEV）分布。GEV具有以下定義的分布函數：
根據這一結果，Pickands（1975）指出，當門檻值接近目标變量的端點µend時，門檻值門檻值的标準化超額的極限分布是廣義Pareto分布（GPD）。也就是說，如果X是一個随機變量，則：
基本用法
随機數和分布函數
首先，讓我們從基本的東西開始。将R用于随機數生成和分布函數。
> rgpd(5, loc = 1, scale = 2, shape = -0.2)
[1] 1.523393 2.946398 2.517602 1.199393 2.541937
> rgpd(6, c(1, -5), 2, -0.2)
[1] 1.3336965 -4.6504749 3.1366697 -0.9330325 3.5152161 -4.4851408
> rgpd(6, 0, c(2, 3), 0)
[1] 3.1139689 6.5900384 0.1886106 0.9797699 3.2638614 5.4755026
> pgpd(c(9, 15, 20), 1, 2, 0.25)
[1] 0.9375000 0.9825149 0.9922927
> qgpd(c(0.25, 0.5, 0.75), 1, 2, 0)
[1] 1.575364 2.386294 3.772589
> dgpd(c(9, 15, 20), 1, 2, 0.25)
[1] 0.015625000 0.003179117 0.001141829
使用選項lower.tail = TRUE或lower.tail = FALSE分别計算不超過或超過機率；
指定分位數是否超過機率分别帶有選項lower.tail = TRUE或lower.tail = FALSE；
指定是分别使用選項log = FALSE還是log = TRUE計算密度或對數密度。
門檻值選擇圖
此外，可以使用Fisher資訊來計算置信區間。
> x <- runif(10000)
> par(mfrow = c(1, 2))
結果如圖所示。我們可以清楚地看到，将門檻值設為0.98是合理的選擇。
可以将置信區間添加到該圖，因為經驗均值可以被認為是正态分布的（中心極限定理）。但是，對于高門檻值，正态性不再成立，此外，通過構造，該圖始終會收斂到點（xmax; 0）。
這是另一個綜合示例。
> x <- rnorm(10000)
plot(x, u.range = c(1, quantile(x, probs = 0.995)), col =
L-矩圖
L-矩是機率分布和資料樣本的摘要統計量。它們類似于普通矩{它們提供位置，離散度，偏度，峰度以及機率分布或資料樣本形狀的其他方面的路徑成本{但是是從有序資料值的線性組合中計算出來的（是以有字首L）。
這是一個簡單的例子。
> x <- c(1 - abs(rnorm(200, 0, 0.2)), rgpd(100, 1, 2, 0.25))
我們發現該圖形在真實資料上的性能通常很差。
色散指數圖
在處理時間序列時，色散指數圖特别有用。EVT指出，超出門檻值的超出部分可以通過GPD近似。但是，EVT必須通過蔔瓦松過程來表示這些超額部分的發生。
對于下一個示例，我們使用POT包中包含的資料集。此外，由于洪水資料是一個時間序列，是以具有很強的自相關性，是以我們必須“提取”極端事件，同時保持事件之間的獨立性。
5, clust.max = TRUE)
> diplot(events, u.range = c(2, 20))
色散指數圖如圖所示。從該圖可以看出，大約5的門檻值是合理的。
拟合GPD
單變量情況
可以根據多個估算器拟合GPD。
MLE是一種特殊情況，因為它是唯一允許變化門檻值的情況。
我們在此給出一些教學示例。
scale shape
mom 1.9538076495 0.2423393
mle 2.0345084386 0.2053905
pwmu 2.0517348996 0.2043644
pwmb 2.0624399910 0.2002131
pickands 2.3693985422 -0.0708419
med 2.2194363549 0.1537701
mdpd 2.0732577511 0.1809110
mple 2.0499646631 0.1960452
ad2r 0.0005539296 27.5964097
MLE，MPLE和MGF估計允許比例或形狀參數。例如，如果我們要拟合指數分布：
> fit(x, thresh = 1, shape = 0, est = "mle")
Estimator: MLE
Deviance: 322.686
AIC: 324.686
Varying Threshold: FALSE
Threshold Call: 1
Number Above: 100
Proportion Above: 1
Estimates
scale
1.847
Standard Error Type: observed
Standard Errors
scale
0.1847
Asymptotic Variance Covariance
scale
scale 0.03410
Optimization Information
Convergence: successful
Function Evaluations: 7
Gradient Evaluations: 1
> fitgpd(x, thresh = 1, scale = 2, est = "mle")
Estimator: MLE
Deviance: 323.3049
AIC: 325.3049
Varying Threshold: FALSE
Threshold Call: 1
Number Above: 100
Proportion Above: 1
Estimates
shape
0.0003398
Standard Error Type: observed
Standard Errors
shape
0.06685
Asymptotic Variance Covariance
shape
shape 0.004469
Optimization Information
Convergence: successful
Function Evaluations: 5
Gradient Evaluations: 1
If now, we want to fit a GPD with a varying threshold, just do:
> x <- rgpd(500, 1:2, 0.3, 0.01)
> fitgpd(x, 1:2, est = "mle")
Estimator: MLE
Deviance: -176.1669
AIC: -172.1669
Varying Threshold: TRUE
Threshold Call: 1:2
Number Above: 500
Proportion Above: 1
Estimates
scale shape
0.3261 -0.0556
Standard Error Type: observed
Standard Errors
scale shape
0.02098 0.04632
Asymptotic Variance Covariance
scale shape
scale 0.0004401 -0.0007338
shape -0.0007338 0.0021451
Optimization Information
Convergence: successful
Function Evaluations: 62
Gradient Evaluations: 11
| scale1 | shape1 | scale2 | shape2 | α |
| 6.784e-02 | 5.303e-02 | 2.993e-02 | 3.718e-02 | 2.001e-06 |
Asymptotic Variance Covariance
scale1 shape1 scale2 shape2 alpha
scale1 4.602e-03 -2.285e-03 1.520e-06 -1.145e-06 -3.074e-11
shape1 -2.285e-03 2.812e-03 -1.337e-07 4.294e-07 -1.843e-11
scale2 1.520e-06 -1.337e-07 8.956e-04 -9.319e-04 8.209e-12
shape2 -1.145e-06 4.294e-07 -9.319e-04 1.382e-03 5.203e-12
alpha -3.074e-11 -1.843e-11 8.209e-12 5.203e-12 4.003e-12
Optimization Information
Convergence: successful
Function Evaluations: 150
Gradient Evaluations: 21
雙變量情況
拟合雙變量POT。所有這些模型均使用最大似然估計量進行拟合。
vgpd(cbind(x, y), c(0, 2), model = "log")
> Mlog
Estimator: MLE
Dependence Model and Strenght:
Model : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = NA
Deviance: 1313.260
AIC: 1323.260
Marginal Threshold: 0 2
Marginal Number Above: 500 500
Marginal Proportion Above: 1 1
Joint Number Above: 500
Joint Proportion Above: 1
Number of events such as {Y1 > u1} U {Y2 > u2}: 500
Estimates
scale1 shape1 scale2 shape2 alpha
0.9814 0.2357 0.5294 -0.2835 0.9993
Standard Errors
在摘要中，我們可以看到lim_u Pr [X_1> u | X_2> u] = 0.02。這是Coles等人的χ統計量。（1999）。對于參數模型，我們有：
對于自變量，χ= 0，而對于完全依存關系，χ=1。在我們的應用中，值0.02表示變量是獨立的{這是顯而易見的。從這個角度來看，可以固定一些參數。
vgpd(cbind(x, y), c(0, 2), model = "log"
Dependence Model and Strenght:
Model : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = NA
Deviance: 1312.842
AIC: 1320.842
Marginal Threshold: 0 2
Marginal Number Above: 500 500
Marginal Proportion Above: 1 1
Joint Number Above: 500
Joint Proportion Above: 1
Number of events such as {Y1 > u1} U {Y2 > u2}: 500
Estimates
scale1 shape1 scale2 shape2
0.9972 0.2453 0.5252 -0.2857
Standard Errors
scale1 shape1 scale2 shape2
0.07058 0.05595 0.02903 0.03490
Asymptotic Variance Covariance
scale1 shape1 scale2 shape2
scale1 4.982e-03 -2.512e-03 -3.004e-13 3.544e-13
shape1 -2.512e-03 3.130e-03 1.961e-13 -2.239e-13
scale2 -3.004e-13 1.961e-13 8.427e-04 -8.542e-04
shape2 3.544e-13 -2.239e-13 -8.542e-04 1.218e-03
Optimization Information
Convergence: successful
Function Evaluations: 35
Gradient Evaluations: 9
注意，由于所有雙變量極值分布都漸近相關，是以Coles等人的χ統計量。（1999）始終等于1。
檢測依賴強度的另一種方法是繪制Pickands的依賴函數：
> pickdep(Mlog)
水準線對應于獨立性，而其他水準線對應于完全相關性。請注意，通過構造，混合模型和非對稱混合模型無法對完美的依賴度變量模組化。
使用馬爾可夫鍊對依賴關系結構進行模組化
超越的馬爾可夫鍊進行超過門檻值的峰分析的經典方法是使GPD拟合最大值。但是，由于僅考慮群集最大值，是以存在資料浪費。主要思想是使用馬爾可夫鍊對依賴關系結構進行模組化，而聯合分布顯然是多元極值分布。這個想法是史密斯等人首先提出的。（1997）。在本節的其餘部分，我們将隻關注一階馬爾可夫鍊。是以，所有超出的可能性為：
對于我們的應用程式，我們模拟具有極值依賴結構的一階馬爾可夫鍊。
mcgpd(mc, 2, "log")
Estimator: MLE
Dependence Model and Strenght:
Model : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = NA
Deviance: 1334.847
AIC: 1340.847
Threshold Call:
Number Above: 998
Proportion Above: 1
Estimates
scale shape alpha
0.8723 0.1400 0.5227
Standard Errors
scale shape alpha
0.08291 0.04730 0.02345
Asymptotic Variance Covariance
scale shape alpha
scale 0.0068737 -0.0016808 -0.0009066
shape -0.0016808 0.0022369 -0.0004412
alpha -0.0009066 -0.0004412 0.0005501
Optimization Information
Convergence: successful
Function Evaluations: 70
Gradient Evaluations: 15
置信區間
拟合統計模型後，通常會給出置信區間。目前，隻有mle，pwmu，pwmb，矩估計量可以計算置信區間。
conf.inf.scale conf.sup.scale
2.110881 2.939282
conf.inf.scale conf.sup.scale
1.633362 3.126838
conf.inf.scale conf.sup.scale
2.138851 3.074945
conf.inf.scale conf.sup.scale
2.149123 3.089090
對于形狀參數置信區間：
conf.inf conf.sup
2.141919 NA
conf.inf conf.sup
0.05757576 0.26363636
分位數的置信區間也可用。
conf.inf conf.sup
2.141919 NA
conf.inf conf.sup
0.05757576 0.26363636
conf.inf conf.sup
8.64712 12.22558
conf.inf conf.sup
8.944444 12.833333
conf.inf conf.sup
8.64712 12.22558
conf.inf conf.sup
8.944444 12.833333
輪廓置信區間函數既傳回置信區間，又繪制輪廓對數似然函數。
模型檢查
要檢查拟合的模型，使用者必須調用函數圖。
> plot(fitted, npy = 1)
圖顯示了執行獲得的圖形視窗。
聚類技術
在處理時間序列時，超過門檻值的峰值可能會出現問題。确實，由于時間序列通常是高度自相關的，是以選擇高于門檻值可能會導緻相關事件。
該函數試圖在滿足獨立性标準的同時識别超過門檻值的峰。為此，該函數至少需要兩個參數：門檻值u和獨立性的時間條件。群集辨別如下：
1.第一次超出會啟動第一個叢集；
2.在門檻值以下的第一個觀察結果将“終止叢集”，除非tim.cond不成立；
3.下一個超過tim.cond的叢集将啟動新叢集；
4.根據需要進行疊代。
一項初步研究表明，如果兩個洪水事件不在8天之内，則可以認為兩個洪水事件是獨立的，請注意，定義tim.cond的機關必須與所分析的資料相同。
傳回一個包含已識别叢集的清單。
clu(dat, u = 2, tim.cond = 8/365)
其他
傳回周期
将傳回周期轉換為非超出機率，反之亦然。它既需要傳回期，也需要非超越機率。此外，必須指定每年平均事件數。
npy retper prob
1 1.8 50 0.9888889
npy retper prob
1 2.2 1.136364 0.6
無偏樣本L矩
函數samlmu計算無偏樣本L矩。
| l_1 | l_2 | t_3 | t_4 | t_5 |
| 0.455381591 | 0.170423740 | 0.043928262 -0.005645249 -0.009310069 | ||
| 3.7.3 |
河流門檻值分析
在本節中，我們提供了對河流門檻值的全面和詳細的分析。
時間序列的移動平均視窗
從初始時間序列ts計算“平均”時間序列。這是通過在初始時間序列上使用長度為d的移動平均視窗來實作的。
> plot(dat, type = "l", col = "blue")
> lines(tsd, col = "green")
執行過程如圖所示。圖描繪了瞬時洪水數量-藍線。
由于這是一個時間序列，是以我們必須選擇一個門檻值以上的獨立事件。首先，我們固定一個相對較低的門檻值以“提取”更多事件。是以，其中一些不是極端事件而是正常事件。這對于為GPD的漸近逼近選擇一個合理的門檻值是必要的。
> par(mfrow = c(2, 2))
> plot(even[, "obs"])
> abline(v = 6, col = "green"
根據圖，門檻值6m3 = s應該是合理的。平均剩餘壽命圖-左上方面闆-表示大約10m3 = s的門檻值應足夠。但是，所選門檻值必須足夠低，以使其上方具有足夠的事件以減小方差，而所選門檻值則不能過低，因為它會增加偏差3。
是以，我們現在可以\重新提取門檻值6m3 = s以上的事件，以獲得對象event1。注意，可以使用極值指數的估計。
> events1 <-365, clust.max = TRUE)
> np
time obs
Min. :1970 Min. : 0.022
1st Qu.:1981 1st Qu.: 0.236
Median :1991 Median : 0.542
Mean :1989 Mean : 1.024
3rd Qu.:1997 3rd Qu.: 1.230
Max. :2004 Max. :44.200
NA's : 1.000
結果給出了估計器的名稱（門檻值，門檻值以上的觀測值的數量和比例，參數估計，标準誤差估計和類型，漸近方差-協方差矩陣和收斂性診斷。
圖顯示了拟合模型的圖形診斷。可以看出，拟合模型“ mle”似乎是合适的。假設我們想知道與100年傳回期相關的傳回水準。
> rp2p
npy retper prob
1 1.707897 100 0.9941448
>
scale
36.44331
考慮到不确定性，圖描繪了與100年傳回期相關的分位數的分布置信區間。
conf.inf conf.sup
25.56533 90.76633
有時有必要知道指定事件的估計傳回周期。讓我們對較大事件進行處理。
> maxEvent
[1] 44.2
shape
0.9974115
> prob
npy retper prob
1 1.707897 226.1982 0.9974115
是以，2000年6月發生的最大事件的重制期大約為240年。
conf.inf conf.sup
25.56533 90.76633
最受歡迎的見解
1.R語言基于ARMA-GARCH-VaR模型拟合和預測實證研究
2.R語言時變參數VAR随機模型
3.R語言時變參數VAR随機模型
4.R語言基于ARMA-GARCH過程的VAR拟合和預測
5.GARCH（1,1），MA以及曆史模拟法的VaR比較
6.R語言時變參數VAR随機模型
7.R語言實作向量自動回歸VAR模型
8.R語言随機搜尋變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
9.R語言VAR模型的不同類型的脈沖響應分析