Compare return levels of fitted GPD using MLE in different R packages

This question is related to this post: https://stackoverflow.com/questions/27524131/calculation-of-return-levels-based-on-a-gpd-in-different-r-packages

I want to constraint "potvalues" data to be in a period of 6 years, this is, 16 observations per year as the number of samples is 96. I want to do the calculation of return levels (mle parameters estimation and GP model) with extRemes package and compare the result with extremeStat package.

Be aware of the parameter truncate=0.4956645 to get exactly a threshold of 50.

Why result of line 51 (d$quant[28, ,drop=FALSE]) is not exactly equal to result of line 45 (rl.extremes2), if extremeStat package is using the same package extRemes with MLE and GP to do the calculation?

th <- 50

# sample data:
potvalues <- c(
  58.5,44.2,49.6,59.3,48.3,60.9,94.5,47.1,45.3,57.6,48.2,46.2,44.2,50.6,42.1,52.7,80.9,
  58.5,51.3,48.4,51.7,71.9,60.1,64.4,43.5,55.5,49.3,58.2,47.5,43.7,45.2,52.8,42.2,46.4,
  96.1,47.5,50.1,42.4,60.9,72.6,51.6,59.4,80.5,63.7,59.9,45.0,66.7,47.6,53.3,43.1,51.0,
  46.2,53.6,59.8,51.7,46.7,42.6,44.5,45.0,50.0,44.0,89.9,44.2,47.8,53.3,43.0,55.7,44.6,
  44.6,54.9,45.1,43.9,78.7,45.5,64.0,42.7,47.4,57.0,105.4,64.3,43.2,50.4,80.2,49.9,71.6,
  47.4,44.1,47.6,55.2,44.4,78.6,50.8,42.4,47.1,43.5,51.4)

#------------------------------------------------------------------------------------------#

#Count events over threshold
excesses = potvalues > th
sum(excesses)

# Data corresponding to a period of 6 years

#-------------------------------------------------------------------------------------------#
# MLE Fitting of GPD - package extRemes
# If fit period is 6 years, then I have 16 obs by year

pot.ext2 <- extRemes::fevd(potvalues, method = "MLE", type="GP", threshold=th, 
                           time.units="16/year")


npy2=16  #pot.ext2$npy
span2=5.9375 #pot.ext2$span

w2 = 96/npy2   #Duration of the fit period (6 years)
lambda2 = sum(excesses)/w2
Tr=c(2,5,10,20,50,100)
myp2 = (1 - (1/(lambda2*Tr)))
myp2 = myp2[myp2>0]
#Get return level using quantile function!
vel1 = extRemes::qevd(myp2, loc = pot.ext2$threshold, scale = pot.ext2$results$par[1], 
            shape = pot.ext2$results$par[2], 
            threshold = pot.ext2$threshold, type = "GP")
vel1
# return levels with 6 years, 16 obs, using return.level function
rl.extremes2 <-  extRemes::return.level(pot.ext2, conf = 0.05,
                             return.period= c(2,5,10,20,50,100))
rl.extremes2 <- as.numeric(rl.extremes2)
rl.extremes2
#------------------------------------------------------------------------------------------#
npy=16
Tr=c(2,5,10,20,50,100)
p = (1 - (1/(npy*Tr)))
d <- extremeStat::distLquantile(potvalues, truncate=0.4956645, probs=p, quiet=TRUE, list=TRUE)
d$quant[28, ,drop=FALSE]

dlf <- extremeStat::distLextreme(potvalues, quiet=TRUE, npy=16, truncate=0.4956645)
dlf$returnlev["threshold",1]
dlf$returnlev[28, , drop=FALSE]

Here the answer of the author of extremeStat

"due to time.units, your shape and scale of the estimated distribution function are different" ... "the probabilities are slightly different, but that doesn't affect the results (see outcommented test)"

Code below shows 1) how extremeStat is doing the calculations of return levels, using extRemes, and 2) how to work directly extRemes using two approaches: extRemes::qevd and extRemes::return.level

