GevModelling            package:fExtremes            R Documentation

_G_e_n_e_r_a_l_i_z_e_d _E_x_t_r_e_m_e _V_a_l_u_e _M_o_d_e_l_l_i_n_g

_D_e_s_c_r_i_p_t_i_o_n:

     A collection and description functions to compute the generalized
     extreme value distribution and to estimate it parameters. The
     functions compute  density, distribution function, quantile
     function  and generate random deviates for the GEV, for the 
     Frechet, Gumbel, and Weibull distributions. To model the GEV three
     types of approaches for parameter  estimation are provided:
     Maximum likelihood estimation, probability weighted moment method,
     and estimation by the MDA approach. MDA includes functions for the
     Pickands, Einmal-Decker-deHaan,  and Hill estimators together with
     several plot  variants. 

     The GEV distribution functions are:

       'dgev'           density of the GEV Distribution,
       'pgev'           probability function of the GEV Distribution,
       'qgev'           quantile function of the GEV Distribution,
       'rgev'           random variates from the GEV Distribution.
       '[dpqr]frechet'  Frechet Distribution,
       '[dpqr]gumbel'   Gumbel Distribution,
       '[dpqr]weibull'  Weibull Distribution,
       '[dpqr]evd'      an alternative call for the GEV Distribution.

     The GEV modelling functions are:

       'gevSim'         generates data from the GEV,
       'gevFit'         fits empirical or simulated data to the distribution,
       'print'          print method for a fitted GEV object,
       'plot'           plot method for a fitted GEV object,
       'summary'        summary method for a fitted GEV object,
       'gevrlevelPlot'  k-block return level with confidence intervals.

     Maximum Domain of Attraction Estimators:

       'hillPlot'     shape parameter and Hill estimate of the tail index,
       'shaparmPlot'  variation of shape parameter with tail depth.

_U_s_a_g_e:

     dgev(x, xi = 1, mu = 0, sigma = 1, log = FALSE)
     pgev(q, xi = 1, mu = 0, sigma = 1, lower.tail = TRUE)
     qgev(p, xi = 1, mu = 0, sigma = 1, lower.tail = TRUE)
     rgev(n, xi = 1, mu = 0, sigma = 1)
     devd(x, loc = 0, scale = 1, shape = 0, log = FALSE)
     pevd(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)
     qevd(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)
     revd(n, loc = 0, scale = 1, shape = 0)

     gevSim(model = list(shape = 0.25, location = 0, scale = 1), n = 1000)
     gevFit(x, block = NA, type = c("mle", "pwm"), gumbel = FALSE, ...)
     ## S3 method for class 'gevFit':
     print(x, ...)
     ## S3 method for class 'gevFit':
     plot(x, which = "all", ...)
     ## S3 method for class 'gevFit':
     summary(object, doplot = TRUE, which = "all", ...)
     gevrlevelPlot(object, k.blocks = 20, add = FALSE, ...)

     hillPlot(x, option = c("alpha", "xi", "quantile"), start = 15,
         end = NA, reverse = FALSE, p = NA, ci = 0.95, autoscale = TRUE, 
         labels = TRUE, ...) 
     shaparmPlot(x, revert = FALSE, standardize = FALSE, tails = 0.01*(1:10), 
         doplot = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE), 
         which = c(TRUE, TRUE, TRUE), doprint = TRUE, both.tails = TRUE, 
         xi.range = c(-0.5, 1.5), alpha.range = c(0, 10)) 

_A_r_g_u_m_e_n_t_s:

     add: [gevrlevelPlot] - 
           whether the return level should be added graphically to a
          time series plot; if 'FALSE' a graph of the profile
          likelihood curve showing the return level and its confidence
          interval is produced. 

alpha.range, xi.range: [saparmPlot] - 
           plotting ranges. 

autoscale: [hillPlot] - 
           whether or not plot should be automatically scaled; if not,
          'xlim' and 'ylim' graphical parameters may be entered. 

   block: [gevFit] - 
           the block size. Only used if 'method="mle"' is selected. A
          numeric value is interpreted as the number of data values in
          each successive block. All the data is used, so the last
          block may not contain 'block' observations. If the 'data' has
          a 'times' attribute containing (in an object of class
          '"POSIXct"', or an object that can be converted to that
          class; see 'as.POSIXct') the times/dates of each observation,
          then 'block' may instead take the character values '"month"',
          '"quarter"', '"semester"' or '"year"'. 

