閱讀全文：http://tecdat.cn/?p=17375

最近我們被客戶要求撰寫關(guān)于馬爾可夫鏈的研究報告，包括一些圖形和統(tǒng)計輸出。

為了幫助客戶使用POT模型，本指南包含有關(guān)使用此模型的實用示例。本文快速介紹了極值理論（EVT）、一些基本示例，最后則通過案例對河流的極值進(jìn)行了具體的統(tǒng)計分析?

EVT的介紹

單變量情況

假設(shè)存在歸一化常數(shù)an> 0和bn使得：

根據(jù)極值類型定理（Fisher和Tippett，1928年），G必須是Fr'echet，Gumbel或負(fù)Weibull分布。Jenkinson（1955）指出，這三個分布可以合并為一個參數(shù)族：廣義極值（GEV）分布。GEV具有以下定義的分布函數(shù)：

根據(jù)這一結(jié)果，Pickands（1975）指出，當(dāng)閾值接近目標(biāo)變量的端點μend時，閾值閾值的標(biāo)準(zhǔn)化超額的極限分布是廣義Pareto分布（GPD）。也就是說，如果X是一個隨機(jī)變量，則：

基本用法

隨機(jī)數(shù)和分布函數(shù)

首先，讓我們從基本的東西開始。將R用于隨機(jī)數(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é)果如圖所示。我們可以清楚地看到，將閾值設(shè)為0.98是合理的選擇。

可以將置信區(qū)間添加到該圖，因為經(jīng)驗均值可以被認(rèn)為是正態(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ù)集。此外，由于洪水?dāng)?shù)據(jù)是一個時間序列，因此具有很強(qiáng)的自相關(guān)性，因此我們必須“提取”極端事件，同時保持事件之間的獨立性。

5,?clust.max?=?TRUE) >?diplot(events,?u.range?=?c(2,?20))

色散指數(shù)圖如圖所示。從該圖可以看出，大約5的閾值是合理的。

點擊標(biāo)題查閱往期內(nèi)容

極值理論 EVT、POT超閾值、GARCH 模型分析股票指數(shù)VaR、條件CVaR：多元化投資組合預(yù)測風(fēng)險測度分析

左右滑動查看更多

01

02

03

04

擬合GPD

單變量情況

可以根據(jù)多個估算器擬合GPD。
MLE是一種特殊情況，因為它是唯一允許變化閾值的情況。
我們在此給出一些教學(xué)示例。

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估計允許比例或形狀參數(shù)。例如，如果我們要擬合指數(shù)分布：

>?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。所有這些模型均使用最大似然估計量進(jìn)行擬合。

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等人的χ統(tǒng)計量。（1999）。對于參數(shù)模型，我們有：

對于自變量，χ= 0，而對于完全依存關(guān)系，χ=1。在我們的應(yīng)用中，值0.02表示變量是獨立的{這是顯而易見的。從這個角度來看，可以固定一些參數(shù)。

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

注意，由于所有雙變量極值分布都漸近相關(guān)，因此Coles等人的χ統(tǒng)計量。（1999）始終等于1。
檢測依賴強(qiáng)度的另一種方法是繪制Pickands的依賴函數(shù)：

>?pickdep(Mlog)

水平線對應(yīng)于獨立性，而其他水平線對應(yīng)于完全相關(guān)性。請注意，通過構(gòu)造，混合模型和非對稱混合模型無法對完美的依賴度變量建模。

使用馬爾可夫鏈對依賴關(guān)系結(jié)構(gòu)進(jìn)行建模

超越的馬爾可夫鏈進(jìn)行超過閾值的峰分析的經(jīng)典方法是使GPD擬合最大值。但是，由于僅考慮群集最大值，因此存在數(shù)據(jù)浪費。主要思想是使用馬爾可夫鏈對依賴關(guān)系結(jié)構(gòu)進(jìn)行建模，而聯(lián)合分布顯然是多元極值分布。這個想法是史密斯等人首先提出的。（1997）。在本節(jié)的其余部分，我們將只關(guān)注一階馬爾可夫鏈。因此，所有超出的可能性為：

對于我們的應(yīng)用程序，我們模擬具有極值依賴結(jié)構(gòu)的一階馬爾可夫鏈。

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

置信區(qū)間

擬合統(tǒng)計模型后，通常會給出置信區(qū)間。當(dāng)前，只有mle，pwmu，pwmb，矩估計量可以計算置信區(qū)間。

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

對于形狀參數(shù)置信區(qū)間：

conf.inf?conf.sup 2.141919?NAconf.inf?conf.sup 0.05757576?0.26363636

分位數(shù)的置信區(qū)間也可用。

conf.inf?conf.sup 2.141919?NAconf.inf?conf.sup 0.05757576?0.26363636

?conf.inf?conf.sup 8.64712?12.22558conf.inf?conf.sup 8.944444?12.833333

conf.inf?conf.sup 8.64712?12.22558conf.inf?conf.sup 8.944444?12.833333

輪廓置信區(qū)間函數(shù)既返回置信區(qū)間，又繪制輪廓對數(shù)似然函數(shù)。

模型檢查

要檢查擬合的模型，用戶必須調(diào)用函數(shù)圖。

>?plot(fitted,?npy?=?1)

圖顯示了執(zhí)行獲得的圖形窗口。

聚類技術(shù)

在處理時間序列時，超過閾值的峰值可能會出現(xiàn)問題。確實，由于時間序列通常是高度自相關(guān)的，因此選擇高于閾值可能會導(dǎo)致相關(guān)事件。
該函數(shù)試圖在滿足獨立性標(biāo)準(zhǔn)的同時識別超過閾值的峰。為此，該函數(shù)至少需要兩個參數(shù)：閾值u和獨立性的時間條件。群集標(biāo)識如下：
1.第一次超出會啟動第一個集群；
2.在閾值以下的第一個觀察結(jié)果將“終止集群”，除非tim.cond不成立；
3.下一個超過tim.cond的集群將啟動新集群；
4.根據(jù)需要進(jìn)行迭代。
一項初步研究表明，如果兩個洪水事件不在8天之內(nèi)，則可以認(rèn)為兩個洪水事件是獨立的，請注意，定義tim.cond的單位必須與所分析的數(shù)據(jù)相同。
返回一個包含已識別集群的列表。

?clu(dat,?u?=?2,?tim.cond?=?8/365)

其他

返回周期

將返回周期轉(zhuǎn)換為非超出概率，反之亦然。它既需要返回期，也需要非超越概率。此外，必須指定每年平均事件數(shù)。

npy?retper?prob 1?1.8?50?0.9888889 npy?retper?prob 1?2.2?1.136364?0.6

無偏樣本L矩

函數(shù)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

河流閾值分析

在本節(jié)中，我們提供了對河流閾值的全面和詳細(xì)的分析。

時間序列的移動平均窗口

從初始時間序列ts計算“平均”時間序列。這是通過在初始時間序列上使用長度為d的移動平均窗口來實現(xiàn)的。

>?plot(dat,?type?=?"l",?col?=?"blue") >?lines(tsd,?col?=?"green")

執(zhí)行過程如圖所示。圖描繪了瞬時洪水?dāng)?shù)量-藍(lán)線。

