共计 9372 个字符,预计需要花费 24 分钟才能阅读完成。
浏览全文: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 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
点击文末 “浏览原文”
获取全文残缺材料。
本文选自《 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 法