both.tails: [shaparmPlot] - 
           a logical, decides whether or not both tails should be 
          investigated. By default TRUE. If FALSE only the lower tail
          will be investigated. 

      ci: [hillPlot] - 
           probability for asymptotic confidence band; for no
          confidence band set 'ci' to zero. 

  doplot: a logical. Should the results be plotted? 
           [shaparmPlot] - 
           a vector of logicals of the same lengths as tails defining
          for wich tail depths plots should be created, by default
          plots will be generated for a tail depth of 5 percent. By
          default 'c(FALSE, FALSE, FALSE, FALSE,  TRUE, FALSE, FALSE,
          FALSE, FALSE, FALSE)'. 

 doprint: [shaparmPlot] - 
           a logical, decides whether or not for all tail depths the
          result for the shape parameter '1/alpha' should be printed. 

  gumbel: [gevFit] - 
           a logical, by default FALSE. To fit a Gumbel model with
          fixed 'shape=0' set 'gumbel=TRUE'. 

k.blocks: [gevrlevelPlot] - 
           specifies the particular return level to be estimated;
          default  set arbitrarily to 20. 

  labels: [hillPlot] - 
           whether or not axes should be labelled. 

loc, scale, shape: 'loc' is the location parameter,  'scale' the scale
          parameter, and 'shape' is the shape parameter. The default
          values are 'loc=0', 'scale=1', and  'shape=0'. 

     log: a logical, if 'TRUE', the log density is returned. 

lower.tail: a logical, if 'TRUE', the default, then probabilities are
          'P[X <= x]', otherwise, 'P[X > x]'. 

   model: [gevSim] - 
           a list with components 'shape', 'location' and  'scale'
          giving the parameters of the GEV distribution. By default the
          shape parameter has the value 0.25, the location is zero and
          the scale is one. To fit random deviates from a Gumbel
          distribution set  'shape=0'. 

       n: [gevSim] - 
           number of generated data points, an integer value. 
           [rgev][revd] - 
                   the number of observations. 

  object: [summary][grlevelPlot] - 
           a fitted object of class '"gevFit"'. 

  option: [hillPlot] - 
           whether 'alpha', 'xi' (1/alpha) or 'quantile' (a quantile
          estimate) should be plotted. 

       p: [qgev][qevs] - 
           a numeric vector of probabilities. [hillPlot] - 
           probability required when option 'quantile' is chosen. 

       q: [pgev][pevs] - 
           a numeric vector of quantiles. 

 reverse: [hillPlot] - 
           whether plot is to be by increasing threshold, 'TRUE', or 
          increasing number of order statistics 'FALSE'. 

  revert: [shaparmPlot] - 
           a logical value, by default FALSE, if set to TRUE the sign
          of the vector will be reverted: 'x = -x'. 

start, end: [hillPlot] - 
           lowest and highest number of order statistics at which to
          plot a point. 

standardize: [shaparmPlot] - 
           a logical value, by default FALSE, if set to  TRUE the
          vector 'x' will be standardized:  'x =
          (x-mean(x))/sqrt(var(x))'. 

   tails: [shaparmPlot] - 
           a numeric vector of tail depths to be considered; by default
          ten values ranging from 0.1 to 1.0 in steps of 0.1
          corresponding to values ranging from 1 to 10 percent. 

    type: a character string denoting the type of parameter estimation,
          either by maximum likelihood estimation '"mle"', the default
          value, or by the probability weighted moment menthod '"pwm"'. 

   which: [shaparmPlot] - 
           a vector of 3 logicals indicating which plots from the 
          three methods will be created. The first entry decides for
          the Pickands estimator, the second for the Hill  estimator,
          and the last for the Deckers-Einmahl-deHaan estimator. By
          default all three will be created. By default 'c(TRUE, TRUE,
          TRUE)'. 
           [plot][summary] - 
           a vector of logicals, one for each plot, denoting which plot
           should be displayed. Alkternatively if 'which="ask"' the
          user will be interactively asked which of the plots should be
          desplayed. By default 'which="all"'. 

       x: [dgev][devd] - 
           a numeric vector of quantiles. 
           [gevFit] - 
           data vector. In the case of 'method="mle"' the
          interpretation  depends on the value of block: if no block
          size is specified then  data are interpreted as block maxima;
          if block size is set, then data are interpreted as raw data
          and block maxima are calculated.  
           [hillPlot][shaparmPlot] - 
           the data from which to calculate the shape parameter, a
          numeric vector. 
           [print][plot] - 
           a fitted object of class '"gevFit"'. 

