In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) : NaNs produced

Question

I am estimating a log-likelihood using optim(). I am having some problems with the eigenvalues that don't let me find a valid hessian matrix and, therefore, the standard errors can't be calculated.

Here are the "warning" messages I get:

Warning messages:
1: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
2: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
3: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
4: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
5: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced

Can somebody help me find the problem?

Here's my code:

tobit.ll <- function(theta,x,y1,y2,y3)
{
library("copula")
library(mvtnorm)
    p <- nrow(theta)
    n<-nrow(x)
    x<-as.matrix(cbind( 1,x))
    km<- as.numeric(ncol(x))
    beta1 <- theta[1:km]
      beta2 <- theta[(km+1):(2*km)]
    beta3 <- theta[(2*km+1):(km*3)]
        s1  <- theta[(3*km+1)]   ### If sigma^2 is given, take sqrt()
            if (s1 < 0)   return(1e-10)
        s2  <- theta[(3*km+2)]   ### If sigma^2 is given, take sqrt()
            if (s2 < 0)   return(1e-10)
        s3  <- theta[(3*km+3)]  
            if (s3 < 0)   return(1e-10) 
        rho12   <- theta[(3*km+4)] 
            if (rho12< -1 || rho12 > 1)   return(NA)  
        rho13   <- theta[(3*km+5)] 
            if (rho13< -1 || rho13> 1)   return(NA)  
    rho23   <- theta[(3*km+6)]  
            if (rho23< -1 || rho23> 1)   return(NA) 
    fy1 <- ifelse(y1 >0, dnorm(y1,x%*%beta1, s1), (1-pnorm((x%*%beta1)/s1)))
    fy2 <- ifelse(y2 >0, dnorm(y2,x%*%beta2, s2), (1-pnorm((x%*%beta2)/s2)))
    fy3 <- ifelse(y3 >0, dnorm(y3,x%*%beta3, s3), (1-pnorm((x%*%beta3)/s3)))
    Fy1 <- ifelse(y1>0, pnorm((y1-x%*%beta1)/s1),(1-pnorm((x%*%beta1)/s1)))
    Fy2 <- ifelse(y2>0, pnorm((y2-x%*%beta2)/s2),(1-pnorm((x%*%beta2)/s2)))
    Fy3 <- ifelse(y3>0, pnorm((y3-x%*%beta3)/s3),(1-pnorm((x%*%beta3)/s3)))
    norm.cop= ellipCopula("normal", param = c(rho12, rho13, rho23),dim = 3, dispstr = "un")
    u <- as.matrix(cbind(Fy1,Fy2,Fy3))
    dcop<- dCopula(u, norm.cop)
    fy1.lg <- ifelse(fy1==0, log(0.0001), log(fy1))
    fy2.lg <- ifelse(fy2==0, log(0.0001), log(fy2))
    fy3.lg <- ifelse(fy3==0, log(0.0001), log(fy3))
    dcop.lg <- ifelse(dcop==0, log(0.0001), log(dcop)) 
    lglik <- sum(fy1.lg+fy2.lg+fy3.lg) +sum(dcop.lg)
    return(-(lglik))
        }

Answer 1

The parameters rho12, rho13, rho23 are apparently unconstrained: there is no guarantee that the correlation matrix is positive semi-definite.

You can force it to be positive semi-definite by "truncating" the negative eigenvalues (but that makes the function you are trying to minimize non-differentiable: this can hamper convergence).

# Sample data
rho12 <- .9
rho13 <- -.9
rho23 <- .9
correlation <- matrix( c( 1,rho12,rho13, rho12,1,rho23, rho13,rho23,1 ), 3, 3 )

# Fix the correlation matrix
e <- eigen(correlation)
e$values    # Not positive semi-definite
correlation <- e$vectors %*% diag(pmax(0,e$values)) %*% t(e$vectors)
correlation <- cov2cor(correlation)

Alternatively, you could parametrize the correlation matrix in a different way (e.g., using angles).