Variance-covariance matrix of the parameter estimates wrongly calculated?

I fitted an hyperbolic distribution to my data with the hyperbFit(mydata,hessian=TRUE) command (package HyperbolicDist). The hessian looks like:

> hyperbfitmymodel$hessian hyperbPi lZeta lDelta mu hyperbPi 536.61654 -23.82800 25.62345 26153.16561 lZeta -23.82800 250.74196 -261.20570 -35.58481 lDelta 25.62345 -261.20570 272.77771 182.75927 mu 26153.16561 -35.58481 182.75927 2028904.75586 Now I want to calculate the variance-covariance matrix of the parameter estimates, according to this page 2: The asymptotic covariance matrix of$\hat{\theta}$is given by the inverse of the negative of the Hessian matrix evaluated at$\hat{\theta}$. I therefore calculate: solve(-hyperbfitalv$hessian)

which gives

hyperbPi         lZeta        lDelta            mu
hyperbPi -5.113433e-03 -0.0091511819 -0.0083271877  6.650321e-05
lZeta    -9.151182e-03 -1.6617499980 -1.5905496996  2.320893e-04
lDelta   -8.327188e-03 -1.5905496996 -1.5261031428  2.169113e-04
mu        6.650321e-05  0.0002320893  0.0002169113 -1.365591e-06

This looks clearly wrong to me, because there are negative values for the variance, but a variance cannot be negative? The covariance yes, but not the variance?

EDIT: The complete output of hyperbFit(mydata,hessian=TRUE):

Data:     mydata
Parameter estimates:
pi       zeta      delta         mu
0.090747   0.204827   0.002035  -0.002494
Likelihood:         756.911
Convergence code:   0
Iterations:         365

2nd EDIT: If I use solve(hyperbfitalv$hessian) I get hyperbPi lZeta lDelta mu hyperbPi 5.113433e-03 0.0091511819 0.0083271877 -6.650321e-05 lZeta 9.151182e-03 1.6617499980 1.5905496996 -2.320893e-04 lDelta 8.327188e-03 1.5905496996 1.5261031428 -2.169113e-04 mu -6.650321e-05 -0.0002320893 -0.0002169113 1.365591e-06 3rd EDIT: The output of summary(hyperbfitalv): Data: mydata Parameter estimates: pi zeta delta mu 0.090747 0.204827 0.002035 -0.002494 ( 0.071508) ( 0.264040) ( 0.002514) ( 0.001169) Likelihood: 756.911 Method: Nelder-Mead Convergence code: 0 Iterations: 365 4th EDIT: Ok, this is the hessian of pi, log(zeta), log(delta), and mu but how can I get the hessian of pi, zeta, delta and mu? • This hessian matrix you've shown is positive definite. Therefore, either the optimization algorithm did not converge or the converged point is a local minimum, not a maximum. I suspect this is a computing error. I suggest supplying the rest of the code you've entered off screen and posting this on stack overflow. – Macro Aug 16 '13 at 17:27 • @Macro Thanks for the answer, but there is no other code? I just used hyperbFit(mydata,hessian=TRUE) and that's it? Before I load the package HyperbolicDist. But nothing more? I added the complete output, not only the hessian. As you can see the convergence code is 0. In the manual of the command "optim" it says: "0 indicates successful completion". (I used default values so in default it is Nelder-Mead of the optim command.) – Jen Bohold Aug 16 '13 at 17:33 • Hi @Jen - there must be some information missing because, according to the help file, hyperbFit doesn't even take an argument called hessian. I get a warning when I try to pass that argument. You calculate the hessian by an appropriate call to the summary function, e.g. summary(fittedmodel, hessian=TRUE). You may try refitting the model with different starting values because I don't think the model converged. Again, I think this is a computing issue so I suggest deleting this and re-posting on stackoverflow. – Macro Aug 16 '13 at 17:37 • @Macro help.rmetrics.org/HyperbolicDist/hyperbFit.html of course there is an hessian argument? Your link does not work for me. I did not use a summary function up to now? I just used the code as it is and it works with the commands and especially with the hessian? – Jen Bohold Aug 16 '13 at 17:39 • OK, the manual says: It is the hessian of " pi, log(zeta), log(delta), and mu", but how can I get the hessian of pi,zeta,delta and mu? So without the logs? – Jen Bohold Aug 16 '13 at 18:18 2 Answers Optim in default setting is doing minimization, see the manual: By default optim performs minimization So the output is already the negative hessian. It should be further noted that: Because the parameters in the call to the optimiser are pi, log(zeta), log(delta), and mu, the delta method is used to obtain the standard errors for zeta and delta. Source here. Seconded to @Jen 's answer. In fact the 5th line in the result of summary(hyperbfitalv) are SE's. They are indeed the square root of the diagonal elements of inverse-hessian solve(hyperbfitalv$hessian).

>>> sqrt(1.365591e-6)#for pi
0.0011685850418347824
>>> sqrt(5.113433e-3)#for mu
0.071508272248740568
>>> sqrt(1.5261031428)*0.002035#for delta
0.0025139483860139073
>>> sqrt(1.6617499980)*0.204827#for zeta
0.26404019669949413

Note that lZeta and lDelta are in fact log(Zeta) and log(Delta). Cheers!

• I think you switched up pi and mu? – Jen Bohold Aug 16 '13 at 18:45