xi, mu, sigma: 'xi' is the shape parameter,  'mu' the location
          parameter, and 'sigma' is the scale parameter. The default
          values are 'xi=1', 'mu=0', and  'sigma=1'. 

     ...: [gevFit] - 
           control parameters optionally passed to the  optimization
          function. Parameters for the optimization function are passed
          to components of the 'control' argument of 'optim'.  
           [hillPlot] - 
           other graphics parameters. 
           [plot][summary] - 
           arguments passed to the plot function. 

_D_e_t_a_i_l_s:

     *Generalized Extreme Value Distribution:* 

      Computes density, distribution function, quantile function and 
     generates random variates for the Generalized Extreme Value 
     Distribution, GEV, for the Frechet, Gumbel, and Weibull 
     distributions. 

     *Parameter Estimation:* 

      'gevFit' estimates the parameters either by the probability 
     weighted moment method, 'method="pwm"' or by maximum log 
     likelihood estimation 'method="mle"'.  As a limiting case the
     Gumbel distribution can be selected. The  summary method produces
     diagnostic plots for fitted GEV or Gumbel  models. 

     *Methods:* 

      'print.gev', 'plot.gev' and 'summary.gev' are print, plot, and
     summary methods for a fitted object of class  'gev'. Concerning
     the summary method, the data are  converted to unit exponentially
     distributed residuals under null  hypothesis that GEV fits. Two
     diagnostics for iid exponential data  are offered. The plot method
     provides two different residual plots  for assessing the fitted
     GEV model. Two diagnostics for  iid exponential data are offered.  

     *Return Level Plot:* 

      'gevrlevelPlot' calculates and plots the k-block return level 
     and 95% confidence interval based on a GEV model for block maxima,
      where 'k' is specified by the user. The k-block return level  is
     that level exceeded once every 'k' blocks, on average. The  GEV
     likelihood is reparameterized in terms of the unknown return 
     level and profile likelihood arguments are used to construct a 
     confidence interval.  

     *Hill Plot:* 

      The function 'hillPlot' investigates the shape parameter and 
     plots the Hill estimate of the tail index of heavy-tailed data, or
      of an associated quantile estimate. This plot is usually
     calculated  from the alpha perspective. For a generalized Pareto
     analysis of  heavy-tailed data using the 'gpdFit' function, it
     helps to  plot the Hill estimates for 'xi'.  

     *Shape Parameter Plot:* 

      The function 'shaparmPlot' investigates the shape parameter and 
     plots for the upper and lower tails the shape parameter as a
     function  of the taildepth. Three approaches are considered, the
     _Pickands_  estimator, the _Hill_ estimator, and the
     _Decker-Einmal-deHaan_ estimator.

_V_a_l_u_e:

     'd*' returns the density, 
      'p*' returns the probability, 
      'q*' returns the quantiles, and 
      'r*' generates random variates. 
      All values are numeric vectors. 

     'gevSim' 
      returns a vector of data points from the simulated series. 

     'gevFit' 
      returns an object of class 'gev' describing the fit. 

     'print.summary' 
      prints a report of the parameter fit. 

     'summary' 
      performs diagnostic analysis. The method provides two different 
     residual plots for assessing the fitted GEV model.   

     'gevrlevelPlot' 
      returns a vector containing the lower 95% bound of the confidence
      interval, the estimated return level and the upper 95% bound.  

     'hillPlot' 
      displays a plot. 

     'shaparmPlot'  
      returns a list with one or two entries, depending on the
     selection of the input variable 'both.tails'. The two  entries
     'upper' and 'lower' determine the position of  the tail. Each of
     the two variables is again a list with entries  'pickands',
     'hill', and 'dehaan'. If one of the  three methods will be
     discarded the printout will display zeroes.

