关于数据挖掘:R语言有极值EVT依赖结构的马尔可夫链MC对洪水极值分析

原文链接：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 模型的不同类型的脉冲响应剖析