A user recently reported a problem with estimation of the variances from the likelihood surface in evol.vcv(). This is done by first computing the Hessian matrix (the matrix of partial second derivatives) of the likelihood surface at the optimum, and then calculating the negative inverse of this matrix. This matrix is the (asymptotic) variance-covariance matrix of our ML parameter estimates. The element-wise square root of the diagonal should then provide the (asymptotic) standard errors.
The user reported problem is that in some rare cases an estimated variance obtained this way is negative. Although via troubleshooting I have managed to eliminate some possible reasons for this problem (for instance, failure to converge or errors in the calculation of the Hessian), I still have not managed to figure out why this is the case. I have not eliminated the possibility that this might have something to do with reparameterizing the likelihood optimization in terms of the Cholesky matrices. [Note that I then compute the Hessian from the likelihood function in terms of the VCV matrix.] Any suggestions (e.g., Carl?) are welcome!