由于這是一個時間序列，因此我們必須選擇一個閾值以上的獨立事件。首先，我們固定一個相對較低的閾值以“提取”更多事件。因此，其中一些不是極端事件而是常規(guī)事件。這對于為GPD的漸近逼近選擇一個合理的閾值是必要的。

>?par(mfrow?=?c(2,?2)) >?plot(even[,?"obs"]) >?abline(v?=?6,?col?=?"green"

根據(jù)圖，閾值6m3 = s應(yīng)該是合理的。平均剩余壽命圖-左上方面板-表示大約10m3 = s的閾值應(yīng)足夠。但是，所選閾值必須足夠低，以使其上方具有足夠的事件以減小方差，而所選閾值則不能過低，因為它會增加偏差3。
因此，我們現(xiàn)在可以\重新提取閾值6m3 = s以上的事件，以獲得對象event1。注意，可以使用極值指數(shù)的估計。

>?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

結(jié)果給出了估計器的名稱（閾值，閾值以上的觀測值的數(shù)量和比例，參數(shù)估計，標(biāo)準(zhǔn)誤差估計和類型，漸近方差-協(xié)方差矩陣和收斂性診斷。
圖顯示了擬合模型的圖形診斷。可以看出，擬合模型“ mle”似乎是合適的。假設(shè)我們想知道與100年返回期相關(guān)的返回水平。

>?rp2p npy?retper?prob 1?1.707897?100?0.9941448 >? scale 36.44331

考慮到不確定性，圖描繪了與100年返回期相關(guān)的分位數(shù)的分布置信區(qū)間。

conf.inf?conf.sup 25.56533?90.76633

有時有必要知道指定事件的估計返回周期。讓我們對較大事件進(jìn)行處理。

>?maxEvent [1]?44.2 shape 0.9974115 >?prob npy?retper?prob 1?1.707897?226.1982?0.9974115

因此，2000年6月發(fā)生的最大事件的重現(xiàn)期大約為240年。

conf.inf?conf.sup 25.56533?90.76633

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《R語言有極值（EVT）依賴結(jié)構(gòu)的馬爾可夫鏈(MC)對洪水極值分析》。

點擊標(biāo)題查閱往期內(nèi)容

R語言極值理論：希爾HILL統(tǒng)計量尾部指數(shù)參數(shù)估計可視化
極值理論 EVT、POT超閾值、GARCH 模型分析股票指數(shù)VaR、條件CVaR：多元化投資組合預(yù)測風(fēng)險測度分析
R語言POT超閾值模型和極值理論EVT分析
R語言極值推斷：廣義帕累托分布GPD使用極大似然估計、輪廓似然估計、Delta法
R語言極值理論EVT：基于GPD模型的火災(zāi)損失分布分析
R語言有極值（EVT）依賴結(jié)構(gòu)的馬爾可夫鏈(MC)對洪水極值分析
R語言POT超閾值模型和極值理論EVT分析
R語言混合正態(tài)分布極大似然估計和EM算法
R語言多項式線性模型：最大似然估計二次曲線
R語言Wald檢驗 vs 似然比檢驗
R語言GARCH-DCC模型和DCC（MVT）建模估計
R語言非參數(shù)方法：使用核回歸平滑估計和K-NN(K近鄰算法)分類預(yù)測心臟病數(shù)據(jù)
matlab實現(xiàn)MCMC的馬爾可夫轉(zhuǎn)換ARMA - GARCH模型估計
R語言基于Bootstrap的線性回歸預(yù)測置信區(qū)間估計方法
R語言隨機(jī)搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
Matlab馬爾可夫鏈蒙特卡羅法（MCMC）估計隨機(jī)波動率（SV，Stochastic Volatility） 模型
Matlab馬爾可夫區(qū)制轉(zhuǎn)換動態(tài)回歸模型估計GDP增長率R語言極值推斷：廣義帕累托分布GPD使用極大似然估計、輪廓似然估計、Delta法