R中的极值分析入门–第2部分:块最大值方法

嘿,

welcome to part 2 of our short introduction to 极值分析 using the extRemes package in R.

你们中有些人可能知道,有两种常见的方法可以进行实际的极值分析。今天,我们将重点介绍这两种方法中的第一种,称为块最大值方法。

这种对(时间)观测值的极端进行建模的方法是基于在一定长度的恒定序列内利用这些观测值的最大值或最小值。对于足够大的数量  n  已建立的区块的数量,这些区块的最终峰值  n  可以使用等长的块来对这些数据进行适当的分配。虽然基本上可以自由选择块大小,但必须在偏差(小块)和方差(大块)之间进行权衡。通常,通常选择序列的长度以对应于某个熟悉的时间段,大多数情况下是一年。年度最大值(或最小值)的最终矢量为calles“年度最大值(最小值)系列” or simply AMS.

根据Fisher-Tippett-Gnedenko定理,块最大值的分布可通过 广义极值分布.

The following code shows a short practical example of fitting a 广义极值分布 to a time series of precipitation 数据 using the extRemes package in R. The sample 数据 set features precipitation 数据 of Bärnkopf (Lower Austria) from 1971 through 2014 and is is provided by the hydrographic services of Austria via H. The function read_ehyd() for importing the 数据 set can be found at 使用R从eHYD读取数据.

# load required packages
library(extRemes)
library(xts)

# get 数据 from H
ehyd_url <- "http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=107540&file=2"
precipitation_xts <- read_ehyd(ehyd_url)

# derive AMS for maximum precipitation
ams <- apply.yearly(precipitation_xts, max)

# maximum-likelihood fitting of the GEV distribution
fit_mle <- fevd(as.vector(ams), method = "MLE", type="GEV")
# diagnostic plots
plot(fit_mle)
# return levels:
rl_mle <- return.level(fit_mle, conf = 0.05, return.period= c(2,5,10,20,50,100))

# fitting of GEV distribution based on L-moments estimation
fit_lmom <- fevd(as.vector(ams), method = "Lmoments", type="GEV")
# diagnostic plots
plot(fit_lmom)
# return levels:
rl_lmom <- return.level(fit_lmom, conf = 0.05, return.period= c(2,5,10,20,50,100))

# return level plots
par(mfcol=c(1,2))
# return level plot w/ MLE
plot(fit_mle, type="rl",
     main="Return Level Plot for Bärnkopf w/ MLE",
     ylim=c(0,200), pch=16)
loc <- as.numeric(return.level(fit_mle, conf = 0.05,return.period=100))
segments(100, 0, 100, loc, col= 'midnightblue',lty=6)
segments(0.01,loc,100, loc, col='midnightblue', lty=6)

# return level plot w/ LMOM
plot(fit_lmom, type="rl",
     main="Return Level Plot for Bärnkopf w/ L-Moments",
     ylim=c(0,200))
loc <- as.numeric(return.level(fit_lmom, conf = 0.05,return.period=100))
segments(100, 0, 100, loc, col= 'midnightblue',lty=6)
segments(0.01,loc,100, loc, col='midnightblue', lty=6)

# 比较收益水平
results <- t(data.frame(mle=as.numeric(rl_mle),
                        lmom=as.numeric(rl_lmom)))
colnames(results) <- c(2,5,10,20,50,100)
round(results,1)

比较收益水平
In this case, both results are quite similar. In most cases, L-moments estimation is more robust than maximum likelihood estimation. In addition to these classical estimation methods, extRemes offers Generalized Maximum Likelihood Estimation (GMLE, 马丁斯和施丁格,2000年)和贝叶斯估计方法(吉兰(Gilleland)和卡兹(Katz),2016年). Moreover, I have made the observation that maximum likelihood estimation works more reliable in other R packages in some cases (e.g. fExtremes, ismev).

参考:Coles,Stuart(2001)。极值统计建模简介。统计中的Springer系列。伦敦:施普林格。

关于作者

马蒂亚斯在维也纳自然资源与生命科学大学学习了环境信息管理,并获得了环境统计博士学位。他的论文的重点是罕见(极端)事件的统计建模,作为对关键基础设施进行漏洞评估的基础。他目前在奥地利国家气象和地球物理服务局(ZAMG)和BOKU大学山区风险工程研究所工作。他目前专注于(统计)不良天气事件和自然灾害以及减少灾害风险的评估。他的主要兴趣是环境现象的统计建模以及用于数据科学,地理信息和遥感的开源工具。