_N_o_t_e:

     *Generalized Extreme Value Distribution:* 

      Here we implemented the notation for the arguments as used by the
     GEV functions in the EVIS package or SPlus/FinMetrics  module.
     Additional arguments to these packages are the 'log'  and the
     'lower.tail' arguments, underlying the code from R's 'evd'
     package. 
      An alternative usage is proposed by the 'evd' package. There the
     following arguments are used: 
      '*gev(x, loc = 0, scale = 1, shape = 0, ...)' 
         What you prefer is a matter of taste. The GEV functions from
     the 'evd' package are renamed from '*gev' to '*evd' so that both
     versions are available.  
         In addition functions for the density, probability, quantiles,
      and the generation of random variates for the  Frechet
     '[dpqr]frechet', Gumbel '[dpqr]gumbel', and Weibull
     '[dpqr]weibull' are also available. 
         If you stay with both versions you can access the help page
     for  'evds''s function 'dgev' by 'help(dgev, package="evd")'. 

     *Generalized Extreme Value Distribution:* 

      If method '"mle"' is selected the parameter fitting in 'gevFit' 
     is passed to the internal function 'gev.mle' or 'gumbel.mle'
     depending on the value of 'gumbel', 'FALSE' or 'TRUE'. On the
     other hand, if method '"pwm"' is selected the parameter  fitting
     in 'gevFit' is passed to the internal function  'gev.pwm' or
     'gumbel.pwm' again depending on the value of  'gumbel', 'FALSE' or
     'TRUE'.

_A_u_t_h_o_r(_s):

     Alec Stephenson for R's 'evd' and 'evir' package, and 
      Diethelm Wuertz for this R-port.

_R_e_f_e_r_e_n_c_e_s:

     Coles S. (2001); _Introduction to Statistical Modelling of Extreme
     Values_, Springer.

     Embrechts, P., Klueppelberg, C., Mikosch, T. (1997); _Modelling
     Extremal Events_,  Springer.

