# Standard errors of hyperbolic distribution estimates using delta-method?

I want to calculate the standard errors of a fitted hyperbolic distribution.

In my notation the density is given by \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 am using the HyperbolicDistr package in R. I estimate the parameters via the following command:

hyperbFit(mydata,hessian=TRUE)


This gives me a wrong parameterization. I change it into my desired parameterization with the hyperbChangePars(from=1,to=2,c(mu,delta,pi,zeta)) command. Then I want to have the standard errors of my estimates, I can get it for the wrong parameterization with the summary command. But this gives me the standard errors for the other parameterization. According to this thread I have to use the delta-method (I do not want to use bootstrap or cross-validation or so).

The hyperbFit code is here. And the hyperbChangePars is here. Therefore I know, that $\mu$ and $\delta$ stay the same. Therefore also the standard errors are the same, right?

For transforming $\pi$ and $\zeta$ into $\alpha$ and $\beta$ I need the relationship between them. According to the code this is done as follows:

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


So how do I have to code the delta-method to get the desired standard errors?

EDIT: I am using these data. I first perform the delta-method according to this thread.

# fit the distribution

hyperbfitdb<-hyperbFit(mydata,hessian=TRUE)
hyperbChangePars(from=1,to=2,hyperbfitdb$Theta) summary(hyperbfitdb)  summary(hyperbfitdb) gives the following output: Data: mydata Parameter estimates: pi zeta delta mu 0.0007014 1.3779503 0.0186331 -0.0001352 ( 0.0938886) ( 0.9795029) ( 0.0101284) ( 0.0035774) Likelihood: 615.992 Method: Nelder-Mead Convergence code: 0 Iterations: 315  and hyperbChangePars(from=1,to=2,hyperbfitdb$Theta) gives the following output:

   alpha.zeta     beta.zeta   delta.delta         mu.mu
73.9516898823  0.0518715378  0.0186331187 -0.0001352342


now I define the variables in the following way:

pi<-0.0007014
lzeta<-log(1.3779503)
ldelta<-log(0.0186331)


I now run the code (second edit) and get the following result:

> se.alpha
[,1]
[1,] 13.18457
> se.beta
[,1]
[1,] 6.94268


Is this correct? I am wondering about the following: If I use a bootstrap-algorithm in the following way:

B = 1000 # number of bootstraps

alpha<-NA
beta<-NA
delta<-NA
mu<-NA