Maybe someone expert in the topic can take a decision of which one of the two ways is the 'more suitable', or 'correct' way to calculate the return levels: bb_extRemes or alexys_exRtemes?... the results are similar anyway.

# CLEAN ------------------------------------------------------------------------
rm(list=ls())

potvalues <- c(
  58.5,44.2,49.6,59.3,48.3,60.9,94.5,47.1,45.3,57.6,48.2,46.2,44.2,50.6,42.1,52.7,80.9,
  58.5,51.3,48.4,51.7,71.9,60.1,64.4,43.5,55.5,49.3,58.2,47.5,43.7,45.2,52.8,42.2,46.4,
  96.1,47.5,50.1,42.4,60.9,72.6,51.6,59.4,80.5,63.7,59.9,45.0,66.7,47.6,53.3,43.1,51.0,
  46.2,53.6,59.8,51.7,46.7,42.6,44.5,45.0,50.0,44.0,89.9,44.2,47.8,53.3,43.0,55.7,44.6,
  44.6,54.9,45.1,43.9,78.7,45.5,64.0,42.7,47.4,57.0,105.4,64.3,43.2,50.4,80.2,49.9,71.6,
  47.4,44.1,47.6,55.2,44.4,78.6,50.8,42.4,47.1,43.5,51.4)
options(scipen=10) # nicer printing of 0.00000001


bb_extRemes <- function(x, truncate, RPs=c(2,5,10,20,50), npy)
  {
  normalthr <- berryFunctions::quantileMean(x[is.finite(x)], truncate)
  z <- extRemes::fevd(x, method="MLE", type="GP", threshold=normalthr)
  scale <- z$results$par["scale"]
  shape <- z$results$par["shape"]
  probs <- 1-1/(RPs*npy)
  probs2 <- (probs-truncate)/(1-truncate) # correct probabilities for smaller sample proportion
  probs2[probs < truncate] <- 0   # avoid negative values
  probs2[probs2==0] <- NA
  probs2[probs2==1] <- NA
  #alexys
  #excesses <- x > normalthr
  #w2 <- length(x)/npy # Duration of the fit period (6 years)
  #lambda2 <- sum(excesses)/w2
  #probs <- 1 - 1/(lambda2*RPs)
  #probs[probs<0] <- NA  
  #alexys
  out <- extRemes::qevd(p=probs2, scale=scale, shape=shape, threshold=z$threshold, type="GP")
  names(out) <- paste0("RP.", RPs)
  out <- list(RL=out, PAR=c(thr=normalthr, scale=scale, shape=shape, probs=probs2))
  out
  }

alexys_exRtemes <- function(x, threshold, RPs=c(2,5,10,20,50), npy)
  {
  unit <- paste0(npy,"/year")
  z <- extRemes::fevd(x, method="MLE", type="GP", threshold=threshold, time.units=unit)
  excesses <- x > threshold
  w2 <- length(x)/npy # Duration of the fit period (6 years)
  lambda2 <- sum(excesses)/w2
  probs <- 1 - 1/(lambda2*RPs)
  probs[probs<0] <- NA
  ### probs     0.9375000         0.9750000         0.9875000         0.9937500         0.9975000 
  ###probs <- c(0.93803727875591, 0.97521491150236, 0.98760745575118, 0.99380372787559, 0.99752149115024)
  scale <- z$results$par[1]
  shape <- z$results$par[2]
  vel1 <- extRemes::qevd(probs, loc=z$threshold, scale=scale, shape=shape, threshold=z$threshold, type="GP")
  out <- extRemes::return.level(z, return.period=RPs)
  n <- names(out)
  out <- as.numeric(out)
  names(out) <- n
  out <- list(RL=out, PAR=c(thr=z$threshold, scale=scale, shape=shape, probs=probs, lambda2=lambda2))
  out
  }
bb_extRemes(potvalues, truncate=0.4956645, npy=16)
alexys_exRtemes(potvalues, threshold=50, npy=16)