extreme: New package for extreme value theory

[David Roodman]

Kit Baum has just posted a new package for fitting extreme value theory (EVT) models in Stata.

My impression is that the excellent textbook by Stuart Coles and the associated ismev package for R dominate most people's understanding of EVT. In addition there are many other advanced R packages for EVT. The new package extreme brings much of ismev's functionality to Stata. Using Maximum Likelihood, it fits the generalized Pareto distribution (GPD) and the generalized extreme value distribution (GEV), including the extension for multiple order statistics (such as the top five daily rainfall values for each year). All model parameters can have covariates, taking variables in the data set as linear determinants; in other words, models can be non-stationary. The package also offers various diagnostic plots. Clickable examples in the help file replicate most of the statistical and graphical results in the Coles book.

The major novelty is an emphasis on small-sample bias correction. (Maximum Likelihood is in general biased, and the bias can be significant in the small samples often used in EVT.) extreme offers analytical bias correction according to the Cox-Snell formula (roughly speaking). The Cox-Snell correction for the GPD is derived in a new paper by Giles, Feng, and Godwin. That for the GEV, including for multiple order stats, is new, and I'm writing it up now. Simulations confirm the efficacy of these corrections. The package also offers (parametric) bootstrap-based bias correction, which is slower but, unlike the Cox-Snell approach, can work well even for substantially negative values of the shape parameter (often labelled xi).

I think the influence of R is interesting. Like many R packages, a single command offers estimation and associated plots. There's also a Mata interface, though I haven't documented it yet, giving access to the estimation functionality from a proper programming language, as in R. This is made possible by Stata's Mata interface for maximum likelihood estimation.

I think extreme functionally subsumes the packages gevfit and gumbelfit.

Install with "ssc install extreme". Comments welcome.

--David

[Nick Cox]

Sounds great. This is just to note that an existing package on SSC, extremes, has quite different aims and content. Let's hope people can distinguish.

[David Roodman]

I have just posted, through my blog some simulation results for the program and its various options for bias correction:

http://davidroodman.com/blog/2015/01...distributions/

I have to say, despite all the effort I put into deriving the analytical bias correction for the generalized extreme value distribution, bootstrapped-based correction is quite competitive.

[Volker Seiler]

Dear David,

Great package, thank you very much. I was wondering whether there is a multivariate version/extension available?

Kind regards,
Volker

[David Roodman]

Sorry, Volker. There isn't. You might check https://www.ral.ucar.edu/~ericg/softextreme.php.

[Didar Khan]

David, please could you guide me how to perform this manually or you could share the programmings, fro estimation of GP using MLE on Stata ([email protected])

[David Roodman]

All the code is included in the download package, though it is not easy reading. https://ideas.repec.org/c/boc/bocode/s457953.html.

[David Roodman]

This is the key section in extreme.mata for computing the likelihoods. For GP, GEV=0.

        GEV =  moptimize_util_userinfo(M, 1)
 
        sig = exp(moptimize_util_xb(M, b, 1+GEV))
        xi  =     moptimize_util_xb(M, b, 2+GEV)
        if (GEV) {
               mu = moptimize_util_xb(M, b, 1)
               x = moptimize_util_userinfo(M, 2) :- mu
               R = cols(x)
               x = x :/ sig
               r =  moptimize_util_userinfo(M, 3)
        } else {
               x = moptimize_util_userinfo(M, 2) :/ sig
                r = R = 1
        }
        gumbel  = rows(xi)>1? _selectindex(abs(xi):<1e-10) : (abs(xi):<1e-10? 1::rows(x) : J(0,0,0) )
        if (rows(gumbel)) xgumbel=x[gumbel,]
        inv_xi = 1:/xi
        lny = ln(y = 1:+xi:*x)
 
        if (GEV) {
               f = (R>1? y[,R] : y):^-inv_xi
               if (rows(gumbel)) f[gumbel] = exp(-(R>1? xgumbel[,R] : xgumbel))
        }
... 
         lng = quadrowsum((inv_xi:+1):*lny, 1)  
         if (rows(gumbel)) lng[gumbel] = quadrowsum(xgumbel)  
         lng = (R>1? -r*ln(sig) : -ln(sig)) :- lng