34条留言

您可以在这篇文章中发表评论。


  • 我想要你的帮助

    穆罕默德·阿里夫(Muhammad Arif) 5年前 回复


    • 当然。您可以指定您的要求吗?

      马蒂亚斯 5年前 回复


  • 请给我发送适合广义pareto分布的r编码,用mle方法估计参数估计

    穆罕默德·阿里夫(Muhammad Arif) 5年前 回复


  • 我的研究主题是“极端降水的非平稳频率分析”。但是我不知道我的足迹。我将如何继续工作?
    请帮我。
    谢谢你的好心!

    穆罕默德·阿里夫(Muhammad Arif) 5年前 回复


    • 我还写了一篇有关处理非平稳性的文章,您可以在此处找到更多信息:
      R中的极值分析入门–第4部分:处理趋势
      The fevd() function from the package extRemes is quite powerful and well documented, I recommend to have a look at its help page.

      马蒂亚斯 5年前 回复


  • 嗨,我目前正在撰写有关Garch-EVT方法的论文,以估计风险价值和预期的缺口。您可以将这种方法的R程序发送给我吗?谢谢

    加比 在4年前 回复


    • 嗨,盖比,

      我尚未使用GARCH-EVT-copulas,如果那样的话’s what you intend to do. I am pretty sure that there is no specific R package targeted at this topic. So you probably have to split your procedure into several steps and deal with each step individually. Stack overflow is most likely a very good place to get some information and see the workflow of other people working with this topic. For a start, also take a look at fGarch (which is part of the Rmetrics-suite, just as the excellent package fExtremes) or rugarch for GARCH modelling.
      In addition, you might want to have a look at the package GEVStableGarch which seemingly employs maximum likelihood extimation of GARCH models with a GEV distribution, maybe this is sufficient for your needs.

      问候,
      马蒂亚斯

      马蒂亚斯 在4年前 回复


      • 感谢您的回复。

        I’我只会使用garch和evt。该过程包括将garch模型应用于返回序列,然后从模型中提取残差。然后,将evt应用于从garch模型提取的残差。但是,我很难从残差中选择阈值,以便对其应用EVT-POT或最大块法。你能帮助我吗?

        问候
        加比

        加比 在4年前 回复


    • I’我感到如此轻松,我很震惊。

      丹佛 在4年前 回复


  • 你好

    您的网站非常有用,使该主题易于理解

    您如何将收益水平图上的点变成绿色?

    米沙 在4年前 回复


    • 你好

      unfortunately, the plotting methods for objects of type fevd are –尽管全面–几乎不可定制。因为我没有’t like the default plotting, I modified the plotting functions plot.fevd.mleplot.fevd.lmoments to support better visualization options within the extRemes framework.

      问候,
      马蒂亚斯

      马蒂亚斯 在4年前 回复


      • 马蒂亚斯

        我仍在学习R的阶段。

        您能否为我提供一些有关如何完成此操作的指导。谢谢

        米沙 在4年前 回复


        • 嗨,Misha,

          I have written the following function for creating a simple return level plot (base graphics) from fevd objects that have been fitted of type GEV.
          您可以控制点颜色和拟合函数的线条以及y轴上的刻度数。

          rlplot_gev <- function(fevd_obj,
                                 pointcolor = "firebrick",
                                 linecolor = "darkgreen",
                                 yticks = seq(0,225,25)){
            rperiods = c(1+1e-15, 1.25, 1.5, 2, 3,4,5, 7.5, 10, 15, 20,
                         30, 40, 50, 60, 75, 80, 100, 120, 200, 250)
            model <- fevd_obj$type
            if(model != "GEV"){
              stop("this plotting function only works for GEV")
            }
            pars <- fevd_obj$results$par
            loc <- pars["location"]
            scale <- pars["scale"]
            shape <- pars["shape"]
            xrl <- -1/(log(1 - 1/rperiods))
            yrl <- rlevd(rperiods, loc = loc, scale = scale, 
                         shape = shape)
            
            # empty plot
            plot(xrl, yrl, type = "n", log = "x", xlab = "Return Period",
                 ylab = "Return Level", lwd = 2, yaxt = "n", xlim = c(1,200),
                 main = "Return Level Plot")
            par(las = 1)
            axis(side = 2, at = yticks)
            
            # background grid
            abline(h = yticks,col = "lightgrey", lty = 2)
            abline(v = c(1,2,5,10,20,50,100,200), col = "lightgrey", lty = 2)
            
            # plot line
            lines(xrl, yrl, col = linecolor, lwd = 2)
            
            # ci lines for MLE base plot
            bds <- ci(fevd_obj, return.period = rperiods)
            lines(xrl, bds[, 1], col = linecolor, lty = 4)
            lines(xrl, bds[, 3], col = linecolor, lty = 4)
            
            # Weibull plotting position of points  
            n <- fevd_obj$n
            xp <- ppoints(n = n, a = 0)
            ytmp <- 数据grabber(fevd_obj)
            y <- c(ytmp[, 1])
            points(-1/log(xp), sort(y), pch = 16, col = pointcolor)
          }
          

          Please note that this is a very simple approach, I'm mainly re-using code from plot.fevd.mle().
          为了为阈值超出法创建类似的图,您必须对GP的参数和极值的绘制位置进行一些调整:

          sdat <- sort(y)
          points(-1/log(xp2)[sdat > u]/x$npy, sdat[sdat > u], pch="+", col = pointcolor)
          

          Using rlplot_gev(mlefit, pointcolor = "darkblue") on the code of my post should yield return level plots similar to the default plot(mlefit, type = "rl", xlim = c(1, 200), ylim = c(25, 225))

          问候,
          马蒂亚斯

          马蒂亚斯 在4年前 回复


  • 嘿!

    我目前正在evds上写我的学士论文。

    如果我以您的代码为例,对您来说还可以吗?

    卢卡斯 在4年前 回复


    • 嗨,卢卡斯,

      是的,当然。但是,如果您有足够的时间,不妨对上面介绍的软件包进行更深入的了解。
      这篇文章实际上是作为对该主题的介绍。

      问候,
      马蒂亚斯

      马蒂亚斯 在4年前 回复


  • 当建模极端事件(例如变量X)对商品价格波动(例如Y)的影响时,是否可以将GEV拟合到X,然后将其作为均值或波动率方程中的外生变量插入?如果没有,我该如何解决这个问题?

    怀古鲁 3年前 回复


    • 你好

      我不太确定我是否了解您打算做什么。
      您要回答的问题是什么?
      首先,对X对Y的影响进行建模听起来更像是标准回归模型。

      马蒂亚斯 3年前 回复


  • 你好!
    我有个问题:
    如何使用R和Gumbel分布来预测特定回报期的排放?
    我怎么知道退货期?
    谢谢!

    阿比尔·哈达德(Abeer Haddad) 3年前 回复


    • 你好

      基本上,您可以坚持上面的示例代码。
      If you specifically want to fit a Gumbel distribution, simply set the type argument in the fevd() function to type="Gumbel".
      Return levels can be obtained by using the function return.level().

      问候,
      马蒂亚斯

      马蒂亚斯 3年前 回复


  • 嗨,Matthias,我也有一个问题:
    我想使用块最大值方法来估计风险价值。如何使用EVT BMM做到这一点?从广义极值分布中,我需要估计形状,比例和loval参数,但是R中的哪个函数可以提供此结果?我应该在长度上使用不同的块“months”, “quartal” and “6 月”。然后,借助参数估算风险价值。

    问候,
    玛蒂娜(Martina)

    玛蒂娜(Martina) 3年前 回复


  • 您好Matthias,

    非常感谢这篇文章,它’对我的项目有很大帮助。

    您知道为什么ML估计可能适用于某些样本数据而不适用于其他样本数据的可能原因是什么?正如我在项目中发现的那样,不得不使用L-Moments。

    非常感谢,
    彼得

    彼得 3年前 回复


    • 嗨,彼得,

      apart from the fact that it is of course possible that MLE may fail to converge (due to the existence of an unbounded likelihood function), you might want to try the implementations of either fExtremes or ismev. I have found that they may produce sensible MLE estimates in case the extRemes implementation fails. I have not yet found the leisure to check the fevd()function thoroughly in this respect, since it is a real monster of a function, which is written in an overly way complex way, imho.

      最好的祝福,
      马蒂亚斯

      马蒂亚斯 3年前 回复


      • 感谢您对Matthias的回复,

        好吧,我不’我不认为这是由于存在无限可能性函数引起的,但是无论如何我的项目的L矩估计都应该是合适的。

        再次感谢,彼得

        彼得 3年前 回复


  • 您好Matthias,
    I’ve tried to make a “阈值超出法的类似收益水平图” (see your comment above concerning sdat < - sort(y)和so on), but it does not work. Would you be so kind and give me some advice?
    谢谢,
    提洛

    提洛 3年前 回复


    • 嗨Tilo,
      基本上,您只需要对上述功能进行一些修改。在接下来的几天中,我将使用基本图形发布我的代码解决方案。我一直想为返回水平图提供一个ggplot解决方案,但是我还没有’到目前为止尚无时间实施。

      马蒂亚斯 3年前 回复


  • 嘿,我’我很兴奋,我发现了这一点,
    但是我似乎无法下载示例数据集,请继续获取

    Error during wrapup: scan() expected 'a real', got 'Lücke'

    when running read_ehyd(ehyd_url) using code you provided at http://www.citystoragesrq.com/read-ehyd/

    我也尝试过使用Rcurl(那个可行,但是我也得到了错误

    Error in if (ipars2["scale"] <= 0) ipars2["scale"] <- 1e-08 : 
      missing value where TRUE/FALSE needed
    In addition: Warning messages:
    1: In mean.default(x) : argument is not numeric or logical: returning NA
    2: In var(x) : NAs introduced by coercion
    

    我什么时候一样'我用自己的数据集'我想这可能是一个格式化问题,如果您仅能提供xts对象的样例或以正确格式格式化数据集的代码,我应该能够找出来,在期待您的答复的同时,请阅读您的R-Help很棒的文档和示例。

    卡梅伦·萨姆森(Cameron Samson) 3年前 回复


    • 我通过使用Marielle Pinheiro和Richard Grotjahn撰写的“极值统计简介”找到了解决方案。 (http://grotjahn.ucdavis.edu/EWEs/extremes_primer_v9_22_15.pdf),它基于Extremes软件包,并且对该软件包及其格式进行了介绍,

      最后是数据格式,我设法使其工作。

      再次感谢您的博客,尽管它对您有所帮助,

      卡梅伦

      卡梅伦 3年前 回复


    • 嗨,卡梅隆,

      您描述的问题是由于编码问题所致。数据文件包含德语单词“Lücke”’ (meaning `gap’ or `blank space’),表示NA值。
      麻烦是由Umlautü引起的—您可以在以下位置找到有关该问题的一些信息 Mojibake在Wikipedia上的文章.
      I have set the target file encoding in the read_ehyd()-function to latin1, this should fix the issue.
      该代码在使用Debian派生工具的RStudio上运行良好,但我认为,如果RStudio编辑器的默认编码为Windows-1252而不是UTF-8,则Windows计算机上的R可能存在一些编码问题。

      问候,
      马蒂亚斯

      马蒂亚斯 3年前 回复


  • Hellow, am currently working on Gumbel and weighted inverse Rayleigh distribution. so using flood 数据 i want to find out the return level through r for the given 数据. sir, please provide me r package 和code please
    问候,
    穆罕默德·伊什法克

    穆罕默德·伊什法克 2年前 回复


  • 嘿,非常感谢您的帖子,它非常有用!

    胡安·奥萨(Juan Ossa) 2年前 回复


  • 您好Matthias,
    我想要一个返回水平图,其中x轴旋转180°或类似的角度,以便我可以针对低放电事件进行绘制。较低的值应具有较高的回报水平。有没有办法做到这一点?

    凯瑟琳娜·穆勒(KatharinaMüller) 2年前 回复


    • 嗨Katharina,

      this use case prominently occurs is some research areas, for instance when working with droughts. You might want to have a look at the lfstatpackage in this context.

      马蒂亚斯 1年前 回复


  • 如何在R中用netcdf数据在全国30年,50年和100年的空间上绘制返回水平图

    帕万 5 月 ago 回复


发表回复

*