浏览全文 ： 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的阈值是正当的。

点击题目查阅往期内容

极值实践 EVT、POT超阈值、GARCH 模型剖析股票指数VaR、条件CVaR：多元化投资组合预测危险测度剖析

左右滑动查看更多

01

02

03

04

拟合GPD

单变量状况

能够依据多个估算器拟合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预计容许比例或形态参数。例如，如果咱们要拟合指数分布：

> 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

scale1

shape1

scale2

shape2

6.784e-02

5.303e-02

2.993e-02

3.718e-02

2.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等人的统计量。（1999）。对于参数模型，咱们有：

对于自变量，= 0，而对于齐全依存关系，=1。在咱们的利用中，值0.02示意变量是独立的{这是不言而喻的。从这个角度来看，能够固定一些参数。

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

留神，因为所有双变量极值散布都渐近相干，因而Coles等人的统计量。（1999）始终等于1。
检测依赖强度的另一种办法是绘制Pickands的依赖函数：

> pickdep(Mlog)

水平线对应于独立性，而其余水平线对应于齐全相关性。请留神，通过结构，混合模型和非对称混合模型无奈对完满的依赖度变量建模。

应用马尔可夫链对依赖关系构造进行建模

超过的马尔可夫链进行超过阈值的峰剖析的经典办法是使GPD拟合最大值。然而，因为仅思考群集最大值，因而存在数据节约。次要思维是应用马尔可夫链对依赖关系构造进行建模，而联结散布显然是多元极值散布。这个想法是史密斯等人首先提出的。（1997）。在本节的其余部分，咱们将只关注一阶马尔可夫链。因而，所有超出的可能性为：

对于咱们的应用程序，咱们模仿具备极值依赖构造的一阶马尔可夫链。

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

置信区间

拟合统计模型后，通常会给出置信区间。以后，只有mle，pwmu，pwmb，矩估计量能够计算置信区间。

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

对于形态参数置信区间：

conf.inf conf.sup2.141919 NA
conf.inf conf.sup0.05757576 0.26363636

分位数的置信区间也可用。

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

轮廓置信区间函数既返回置信区间，又绘制轮廓对数似然函数。

模型查看

要查看拟合的模型，用户必须调用函数图。

> 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 prob1 1.8 50 0.9888889npy retper prob1 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)> 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

后果给出了预计器的名称（阈值，阈值以上的观测值的数量和比例，参数估计，标准误差预计和类型，渐近方差-协方差矩阵和收敛性诊断。
图显示了拟合模型的图形诊断。能够看出，拟合模型“ mle”仿佛是适合的。假如咱们想晓得与100年返回期相干的返回程度。

> rp2pnpy retper prob1 1.707897 100 0.9941448> scale36.44331

思考到不确定性，图描述了与100年返回期相干的分位数的散布置信区间。

conf.inf conf.sup25.56533 90.76633

有时有必要晓得指定事件的预计返回周期。让咱们对较大事件进行解决。

> maxEvent[1] 44.2shape0.9974115> probnpy retper prob1 1.707897 226.1982 0.9974115

因而，2000年6月产生的最大事件的重现期大概为240年。

conf.inf conf.sup25.56533 90.76633

点击文末 “浏览原文”

获取全文残缺材料。

本文选自《R语言有极值（EVT）依赖构造的马尔可夫链(MC)对洪水极值剖析》。

点击题目查阅往期内容

R语言极值实践：希尔HILL统计量尾部指数参数估计可视化
极值实践 EVT、POT超阈值、GARCH 模型剖析股票指数VaR、条件CVaR：多元化投资组合预测危险测度剖析
R语言POT超阈值模型和极值实践EVT剖析
R语言极值推断：狭义帕累托散布GPD应用极大似然预计、轮廓似然预计、Delta法
R语言极值实践EVT：基于GPD模型的火灾损失散布剖析
R语言有极值（EVT）依赖构造的马尔可夫链(MC)对洪水极值剖析
R语言POT超阈值模型和极值实践EVT剖析
R语言混合正态分布极大似然预计和EM算法
R语言多项式线性模型：最大似然预计二次曲线
R语言Wald测验 vs 似然比测验
R语言GARCH-DCC模型和DCC（MVT）建模预计
R语言非参数办法：应用核回归平滑预计和K-NN(K近邻算法)分类预测心脏病数据
matlab实现MCMC的马尔可夫转换ARMA - GARCH模型预计
R语言基于Bootstrap的线性回归预测置信区间预计办法
R语言随机搜寻变量抉择SSVS预计贝叶斯向量自回归（BVAR）模型
Matlab马尔可夫链蒙特卡罗法（MCMC）预计随机稳定率（SV，Stochastic Volatility） 模型
Matlab马尔可夫区制转换动静回归模型预计GDP增长率R语言极值推断：狭义帕累托散布GPD应用极大似然预计、轮廓似然预计、Delta法