R語言有極值（EVT）依賴結(jié)構(gòu)的馬爾可夫鏈(MC)對洪水極值分析
?原文鏈接：http://tecdat.cn/?p=17375
?
為了幫助客戶使用POT模型，本指南包含有關(guān)使用此模型的實用示例。本文快速介紹了極值理論（EVT）、一些基本示例，最后則通過案例對河流的極值進行了具體的統(tǒng)計分析。
EVT的介紹
單變量情況
假設存在歸一化常數(shù)an> 0和bn使得：
?
根據(jù)極值類型定理（Fisher和Tippett，1928年），G必須是Fr'echet，Gumbel或負Weibull分布。Jenkinson（1955）指出，這三個分布可以合并為一個參數(shù)族：廣義極值（GEV）分布。GEV具有以下定義的分布函數(shù)：
根據(jù)這一結(jié)果，Pickands（1975）指出，當閾值接近目標變量的端點μend時，閾值閾值的標準化超額的極限分布是廣義Pareto分布（GPD）。也就是說，如果X是一個隨機變量，則：
基本用法
隨機數(shù)和分布函數(shù)
首先，讓我們從基本的東西開始。將R用于隨機數(shù)生成和分布函數(shù)。
> 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分別計算不超過或超過概率；
指定分位數(shù)是否超過概率分別帶有選項lower.tail = TRUE或lower.tail = FALSE；
指定是分別使用選項log = FALSE還是log = TRUE計算密度或?qū)?shù)密度。
?
閾值選擇圖
?
此外，可以使用Fisher信息來計算置信區(qū)間。
?
> x <- runif(10000)> par(mfrow = c(1, 2))
結(jié)果如圖所示。我們可以清楚地看到，將閾值設為0.98是合理的選擇。
?
?可以將置信區(qū)間添加到該圖，因為經(jīng)驗均值可以被認為是正態(tài)分布的（中心極限定理）。但是，對于高閾值，正態(tài)性不再成立，此外，通過構(gòu)造，該圖始終會收斂到點（xmax; 0）。
這是另一個綜合示例。
> x <- rnorm(10000)plot(x, u.range = c(1, quantile(x, probs = 0.995)), col =
?L-矩圖
?
L-矩是概率分布和數(shù)據(jù)樣本的摘要統(tǒng)計量。它們類似于普通矩{它們提供位置，離散度，偏度，峰度以及概率分布或數(shù)據(jù)樣本形狀的其他方面的度量值{但是是從有序數(shù)據(jù)值的線性組合中計算出來的（因此有前綴L）。
?
這是一個簡單的例子。
> x <- c(1 - abs(rnorm(200, 0, 0.2)), rgpd(100, 1, 2, 0.25))
我們發(fā)現(xiàn)該圖形在真實數(shù)據(jù)上的性能通常很差。
?
色散指數(shù)圖
?
在處理時間序列時，色散指數(shù)圖特別有用。EVT指出，超出閾值的超出部分可以通過GPD近似。但是，EVT必須通過泊松過程來表示這些超額部分的發(fā)生。
?
對于下一個示例，我們使用POT包中包含的數(shù)據(jù)集。此外，由于洪水數(shù)據(jù)是一個時間序列，因此具有很強的自相關(guān)性，因此我們必須“提取”極端事件，同時保持事件之間的獨立性。
5, clust.max = TRUE)> diplot(events, u.range = c(2, 20))
色散指數(shù)圖如圖所示。從該圖可以看出，大約5的閾值是合理的。
?
?
?
擬合GPD
單變量情況
?
可以根據(jù)多個估算器擬合GPD。
MLE是一種特殊情況，因為它是唯一允許變化閾值的情況。
我們在此給出一些教學示例。
scale shapemom 1.9538076495 0.2423393mle 2.0345084386 0.2053905pwmu 2.0517348996 0.2043644pwmb 2.0624399910 0.2002131pickands 2.3693985422 -0.0708419med 2.2194363549 0.1537701mdpd 2.0732577511 0.1809110mple 2.0499646631 0.1960452ad2r 0.0005539296 27.5964097
MLE，MPLE和MGF估計允許比例或形狀參數(shù)。例如，如果我們要擬合指數(shù)分布：
> fit(x, thresh = 1, shape = 0, est = "mle")Estimator: MLEDeviance: 322.686AIC: 324.686Varying Threshold: FALSEThreshold Call: 1Number Above: 100Proportion Above: 1Estimatesscale1.847Standard Error Type: observedStandard Errorsscale0.1847Asymptotic Variance Covariancescalescale 0.03410Optimization InformationConvergence: successfulFunction Evaluations: 7Gradient Evaluations: 1> fitgpd(x, thresh = 1, scale = 2, est = "mle")Estimator: MLEDeviance: 323.3049AIC: 325.3049Varying Threshold: FALSEThreshold Call: 1Number Above: 100Proportion Above: 1Estimatesshape0.0003398Standard Error Type: observedStandard Errorsshape0.06685Asymptotic Variance Covarianceshapeshape 0.004469Optimization InformationConvergence: successfulFunction Evaluations: 5Gradient Evaluations: 1If 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: MLEDeviance: -176.1669AIC: -172.1669Varying Threshold: TRUEThreshold Call: 1:2Number Above: 500Proportion Above: 1Estimatesscale shape0.3261 -0.0556Standard Error Type: observedStandard Errorsscale shape0.02098 0.04632Asymptotic Variance Covariancescale shapescale 0.0004401 -0.0007338shape -0.0007338 0.0021451Optimization InformationConvergence: successfulFunction Evaluations: 62Gradient Evaluations: 11
scale1shape1scale2shape2α6.784e-025.303e-022.993e-023.718e-022.001e-06
?
Asymptotic Variance Covariancescale1 shape1 scale2 shape2 alphascale1 4.602e-03 -2.285e-03 1.520e-06 -1.145e-06 -3.074e-11shape1 -2.285e-03 2.812e-03 -1.337e-07 4.294e-07 -1.843e-11scale2 1.520e-06 -1.337e-07 8.956e-04 -9.319e-04 8.209e-12shape2 -1.145e-06 4.294e-07 -9.319e-04 1.382e-03 5.203e-12alpha -3.074e-11 -1.843e-11 8.209e-12 5.203e-12 4.003e-12Optimization InformationConvergence: successfulFunction Evaluations: 150Gradient Evaluations: 21
雙變量情況
擬合雙變量POT。所有這些模型均使用最大似然估計量進行擬合。
vgpd(cbind(x, y), c(0, 2), model = "log")> MlogEstimator: MLEDependence Model and Strenght:Model : Logisticlim_u Pr[ X_1 > u | X_2 > u] = NADeviance: 1313.260AIC: 1323.260Marginal Threshold: 0 2Marginal Number Above: 500 500Marginal Proportion Above: 1 1Joint Number Above: 500Joint Proportion Above: 1Number of events such as {Y1 > u1} U {Y2 > u2}: 500Estimatesscale1 shape1 scale2 shape2 alpha0.9814 0.2357 0.5294 -0.2835 0.9993Standard Errors
在摘要中，我們可以看到lim_u Pr [X_1> u | X_2> u] = 0.02。這是Coles等人的χ統(tǒng)計量。（1999）。對于參數(shù)模型，我們有：
對于自變量，χ= 0，而對于完全依存關(guān)系，χ=1。在我們的應用中，值0.02表示變量是獨立的{這是顯而易見的。從這個角度來看，可以固定一些參數(shù)。
vgpd(cbind(x, y), c(0, 2), model = "log"Dependence Model and Strenght:Model : Logisticlim_u Pr[ X_1 > u | X_2 > u] = NADeviance: 1312.842AIC: 1320.842Marginal Threshold: 0 2Marginal Number Above: 500 500Marginal Proportion Above: 1 1Joint Number Above: 500Joint Proportion Above: 1Number of events such as {Y1 > u1} U {Y2 > u2}: 500Estimatesscale1 shape1 scale2 shape20.9972 0.2453 0.5252 -0.2857Standard Errorsscale1 shape1 scale2 shape20.07058 0.05595 0.02903 0.03490Asymptotic Variance Covariance
?
scale1 shape1 scale2 shape2scale1 4.982e-03 -2.512e-03 -3.004e-13 3.544e-13shape1 -2.512e-03 3.130e-03 1.961e-13 -2.239e-13scale2 -3.004e-13 1.961e-13 8.427e-04 -8.542e-04shape2 3.544e-13 -2.239e-13 -8.542e-04 1.218e-03Optimization InformationConvergence: successfulFunction Evaluations: 35Gradient Evaluations: 9
注意，由于所有雙變量極值分布都漸近相關(guān)，因此Coles等人的χ統(tǒng)計量。（1999）始終等于1。
檢測依賴強度的另一種方法是繪制Pickands的依賴函數(shù)：
> pickdep(Mlog)
水平線對應于獨立性，而其他水平線對應于完全相關(guān)性。請注意，通過構(gòu)造，混合模型和非對稱混合模型無法對完美的依賴度變量建模。
使用馬爾可夫鏈對依賴關(guān)系結(jié)構(gòu)進行建模
超越的馬爾可夫鏈進行超過閾值的峰分析的經(jīng)典方法是使GPD擬合最大值。但是，由于僅考慮群集最大值，因此存在數(shù)據(jù)浪費。主要思想是使用馬爾可夫鏈對依賴關(guān)系結(jié)構(gòu)進行建模，而聯(lián)合分布顯然是多元極值分布。這個想法是史密斯等人首先提出的。（1997）。在本節(jié)的其余部分，我們將只關(guān)注一階馬爾可夫鏈。因此，所有超出的可能性為：
?對于我們的應用程序，我們模擬具有極值依賴結(jié)構(gòu)的一階馬爾可夫鏈。
mcgpd(mc, 2, "log")Estimator: MLEDependence Model and Strenght:Model : Logisticlim_u Pr[ X_1 > u | X_2 > u] = NADeviance: 1334.847AIC: 1340.847Threshold Call:Number Above: 998Proportion Above: 1Estimatesscale shape alpha0.8723 0.1400 0.5227Standard Errorsscale shape alpha0.08291 0.04730 0.02345Asymptotic Variance Covariancescale shape alphascale 0.0068737 -0.0016808 -0.0009066shape -0.0016808 0.0022369 -0.0004412alpha -0.0009066 -0.0004412 0.0005501Optimization InformationConvergence: successfulFunction Evaluations: 70Gradient Evaluations: 15
置信區(qū)間
擬合統(tǒng)計模型后，通常會給出置信區(qū)間。當前，只有mle，pwmu，pwmb，矩估計量可以計算置信區(qū)間。
conf.inf.scale conf.sup.scale2.110881 2.939282conf.inf.scale conf.sup.scale1.633362 3.126838conf.inf.scale conf.sup.scale2.138851 3.074945conf.inf.scale conf.sup.scale2.149123 3.089090
對于形狀參數(shù)置信區(qū)間：
conf.inf conf.sup2.141919 NA
?
conf.inf conf.sup0.05757576 0.26363636
分位數(shù)的置信區(qū)間也可用。
conf.inf conf.sup2.141919 NA
?
conf.inf conf.sup0.05757576 0.26363636
?conf.inf conf.sup8.64712 12.22558
?
conf.inf conf.sup8.944444 12.833333
conf.inf conf.sup8.64712 12.22558
?
conf.inf conf.sup8.944444 12.833333
輪廓置信區(qū)間函數(shù)既返回置信區(qū)間，又繪制輪廓對數(shù)似然函數(shù)。
模型檢查
要檢查擬合的模型，用戶必須調(diào)用函數(shù)圖。
> plot(fitted, npy = 1)
圖顯示了執(zhí)行獲得的圖形窗口。
?
?
聚類技術(shù)
在處理時間序列時，超過閾值的峰值可能會出現(xiàn)問題。確實，由于時間序列通常是高度自相關(guān)的，因此選擇高于閾值可能會導致相關(guān)事件。
該函數(shù)試圖在滿足獨立性標準的同時識別超過閾值的峰。為此，該函數(shù)至少需要兩個參數(shù)：閾值u和獨立性的時間條件。群集標識如下：
1.第一次超出會啟動第一個集群；
2.在閾值以下的第一個觀察結(jié)果將“終止集群”，除非tim.cond不成立；
3.下一個超過tim.cond的集群將啟動新集群；
4.根據(jù)需要進行迭代。
一項初步研究表明，如果兩個洪水事件不在8天之內(nèi)，則可以認為兩個洪水事件是獨立的，請注意，定義tim.cond的單位必須與所分析的數(shù)據(jù)相同。
返回一個包含已識別集群的列表。
clu(dat, u = 2, tim.cond = 8/365)
其他
返回周期
將返回周期轉(zhuǎn)換為非超出概率，反之亦然。它既需要返回期，也需要非超越概率。此外，必須指定每年平均事件數(shù)。
npy retper prob1 1.8 50 0.9888889npy retper prob1 2.2 1.136364 0.6
?
無偏樣本L矩
函數(shù)samlmu計算無偏樣本L矩。
l_1l_2t_3t_4t_50.4553815910.1704237400.043928262 -0.005645249 -0.0093100693.7.3?
?
?
河流閾值分析
在本節(jié)中，我們提供了對河流閾值的全面和詳細的分析。
時間序列的移動平均窗口
從初始時間序列ts計算“平均”時間序列。這是通過在初始時間序列上使用長度為d的移動平均窗口來實現(xiàn)的。
> plot(dat, type = "l", col = "blue")> lines(tsd, col = "green")
執(zhí)行過程如圖所示。圖描繪了瞬時洪水數(shù)量-藍線。
由于這是一個時間序列，因此我們必須選擇一個閾值以上的獨立事件。首先，我們固定一個相對較低的閾值以“提取”更多事件。因此，其中一些不是極端事件而是常規(guī)事件。這對于為GPD的漸近逼近選擇一個合理的閾值是必要的。
?
> par(mfrow = c(2, 2))> plot(even[, "obs"])> abline(v = 6, col = "green"
根據(jù)圖，閾值6m3 = s應該是合理的。平均剩余壽命圖-左上方面板-表示大約10m3 = s的閾值應足夠。但是，所選閾值必須足夠低，以使其上方具有足夠的事件以減小方差，而所選閾值則不能過低，因為它會增加偏差3。
因此，我們現(xiàn)在可以\重新提取閾值6m3 = s以上的事件，以獲得對象event1。注意，可以使用極值指數(shù)的估計。
> events1 <-365, clust.max = TRUE)> nptime obsMin. :1970 Min. : 0.0221st Qu.:1981 1st Qu.: 0.236Median :1991 Median : 0.542Mean :1989 Mean : 1.0243rd Qu.:1997 3rd Qu.: 1.230Max. :2004 Max. :44.200NA's : 1.000
?
結(jié)果給出了估計器的名稱（閾值，閾值以上的觀測值的數(shù)量和比例，參數(shù)估計，標準誤差估計和類型，漸近方差-協(xié)方差矩陣和收斂性診斷。
圖顯示了擬合模型的圖形診斷?？梢钥闯?，擬合模型“ mle”似乎是合適的。假設我們想知道與100年返回期相關(guān)的返回水平。
> rp2pnpy retper prob1 1.707897 100 0.9941448>scale36.44331
考慮到不確定性，圖描繪了與100年返回期相關(guān)的分位數(shù)的分布置信區(qū)間。
?
conf.inf conf.sup25.56533 90.76633
有時有必要知道指定事件的估計返回周期。讓我們對較大事件進行處理。
> maxEvent[1] 44.2shape0.9974115> probnpy retper prob1 1.707897 226.1982 0.9974115
因此，2000年6月發(fā)生的最大事件的重現(xiàn)期大約為240年。
?
conf.inf conf.sup25.56533 90.76633
最受歡迎的見解
1.R語言基于ARMA-GARCH-VaR模型擬合和預測實證研究
2.R語言時變參數(shù)VAR隨機模型
3.R語言時變參數(shù)VAR隨機模型
4.R語言基于ARMA-GARCH過程的VAR擬合和預測
5.GARCH（1,1），MA以及歷史模擬法的VaR比較
6.R語言時變參數(shù)VAR隨機模型
7.R語言實現(xiàn)向量自動回歸VAR模型
8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
9.R語言VAR模型的不同類型的脈沖響應分析
