Renext (version 3.0-0)

fGEV.MAX: Fit a GEV distribution from block maxima or r-largest order statistics using an aggregated Renewal POT process

Description

Fit a GEV distribution from block maxima or $r$-largest order statistics using an aggregated Renewal POT process.

Usage

fGEV.MAX(MAX.data, MAX.effDuration,
         MAX = NULL,
         control = list(maxit = 300, fnscale = -1),
         scaleData = TRUE,
         trace = 0)

Arguments

MAX.data

A list of block maxima or r-largest statistics as in Renouv.

MAX.effDuration

A vector of durations as in Renouv. The durations must be identical in order to have a common GEV distribution for the maxima.

MAX

A compact representation of the needed data as a list. This is typically created by using the (non exported) Renext:::maxkeMAXdata function.

control

List to control the optimisation in optim.

scaleData

Logical. If TRUE, the data in MAX.data are scaled before being used in the likelihood. The scaling operation is carried on the excesses (observations minus threshold).

trace

Integer level of verbosity during execution. With the value 0, nothing is printed.

Value

• A list
• estimate
  Named vector of the estimated parameters for the GEV distribution of maxima.
• opt
  Result of the optimisation.
• loglik
  Identical to opt$value. This is the maximised log-likelihood for the renewal POT model. It differs from the usual log-likelihood for GEV block maxima by a constant (w.r.t. parameters) depending on the number of blocks and the number of order statistics used.
• sd
  Standard deviation of the estimates (approximation based on the ML theory).
• cov
  Covariance matrix of the estimates (approximation based on the ML theory).

Details

The data are assumed to provide maxima or $r$-largest statistics arising from an aggregated renewal POT model with unknown event rate $\lambda$ and unknown two-parameter Generalised Pareto Distribution for the excesses. A threshold $u$ is fixed below the given data and the three unknown parameters lambda, scale and shape of the POT model are found by maximising the likelihood. Then a vector of the three parameters for the GEV distribution is computed by transformation. The covariance matrix and standard deviations are computed as well using the jacobian matrix of the transformation.

The maximisation if for the log-likelihood with the rate lambda concentrated out, so it is a two-parameters optimisation.

References

See the Renext Computing Details document.

See Also

Renouv.

Examples

##====================================================================
## block maxima: simulated data and comparison with  the 'fgev'
## function from the 'evd' package
##====================================================================
set.seed(1234)
u <- 10
nBlocks <- 30
nSim <- 100   ## number of samples 
Par <- array(NA, dim = c(nSim, 3, 2),
             dimnames = list(NULL, c("loc", "scale", "shape"), c("MAX", "evd")))
LL <- array(NA, dim = c(nSim, 2),
            dimnames = list(NULL, c("MAX", "evd")))

for (i in 1:nSim) {
  rd <- rRendata(threshold = u,
                 effDuration = 1,
                 lambda = 12,
                 MAX.effDuration = rep(1, nBlocks),
                 MAX.r = rep(1, nBlocks),
                 distname.y = "exp", par.y = c(rate = 1 / 100))

  MAX <- Renext:::makeMAXdata(rd)
  fit.MAX <- fGEV.MAX(MAX = MAX)
  fit.evd <- fgev(x = unlist(MAX$data))
  Par[i, , "MAX"] <- fit.MAX$estimate
  Par[i, , "evd"] <- fit.evd$estimate
  LL[i, "MAX"] <- fit.MAX$loglik
  LL[i, "evd"] <- logLik(fit.evd)
}

##====================================================================
## r-largest: use 'ismev::rlarg.fit' on the venice data set.
## NB 'venice' is taken from the 'evd' package here.
##====================================================================
require(ismev);
fit1 <- ismev::rlarg.fit(venice)

## transform data: each row is a block
MAX.data <- as.list(as.data.frame(t(venice)))
## remove the NA imposed by the rectangular matrix format
MAX.data <- lapply(MAX.data, function(x) x[!is.na(x)])
MAX.effDuration <- rep(1, length(MAX.data))

fit2 <- fGEV.MAX(MAX.data = MAX.data,
                 MAX.effDuration = MAX.effDuration)

## estimates
est <- cbind(ismev = fit1$mle, RenextLab = fit2$estimate)
print(est)
# covariance
covs <- array(dim = c(2, 3, 3),
              dimnames = list(c("ismev", "RenextLab"),
                colnames(fit2$cov), colnames(fit2$cov)))
                
covs["ismev", , ] <- fit1$cov
covs["RenextLab", , ] <- fit2$cov
print(covs)