_E_x_a_m_p_l_e_s:

     ## SOURCE("fExtremes.52A-GevModelling")

     ## *gev  -
        xmpExtremes("\nStart: GEV Frechet >")
        # Create and plot 1000 GEV/Frechet distributed rdv:
        par(mfrow = c(3, 3))
        r = rgev(n = 1000, xi = 1)
        plot(r, type = "l", main = "GEV/Frechet Series")
        ## Plot empirical density and compare with true density:
        ## Omit values greater than 500 from plot
        hist(r[r<10], n = 25, probability = TRUE, xlab = "r", 
          xlim = c(-5, 5), ylim = c(0, 1.1), main = "Density")
        x = seq(-5, 5, by=0.01)
        lines(x, dgev(x, xi = 1), col = 2)
        ## Plot df and compare with true df:
        plot(sort(r), (1:length(r)/length(r)), 
          xlim = c(-3, 6), ylim = c(0, 1.1),
          cex = 0.5, ylab = "p", xlab = "q", main = "Probability")
        q = seq(-5,5, by=0.1)
        lines(q, pgev(q, xi=1), col=2)
     ## Compute quantiles, a test:
        qgev(pgev(seq(-5, 5, 0.25), xi = 1), xi = 1) 
        
     ## *gev -
        xmpExtremes("\nNext: GEV Gumbel >")
        # Create and plot 1000 Gumbel distributed rdv:
        ##> r = rgev(n = 1000, xi = 0)
        ##> plot(r, type = "l", main = "Gumbel Series")
        ## Plot empirical density and compare with true density:
        ##>hist(r[abs(r)<10], nclass = 25, freq = FALSE, xlab = "r", 
        ##>   xlim = c(-5,5), ylim = c(0,1.1), main = "Density")
        ##>x = seq(-5, 5, by = 0.01)
        ##>lines(x, dgev(x, xi = 0), col=2)
        ## Plot df and compare with true df:
        ##>plot(sort(r), (1:length(r)/length(r)), 
        ##>  xlim = c(-3, 6), ylim = c(0, 1.1),
        ##>   cex=0.5, ylab = "p", xlab="q", main="Probability")
        ##>q = seq(-5, 5, by = 0.1)
        ##>lines(q, pgev(q, xi = 0), col = 2)
        ## Compute quantiles, a test:
        ##>qgev(pgev(seq(-5, 5, 0.25), xi = 0), xi = 0)   

     ## *gev -
        xmpExtremes("\nNext: GEV Weibull >")
        # Create and plot 1000 Weibull distributed rdv:
        r = rgev(n = 1000, xi = -1)
        plot(r, type = "l", main = "Weibull Series")
        ## Plot empirical density and compare with true density:
        hist(r[abs(r)<10], nclass = 25, freq = FALSE, xlab = "r", 
          xlim=c(-5,5), ylim=c(0,1.1), main="Density")
        x = seq(-5, 5, by=0.01)
        lines(x, dgev(x, xi = -1), col = 2)
        ## Plot df and compare with true df:
        plot(sort(r), (1:length(r)/length(r)), 
          xlim = c(-3, 6), ylim = c(0, 1.1),
          cex = 0.5, ylab = "p", xlab = "q", main = "Probability")
        q=seq(-5, 5, by = 0.1)
        lines(q, pgev(q, xi = -1), col = 2)
        ## Compute quantiles, a test:
        qgev(pgev(seq(-5, 5, 0.25), xi = -1), xi = -1)   

     ## gevSim -
     ## gevFit -
        # Simulate GEV Data:
        xmpExtremes("\nStart: Simulte GEV Sample >")
        # Use default length n=1000
        x = gevSim(model = list(shape = 0.25, location =0 , scale = 1))
        # Fit GEV Data by Probability Weighted Moments:
        fit = gevFit(x, type = "pwm") 
        print(fit)
        # Summarize Results:
        par(mfcol = c(3, 2))
        summary(fit)
        
     ## gevFit -
        # Fit GEV Data by Max Log Likelihood Method:
        xmpExtremes("\nNext: Estimate Parameters >")
        fit = gevFit(x, type = "mle") 
        print(fit)
        # Summarize Results:
        summary(fit) 
         
     ## gevSim -
     ## gevFit -
        # Simulate Gumbel Data:
        xmpExtremes("\nNext: Simulte Gumbel Sample >")
        # Use default length n=1000
        ##> x = gevSim(model = list(shape = 0, location = 0, scale = 1))
        # Fit Gumbel Data by Probability Weighted Moments:
        ##> fit = gevFit(x, type = "pwm", gumbel = TRUE) 
        ##> print(fit)
        # Summarize Results:
        ##> par(mfcol = c(3, 2))
        ##> summary(fit)  
        
     ## Fit Gumbel Data by Max Log Likelihood Method:
        xmpExtremes("\nNext: Estimate Parameters >")
        ##> fit = gevFit(x, type = "mle", gumbel = TRUE) 
        ##> print(fit)
        # Summarize Results:
        ##> summary(fit)     
        ##> xmpExtremes("Press any key to continue >") 
        
     ## Return levels based on GEV Fit:
        # BMW Stock Data:
        xmpExtremes("\nNext: Compute BMW Return Levels >")
        par(mfrow = c(2, 1))
        data(bmw)
        # Fit GEV to monthly Block Maxima:
        fit = gevFit(-bmw, block = "month") 
        # Calculate the 40 month return level
        gevrlevelPlot(fit, k.block = 40, main = "BMW: Return Levels")   
        
     ## Return levels based on GEV Fit:
        xmpExtremes("\nNext: Compute SIEMENS Return Levels >")
        # Siemens Stock Data:
        data(siemens)
        # Fit GEV to monthly Block Maxima:
        fit = gevFit(-siemens, block = "month") 
        # Calculate the 40 month return level
        gevrlevelPlot(fit, k.block = 40, main = "SIEMENS: Return Levels")
        
     ## Interactive Plot:
        ##> par(mfrow = c(1, 1), ask = TRUE)
        ##> plot(fit)
        
     ## hillPlot -
        xmpExtremes("\nStart: Hill Estimator >")
        # Hill plot of heavy-tailed Danish fire insurance data 
        # and BMW stock data for estimated 0.999 quantile
        par(mfrow = c(2, 2))
        data(bmw)
        hillPlot(bmw)
        hillPlot(bmw, option = "quantile", end = 500, p = 0.999)
        data(danish)
        hillPlot(danish)
        hillPlot(danish, option = "quantile", end = 500, p = 0.999)
        
     ## shaparmPlot -
        xmpExtremes("\nNext: Shape Parameter Plots >")
        par(mfcol = c(3, 2), cex = 0.6)
        data(bmw)
        shaparmPlot(bmw)

     ## shaparmPlot -
        xmpExtremes("\nNext: Simulated Frechet Data >")
        par(mfcol = c(3, 2), cex = 0.6)
        set.seed(4711)
        x = rgev(10000, xi = 1/4)
        shaparmPlot(x, revert = TRUE, both.tails = FALSE)
        lines(c(0.01, 0.1), c(4, 4), col = "steelblue3") # True Value

