Standard errors of hyperbFit?

Question

I want to fit a hyperbolid distribution according to my notation:

\begin{align*} H(l;\alpha,\beta,\mu,\delta)&=\frac{\sqrt{\alpha^2-\beta^2}}{2\alpha \delta K_1 (\delta\sqrt{\alpha^2-\beta^2})} exp\left(-\alpha\sqrt{\delta^2+(l-\mu)^2}+\beta(l-\mu)\right) \end{align*}

I use the HyperblicDistr package, my R for the first step, the estimation of the parameters (data):

hyperbFit(dat,hessian=TRUE)

The output is

Now I transform the parameters to my notation:

hyperbChangePars(from=1,to=2,c(-0.002494,0.002035,0.090747,0.204827))

Which gives

Now, I want to calculate the standard errors of my parameters, if I do:

summary(hyperbFit(dat,hessian=TRUE))

I get

But these are the standard errors in case of the parameter notation I don't want. So my first question:

1. To which parameters (in my notation) do the standard errors belong? So what is the standard error of $\hat{\alpha}, \hat{\beta}, \hat{\mu}$ and $\hat{\delta}$?
2. The values stay the same, right, so they do not change or?

Answer 1

1. Information about the relationship between $(\alpha,\beta,\mu,\delta)$ and $(\pi,\zeta,\delta,\mu)$ can be found in the following link:

https://github.com/sjp/GeneralizedHyperbolic/blob/master/R/hyperbFit.R

hyperbPi <- beta / sqrt(alpha^2 - beta^2)
zeta <- delta * sqrt(alpha^2 - beta^2)

Since the output of hyperbFit corresponds to the MLE, then they can be inverted directly using the aforementioned reparameterisation.

Although it is possible to calculate the corresponding errors for $(\alpha,\beta,\mu,\delta)$ using the delta method (if you do not know it, google will tell you), I think it requires a significant effort. I would rather go for using bootstrap as follows:

• Obtain B bootstrap samples of the estimators of $(\pi,\zeta,\delta,\mu)$.
• Transform those values into estimators of $(\alpha,\beta,\mu,\delta)$ using the reparameterisation mentioned above (you have to calculate it, recall that this is your thesis).
• Obtain bootstrap intervals using the sample of estimators of $(\alpha,\beta,\mu,\delta)$.

The following R code shows how to obtain the intervals for $(\pi,\zeta,\delta,\mu)$. Once you calculate the reparameterisation you can easily obtain those for the parameters of interest.

n = length(dat)
B=10000

sampB <- data.frame(matrix(ncol = 4, nrow = B))

for(i in 1:B){
print(i)
dat1 <-sample(dat,n,rep=T)
sampB[i,] <- hyperbFit(dat1,hessian=TRUE)$Theta
}

# quantile bootstrap interval for pi
c(quantile(sampB[,1],0.025),quantile(sampB[,1],0.975))

# quantile bootstrap interval for zeta
c(quantile(sampB[,2],0.025),quantile(sampB[,2],0.975))

# quantile bootstrap interval for delta
c(quantile(sampB[,3],0.025),quantile(sampB[,3],0.975))

# quantile bootstrap interval for mu
c(quantile(sampB[,4],0.025),quantile(sampB[,4],0.975))

As it was already mentioned in one of the answers to your questions, the MLE of this family of distributions converge slowly. These issues are inherited by the bootstrap algorithm. Any method that gives you stable results is worthy of suspicion.

Moreover, the use of standard errors in this family is completely unreliable since the distributions of the estimators are quite asymmetric as shown by the histogram of the bootstrap samples

hist(sampB[,1])
hist(sampB[,2])
hist(sampB[,3])
hist(sampB[,4])