The vcov function cannot be applied? I originally asked a question about the delta-method in the context of the hyperbolic distribution. I got an answer there, which is useful, except that it says I should apply the vcov function to my fitted function what is unfortunately not possible. The vcov was there to give me the variance-covariance matrix of the parameter $\zeta,\pi,\delta$. Now I want to calculate this variance-covariance manually. According to this question of me, I know how to get the variance-covariance matrix out of the hessian. Unfortunately according to this this is the hessian for pi, log(zeta), log(delta), and mu and not for the desired pi, zeta, delta, and mu.
So my basic problem is how to get the variance-covariance matrix of the parameter estimates of the hyperbolic distribution (in my original post this was the Sigma solved with the vcov function which does not work)? I know how to get the variance-covariance matrix out of the hessian, but I have the wrong hessian. So what should I do now? Is it possible to transform the hessian from pi, log(zeta), log(delta), and mu into the hessian of pi, zeta, delta, and mu and how is this done?
Or is there an easier way to get vcov running?
I am really stuck here and don't know how to continue.
My main question is: How to get the answer of COOLSerdash in this post working?
The vcov function is not working and also the hessian cannot be used, see my description of the error message and my comments to the answer.
 A: The reason that parameters passed to the optimizer are pi, log(zeta), log(delta), mu not zeta and delta is mostly likely to constrain the optimizer in R+ for zeta and delta. If you need the Hessian of pi, zeta, delta and mu rather than pi, log(zeta), log(delta), and mu, you can to write your own function, which is parametrized differently (in pi, zeta, delta and mu). Then calculate the 2nd derivatives of your new function at the MLE, which is given by thehyperbFit result. 
But why? Are you trying to get the confidence limit for zeta and delta?
A: Actually you already have the S.E. for zeta and delta, calculated from log(zeta) and log(delta), in summary(hyperbfitalv): Variance-covariance matrix of the parameter estimates wrongly calculated?
Let's look at it more closely and see how it was done: lzeta is log(zeta):
What you already have, from solve(hyperbfitalv$hessian), is:
$$
\sigma _{lzeta}^{2}
$$
Now you what to calculate:
$$
\sigma _{zeta}^{2}
$$
Notice: 
$$
zeta=e^{lzeta}
$$
So, applying Delta method:
$$
\sigma _{zeta}^{2}\approx (\frac{d\,zeta}{d\, lzeta})^2\cdot \sigma _{lzeta}^{2}=(\frac{d\,e^{lzeta}}{d\, lzeta})^2\cdot \sigma _{lzeta}^{2}=(e^{lzeta})^2\cdot\sigma _{lzeta}^{2}=zeta^2\cdot\sigma _{lzeta}^{2}
$$
which essentially is how the S.E. of zeta and delta were calculated in the 5th line of summary(hyperbfitalv) output in Variance-covariance matrix of the parameter estimates wrongly calculated?, such as:
>>> sqrt(1.5261031428)*0.002035#for delta
0.0025139483860139073

@COOLSerdash 's answer is great, but he probably didn't realize you already have 
$$
\sigma _{lzeta}^{2}
$$