# Bootstrap
for(i in 1:B){
print(i)
subsample = sample(mydata,rep=T)
hyperboot <- hyperbFit(subsample,hessian=FALSE)
hyperboottransfparam<- hyperbChangePars(from=1,to=2,hyperboot$Theta) alpha[i] = hyperboottransfparam[1] beta[i] = hyperboottransfparam[2] delta[i] = hyperboottransfparam[3] mu[i] = hyperboottransfparam[4] } # hist(beta,breaks=100,xlim=c(-200,200)) sd(alpha) sd(beta) sd(delta) sd(mu)  I get 119.6 for sd(alpha) and 35.85 for sd(beta). The results are very different? Is there a mistake or what is the problem here? ## 2 Answers In the following solution, I assume hyperbPi to be$\pi$. Also, the variances used in the approximations below are simply the squared standard errors calculated by summary after hyperbFit, so$\mathrm{Var}(X)=\mathrm{SE}(X)^2$. In order to calculate the approximation using the delta-method, we need the partial derivatives of the transformation function s$g_{\alpha}(\zeta, \pi, \delta)$and$g_{\beta}(\zeta, \pi, \delta)$. The transformation functions for$\alpha$and$\betaare given by: \begin{align} g_{\alpha}(\zeta, \pi, \delta) &=\frac{\zeta\sqrt{1 + \pi^{2}}}{\delta}\\ g_{\beta}(\zeta, \pi, \delta) &= \frac{\zeta\pi}{\delta}\\ \end{align} The partial derivatives of the transformation function for\alphaare then: \begin{align} \frac{\partial}{\partial \zeta} g_{\alpha}(\zeta, \pi, \delta) &=\frac{\sqrt{1+\pi^{2}}}{\delta}\\ \frac{\partial}{\partial \pi} g_{\alpha}(\zeta, \pi, \delta) &= \frac{\pi\zeta}{\sqrt{1+\pi^{2}}\delta }\\ \frac{\partial}{\partial \delta} g_{\alpha}(\zeta, \pi, \delta) &= -\frac{\sqrt{1+\pi^{2}}\zeta}{\delta^{2}}\\ \end{align} The partial derivatives of the transformation function for\betaare: \begin{align} \frac{\partial}{\partial \zeta} g_{\beta}(\zeta, \pi, \delta) &=\frac{\pi}{\delta}\\ \frac{\partial}{\partial \pi} g_{\beta}(\zeta, \pi, \delta) &= \frac{\zeta}{\delta }\\ \frac{\partial}{\partial \delta} g_{\beta}(\zeta, \pi, \delta) &= -\frac{\pi\zeta}{\delta^{2}}\\ \end{align} Applying the delta-method to the transformations, we get the following approximation for the variance of\alpha$(take square roots to get the standard errors): $$\mathrm{Var}(\alpha)\approx \frac{1+\pi^{2}}{\delta^{2}}\cdot \mathrm{Var}(\zeta)+\frac{\pi^{2}\zeta^{2}}{(1+\pi^{2})\delta^{2}}\cdot \mathrm{Var}(\pi) + \frac{(1+\pi^{2})\zeta^{2}}{\delta^{4}}\cdot \mathrm{Var}(\delta) + \\ 2\times \left[ \frac{\pi\zeta}{\delta^{2}}\cdot \mathrm{Cov}(\pi,\zeta) - \frac{(1+\pi^{2})\zeta}{\delta^{3}}\cdot \mathrm{Cov}(\delta,\zeta)- \frac{\pi\zeta^{2}}{\delta^{3}}\cdot \mathrm{Cov}(\delta,\pi)\right]$$ The approximated variance of$\beta$is: $$\mathrm{Var}(\beta)\approx \frac{\pi^{2}}{\delta^{2}}\cdot \mathrm{Var}(\zeta) + \frac{\zeta^{2}}{\delta^{2}}\cdot \mathrm{Var}(\pi) + \frac{\pi^{2}\zeta^{2}}{\delta^{4}}\cdot \mathrm{Var}(\delta) + \\ 2\times \left[ \frac{\pi\zeta}{\delta^{2}}\cdot \mathrm{Cov}(\pi,\zeta) - \frac{\pi^{2}\zeta}{\delta^{3}}\cdot \mathrm{Cov}(\delta, \zeta) - \frac{\pi\zeta^{2}}{\delta^{3}}\cdot \mathrm{Cov}(\pi, \delta) \right]$$ Coding in R The fastest way to calculate the above approximations is using matrices. Denote$D$the row vector containing the partial derivatives of the transformation function for$\alpha$or$\beta$with respect to$\zeta, \pi, \delta$. Further, denote$\Sigma$the$3\times 3$variance-covariance matrix of$\zeta, \pi, \delta$. The covariance matrix can be retrieved by typing vcov(my.hyperbFit) where my.hyperbFit is the fitted function. The above approximation of the variance of$\alpha$is then $$\mathrm{Var}(\alpha)\approx D_{\alpha}\Sigma D_{\alpha}^\top$$ The same is true for the approximation of the variance of$\beta$. In R, this can be easily coded like this: #----------------------------------------------------------------------------- # The row vector D of the partial derivatives for alpha #----------------------------------------------------------------------------- D.alpha <- matrix( c( sqrt(1+pi^2)/delta, # differentiate wrt zeta ((pi*zeta)/(sqrt(1+pi^2)*delta)), # differentiate wrt pi -(sqrt(1+pi^2)*zeta)/(delta^2) # differentiate wrt delta ), ncol=3) #----------------------------------------------------------------------------- # The row vector D of the partial derivatives for beta #----------------------------------------------------------------------------- D.beta <- matrix( c( (pi/delta), # differentiate wrt zeta (zeta/delta), # differentiate wrt pi -((pi*zeta)/delta^2) # differentiate wrt delta ), ncol=3) #----------------------------------------------------------------------------- # Calculate the approximations of the variances for alpha and beta # "sigma" denotes the 3x3 covariance matrix #----------------------------------------------------------------------------- var.alpha <- D.alpha %*% sigma %*% t(D.alpha) var.beta <- D.beta %*% sigma %*% t(D.beta) #----------------------------------------------------------------------------- # The standard errors are the square roots of the variances #----------------------------------------------------------------------------- se.alpha <- sqrt(var.alpha) se.beta <- sqrt(var.beta)  ## Using$\log(\zeta)$and$\log(\delta)$If the standard errors/variances are only available for$\zeta^{*}=\log(\zeta)$and$\delta^{*}=\log(\delta)$instead of$\zeta$and$\delta, the transformation functions change to: \begin{align} g_{\alpha}(\zeta^{*}, \pi, \delta^{*}) &=\frac{\exp(\zeta^{*})\sqrt{1 + \pi^{2}}}{\exp(\zeta^{*})}\\ g_{\beta}(\zeta^{*}, \pi, \delta^{*}) &= \frac{\exp(\zeta^{*})\pi}{\exp(\delta^{*})}\\ \end{align} The partial derivatives of the transformation function for\alphaare then: \begin{align} \frac{\partial}{\partial \zeta^{*}} g_{\alpha}(\zeta^{*}, \pi, \delta^{*}) &=\sqrt{1+\pi^{2}}\exp(-\delta^{*}+\zeta^{*})\\ \frac{\partial}{\partial \pi} g_{\alpha}(\zeta^{*}, \pi, \delta^{*}) &=\frac{\pi\exp(-\delta^{*}+\zeta^{*})}{\sqrt{1+\pi^{2}}} \\ \frac{\partial}{\partial \delta^{*}} g_{\alpha}(\zeta^{*}, \pi, \delta^{*}) &=-\sqrt{1+\pi^{2}}\exp(-\delta^{*}+\zeta^{*})\\ \end{align} The partial derivatives of the transformation function for\betaare: \begin{align} \frac{\partial}{\partial \zeta^{*}} g_{\beta}(\zeta^{*}, \pi, \delta^{*}) &=\pi\exp(-\delta^{*}+\zeta^{*})\\ \frac{\partial}{\partial \pi} g_{\beta}(\zeta^{*}, \pi, \delta^{*}) &=\exp(-\delta^{*}+\zeta^{*})\\ \frac{\partial}{\partial \delta^{*}} g_{\beta}(\zeta^{*}, \pi, \delta^{*}) &=-\pi\exp(-\delta^{*}+\zeta^{*})\\ \end{align} Applying the delta-method to the transformations, we get the following approximation for the variance of\alpha$: $$\mathrm{Var}(\alpha)\approx (1+\pi^{2})\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Var}(\zeta^{*})+\frac{\pi^{2}\exp(-2\delta^{*}+2\zeta^{*})}{1+\pi^{2}}\cdot \mathrm{Var}(\pi) + (1+\pi^{2})\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Var}(\delta^{*}) + \\ 2\times \left[ \pi\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Cov}(\pi,\zeta^{*}) - (1+\pi^{2})\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Cov}(\delta^{*},\zeta^{*}) - \pi\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Cov}(\delta^{*},\pi)\right]$$ The approximated variance of$\beta$is: $$\mathrm{Var}(\beta)\approx \pi^{2}\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Var}(\zeta^{*})+\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Var}(\pi) + \pi^{2}\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Var}(\delta^{*}) + \\ 2\times \left[\pi\exp(-2\delta^{*}+2\zeta^{*}) \cdot \mathrm{Cov}(\pi,\zeta^{*}) -\pi^{2}\exp(-2\delta^{*}+2\zeta^{*})\cdot \mathrm{Cov}(\delta^{*},\zeta^{*}) -\pi\exp(-2\delta^{*}+2\zeta^{*}) \cdot \mathrm{Cov}(\delta^{*},\pi)\right]$$ Coding in R 2 This time, sigma denotes the covariance matrix but including the variances and covariances for$\zeta^{*}=\log(\zeta)$and$\delta^{*}=\log(\delta)$instead of$\zeta$and$\delta$. #----------------------------------------------------------------------------- # The row vector D of the partial derivatives for alpha #----------------------------------------------------------------------------- D.alpha <- matrix( c( sqrt(1+pi^2)*exp(-ldelta + lzeta), # differentiate wrt lzeta ((pi*exp(-ldelta + lzeta))/(sqrt(1+pi^2))), # differentiate wrt pi (-sqrt(1+pi^2)*exp(-ldelta + lzeta)) # differentiate wrt ldelta ), ncol=3) #----------------------------------------------------------------------------- # The row vector D of the partial derivatives for beta #----------------------------------------------------------------------------- D.beta <- matrix( c( (pi*exp(-ldelta + lzeta)), # differentiate wrt lzeta exp(-ldelta + lzeta), # differentiate wrt pi (-pi*exp(-ldelta + lzeta)) # differentiate wrt ldelta ), ncol=3) #----------------------------------------------------------------------------- # Calculate the approximations of the variances for alpha and beta # "sigma" denotes the 3x3 covariance matrix with log(delta) and log(zeta) #----------------------------------------------------------------------------- var.alpha <- D.alpha %*% sigma %*% t(D.alpha) var.beta <- D.beta %*% sigma %*% t(D.beta) #----------------------------------------------------------------------------- # The standard errors are the square roots of the variances #----------------------------------------------------------------------------- se.alpha <- sqrt(var.alpha) se.beta <- sqrt(var.beta)  • @BenBohold The terms inside the brackets before the covariance terms are the products of the respective partial derivatives of the transformation functions. For example: The term before$\mathrm{Cov}(\pi, \zeta)$is the product of the partial derivative wrt$\pi$multiplied by the partial derivative wrt$\zeta$. In the case of$\beta$this would be:$\zeta/\delta \times \pi/\delta = (\zeta\cdot\pi)/\delta^{2}$. Aug 16 '13 at 7:42 • @BenBohold Strange, but no problem. Try to calculate the covariance matrix in this way: varcov <- solve(hyperbfitalv$hessian). Does this work? After that, you will have to select the sub-matrix containing only $\pi, \zeta, \delta$. The easiest way I could help you would be if you provided a fully working example with data (you don't have to provide all your data). Aug 16 '13 at 19:33
• A big thanks for your answer, but this is EXACTLY the problem, because the parameterization of this hessian is for log(delta) and log(zeta) and not for delta and zeta! See my follow-up posts: stats.stackexchange.com/questions/67595/… and especially the answer of CT Zhu here stats.stackexchange.com/questions/67602/… Aug 16 '13 at 19:46
• one has to get the hessian of pi, log(zeta), log(delta), and mu into the hessian of pi, zeta, delta, and mu. Do you know how this could be done? Aug 16 '13 at 19:48
• I also tried to do it with the "wrong" hessian, so with log(delta) and log(zeta), after this I selected the sub-matrix and did the calculations. The results were not correct, because the values were way too large, so like 60 000 or so. Aug 16 '13 at 19:53

Possible duplicate: Standard errors of hyperbFit?

I could bet some accounts belong to the same person ...

• I already linked that thread in my post! And if you really read my question you will see, that I DO NOT WANT TO use the bootstrap algorithm, but I am asking about the delta-method. Aug 12 '13 at 16:36