浏览全文：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.5964097MLE，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: 11scale1
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 Covariancescale1 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.833333conf.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法