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:


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 
Method:             Nelder-Mead 
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?

  • $\begingroup$ 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. $\endgroup$
    – Macro
    Aug 16 '13 at 17:27
  • $\begingroup$ @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.) $\endgroup$
    – Jen Bohold
    Aug 16 '13 at 17:33
  • $\begingroup$ 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. $\endgroup$
    – Macro
    Aug 16 '13 at 17:37
  • $\begingroup$ @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? $\endgroup$
    – Jen Bohold
    Aug 16 '13 at 17:39
  • 4
    $\begingroup$ 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? $\endgroup$
    – Jen Bohold
    Aug 16 '13 at 18:18

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
>>> sqrt(5.113433e-3)#for mu
>>> sqrt(1.5261031428)*0.002035#for delta
>>> sqrt(1.6617499980)*0.204827#for zeta

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.