Wrong Hessian from optim in R

Question

I am doing some Extreme Values analysis. I don't want to use the fevd package for a variety of reasons (the first I want to be able to tweak some things that I cannot do otherwise). I wrote my own code. It is mostly very simple, and I thought I had solved everything. But for some parameter combinations, the Hessian coming out of my log-likelihood analysis (based on optim ) will not be correct.

Going over one step at the time. My code - or selected part of it - looks like this:

# routines for non stationary
Log_lik_GEV <- function(dataIN,scaleIN,shapeIN,locationIN){
    # simply calculate the negative log likelihood value for a set of X and parameters, for the GPD
    #xi, mu, sigma  - xi is the shape parameter, mu the location parameter, and sigma is the scale parameter.
    # shape = xi
    # location = mu
    # scale = beta
    library(fExtremes)
    #dgev   Density of the GEV Distribution, dgev(x, xi = 1, mu = 0, sigma = 1)

    LLvalues <- dgev(dataIN, xi = shapeIN, mu = locationIN, beta = scaleIN) 
    NLL <- -sum(log(LLvalues[is.finite(LLvalues)]))
    return(NLL)
}


function_MLE <- function(par , dataIN){
    scoreLL <- 0
    shape_param <- par[1]
    scale_param <- par[2]
    location_param <- par[3]
    scoreLL <- Log_lik_GEV(dataIN, scale_param, shape_param, location_param)
    if (abs(shape_param) > 0.3) scoreLL <- scoreLL*10000000
    if ((scale_param) <= 0) {
        scale_param <- abs(scale_param)
        par[2] <- abs(scale_param)
        scoreLL <- scoreLL*1000000000
    }
    sum(scoreLL) 
}

kernel_estimation <- function(dati_AM, shape_o, scale_o, location_o) {

     paramOUT <- optim(par = c(shape_o, scale_o, location_o), fn = function_MLE, dataIN = dati_AM, control = list(maxit = 3000, reltol = 0.00000001), hessian = TRUE)

     # calculation std errors
     covmat <- solve(paramOUT$hessian)
     stde <- sqrt(diag(covmat))
     print(covmat)

     print('')

     result <- list(shape_gev =paramOUT$par[1], scale_gev = paramOUT$par[2],location_gev =paramOUT$par[3], var_covar = covmat)

     return(result)
}

Everything works great, in some cases. If I run my routines and the fevd routines, I get exactly the same results. In some cases (in my specific case when shape=-0.29 so strongly negative/weibull), my routine will give negative variances and funky hessians. It is not always wrong, but some parameter combinations are clearly not giving valid hessian (Note: the parameters are still estimated correctly, meaning are identical to the fevd results, but the covariance matrix is completely off).

I found this post that compared the hessian from two procedures, and indeed optim seems to be flaky. However, if I simply substitute maxLik in my routine, it just doesn't converge at all (even in those cases when the convergence was happening).

 paramOUT = maxLik(function_MLE, start =c(shape_o, scale_o, location_o),   
                          dataIN=dati_AM, method ='NR' )

I tried to give different initial values - even the correct ones - but it just doesn't converge.

I am not supplying data because I think that the optim routine is used correctly in my example. Simply, the numerical results are not stable for some parameter combination. My question is:

1) Am I missing something in the way I use maxLik?

2) Are there other optimization routines, besides maxLik, from which I can extract the hessian?

thanks