nonlinear optimization - MLE using nlminb in R - understand/debug certain errors

Question

This is my first question here, so I will try to make it as well written as possible. Please be overbearing should I make a silly mistake.

Briefly, I am trying to do a maximum likelihood estimation where I need to estimate 5 parameters. The general form of the problem I want to solve is as follows: A weighted average of three copulas, each with one parameter to be estimated, where the weights are nonnegative and sum to 1 and also need to be estimated.

There are packages in R for doing MLE on single copulas or on a weighted average of copulas with fixed weights. However, to the best of my knowledge, no packages exist to directly solve the problem I outlined above. Therefore I am trying to code the problem myself. There is one particular type of error I am having trouble tracing to its source. Below I have tried to give a minimal reproducible example where only one parameter needs to be estimated.

library(copula)

set.seed(150)
x <- rCopula(100, claytonCopula(250)) 

# Copula density
clayton_density <- function(x, theta){
  dCopula(x, claytonCopula(theta))
}

# Negative log-likelihood function
nll.clayton <- function(theta){
  theta_trans <- -1 + exp(theta) # admissible theta values for Clayton copula
  nll <- -sum(log(clayton_density(x, theta_trans)))
  return(nll)
}

# Initial guess for optimization
guess <- function(x){
  
  init <- rep(NA, 1)
  tau.n <- cor(x[,1], x[,2], method = "kendall")
  
  # Guess using method of moments 
  itau <- iTau(claytonCopula(), tau = tau.n)
  
  # In case itau is negative, we need a conditional statement
  # Use log because it is (almost) inverse of theta transformation above
  if (itau <= 0) {
    init[1] <- log(0.1) # Ensures positive initial guess
  }
  else {
    init[1] <- log(itau)
  }
}


estimate <- nlminb(guess(x), nll.clayton)
(parameter <- -1 + exp(estimate$par)) # Retrieve estimated parameter

fitCopula(claytonCopula(), x) # Compare with fitCopula function

This works great when simulating data with small values of the copula parameter, and gives almost exactly the same answer as fitCopula() every time.

For large values of the copula parameter, such as 250, the following error shows up when I run the line with nlminb():"Error in .local(u, copula, log, ...) : parameter is NA Called from: .local(u, copula, log, ...) Error during wrapup: unimplemented type (29) in 'eval'"

When I run fitCopula(), the optimization is finished, but this message pops up: "Warning message: In dlogcdtheta(copula, u) : dlogcdtheta() returned NaN in column(s) 1 for this explicit copula; falling back to numeric derivative for those columns"

I have been able to find out using debug() that somewhere in the optimization process of nlminb, the parameter of interest is assigned the value NaN, which then yields this error when dCopula() is called. However, I do not know at which iteration it happens, and what nlminb() is doing when it happens. I suspect that perhaps at some iteration, the objective function is evaluated at Inf/-Inf, but I do not know what nlminb() does next. Also, something similar seems to happen with fitCopula(), but the optimization is still carried out to the end, only with the abovementioned warning.

I would really appreciate any help in understanding what is going on, how I might debug it myself and/or how I can deal with the problem. As might be evident from the question, I do not have a strong background in coding. Thank you so much in advance to anyone that takes the time to consider this problem.

Update: When I run dCopula(x, claytonCopula(-1+exp(guess(x)))) or equivalently clayton_density(x, -1+exp(guess(x))), it becomes apparent that the density evaluates to 0 at several datapoints. Unfortunately, creating pseudobservations by using x <- pobs(x) does not solve the problem, which can be see by repeating dCopula(x, claytonCopula(-1+exp(guess(x)))). The result is that when applying the logarithm function, we get several -Inf evaluations, which of course implies that the whole negative log-likelihood function evaluates to Inf, as can be seen by running nll.clayton(guess(x)). Hence, in addition to the above queries, any tips on handling log(0) when doing MLE numerically is welcome and appreciated.

Second update Editing the second line in nll.clayton as follows seems to work okay:

nll <- -sum(log(clayton_density(x, theta_trans) + 1e-8))

However, I do not know if this is a "good" way to circumvent the problem, in the sense that it does not introduce potential for large errors (though it would surprise me if it did).

question from:https://stackoverflow.com/questions/65869827/mle-using-nlminb-in-r-understand-debug-certain-errors

Answer 1

Waitting for answers