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 $\beta$ are 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 $\alpha$ are 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 $\beta$ are:
$$
\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 $\alpha$ are 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 $\beta$ are:
$$
\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)
