Estimating parameter covariance matrix when fitting a Generalized Extreme Value (GEV) model using 'extRemes' in R

Question

My question is similar to Error with fitting a Generalized Extreme Value (GEV) using `extRemes` in R?. However, I am fitting non-stationary Generalized Extreme Value (GEV) distribution, i.e., when the location parameter depends on covariates and is no longer a constant. I am also using extRemes package in R. The solution presented in the above post does not work in my case as the function (EnvStats::egevd) used in the solution does not work for non-stationary GEVs.

Since the origincal data is large, here is a subset of my data which gives me a warning that follows as a result of which it does not produce a matrix giving the parameter covariances.

library(extRemes)
y <-  c(4.844187, 4.844187, 4.174387, 4.744932, 4.094345, 4.174387, 4.343805, 4.624973,
  4.094345, 4.094345, 4.343805, 4.174387, 4.499810, 4.094345, 4.499810, 4.343805,
  4.653960, 4.499810, 4.941642, 4.624973, 4.499810, 4.744932, 4.317488, 4.499810,
  4.605170, 4.844187, 4.744932, 4.499810, 4.094345, 4.442651, 4.653960, 4.317488,
  4.744932, 4.174387, 4.605170, 4.700480, 4.744932, 4.553877, 4.248495, 4.094345,
  4.442651, 4.553877, 4.317488, 4.382027, 4.248495, 4.382027, 4.499810, 4.382027,
  4.382027, 4.382027, 4.248495, 4.248495, 4.248495, 4.382027, 4.248495, 4.744932,
  4.094345, 4.700480, 4.744932, 4.499810, 4.317488, 4.382027, 4.382027, 4.499810,
  4.248495, 4.499810, 4.382027, 4.700480, 4.653960, 4.744932, 4.605170, 4.094345,
  4.700480, 4.700480, 4.499810, 4.700480, 4.605170, 4.700480, 4.653960, 4.382027,
  4.499810, 4.174387, 4.174387, 4.174387, 4.499810, 4.382027, 4.553877, 4.499810,
  4.094345, 4.174387, 4.442651, 4.094345, 4.382027, 4.317488, 4.317488, 4.382027,
  4.442651, 4.442651, 4.174387, 4.094345)

x <-  c(4.535284, 4.556505, 3.305054, 4.293878, 3.228826, 3.375880, 3.562466, 4.047428,
   3.719651, 2.957511, 3.228826, 3.375880, 3.896909, 3.719651, 4.011868, 3.442019,
  4.208417, 3.812203, 4.617593, 4.011868, 3.766997, 4.397531, 3.504055, 4.114964,
  4.081766, 4.535284, 4.237723, 3.936716, 3.812203, 3.896909, 4.147095, 3.789855,
  4.293878, 3.305054, 4.147095, 4.147095, 4.293878, 3.896909, 3.644144, 3.562466,
  3.617652, 3.974998, 3.562466, 3.719651, 3.617652, 3.896909, 3.974998, 3.896909,
  4.147095, 3.766997, 3.695110, 3.695110, 3.695110, 3.896909, 3.504055, 4.480174,
  3.504055, 4.433789, 4.587515, 4.147095, 3.695110, 3.766997, 3.766997, 3.974998,
  3.695110, 4.081766, 3.896909, 4.421848, 4.347047, 4.587515, 4.421848, 3.504055,
  4.421848, 4.421848, 4.064744, 4.293878, 4.293878, 4.421848, 4.293878, 3.896909,
  4.193435, 3.766997, 3.504055, 3.644144, 4.064744, 3.974998, 4.147095, 4.064744,
  3.504055, 3.834061, 4.147095, 3.695110, 4.081766, 3.974998, 4.147095, 4.147095,
  4.266195, 4.223177, 3.812203, 3.504055)

fac <- as.factor(c(rep("a", 45), rep("b", 55)))

# Data 
dat <- data.frame(y = y, x= x, fac = fac)

# Non-stationary GEV model:
gev_mod <- fevd(dat$y, data = dat, location.fun = ~x*fac + I(x^2)*fac, 
    scale.fun=~1, shape.fun =~1, type = 'GEV')

Here is the output with the warning. Usually when you display the summary or the model output you get a matrix for covariances of parameters but in this case we don't. Although we get the parameter estimates but I need the parameter covariance matrix for inference purposes and I am not able to figure out where the problem lies. The default estimation method used is MLE, but I tried other methods like the Generalized MLE (GMLE) as well and still had the same issue. I am just starting out in the field and any help on how to fix this issue is much appreciated.

Warning messages:
1: In log(z) : NaNs produced
2: In log(z) : NaNs produced
3: In log(z) : NaNs produced
4: In log(z) : NaNs produced
5: In log(z) : NaNs produced
6: In log(z) : NaNs produced
7: In log(z) : NaNs produced
8: In log(z) : NaNs produced
> gev_mod

fevd(x = dat$y, data = dat, location.fun = ~x * fac + I(x^2) * 
    fac, scale.fun = ~1, shape.fun = ~1, type = "GEV")

[1] "Estimation Method used: MLE"


 Negative Log-Likelihood Value:  -114.3946 


 Estimated parameters:
       mu0        mu1        mu2        mu3        mu4        mu5      scale 
 5.5752723 -1.0508929 -2.3604995  0.1951670  1.0569155 -0.1229494  0.1034723 
     shape 
-0.7972565 

 AIC = -212.7893 

 BIC = -191.9479 
> summary(gev_mod$cov.theta)
Length  Class   Mode 
     0   NULL   NULL 

Answer 1

I was able to figure out a way to solve the problem. I looked at the source code for the fevd function in the extRemes package and noticed that log(z) appeared in the grlevd function which is used to estimate the gradient of likelihood needed for estimation. I used the trace() function in R to change the log(z) to log1p(z) and that seemed to have worked. Now I don't get any warnings and also get the estimated parameter covariance matrix. Just a side note, having worked with EVDs over the past couple of days, it looks like convergence in optimization could be tricky, especially for the non-stationary models.