Expected value of the SHASH distribution? The sinh-arcsinh (SHASH) distribution has a pdf as follows:
$f(x) = {\delta cosh(\omega)\over \sqrt{1+({x-\theta \over \sigma})^2}}\phi[sinh(\omega)]$
where $\omega=\gamma+\delta sinh^{-1}({x-\theta \over \sigma})$
and $\phi()$ is the standard normal pdf.
Note:  $-\infty<\gamma,x,\theta<\infty; 0<\delta, \sigma$
Does anyone know what the expected value and variance of it is?  I can see it is closely related to the Johnson SU distribution which has the mean, variance, median, and quantiles analytically defined.  It also has the Normal distribution as a special case when $\gamma = 0$ and $ \delta = 1$.  I know some distributions don't have these available in closed form analytic expressions....but I don't know if that's the case here.  If available, it would be very helpful if someone knows what they are!
Update:  It appears the moments are published, but with the caveat that they depend on the modified Bessel function of the second kind ($P_\nu)$ below.  E.g.
$E(X_{\gamma, \delta})=-sinh(\gamma/ \delta)P_{1/\delta}$
$Var(X_{\gamma, \delta})={1 \over 2}(cosh({2\gamma \over \delta})P_{2/\delta} - 1) - \mu_{\gamma, \delta}^2$
However -- this doesn't take into account the center, scale in the full version ($\theta, \sigma$).  Can someone help modify the mean and var with center and scale included?
 A: Indeed, the moments of this distribution are already calculated in the original paper
https://www.jstor.org/stable/27798865
For the case with location and scale parameters, you just need to use the usual properties of the location-scale family. Let $Z = \mu + \sigma X$, where $X \sim SHASH(\gamma,\delta)$, then
$$E[Z] = \mu + \sigma E[X].$$
Also,
$$Var[Z] = \sigma^2 Var[X].$$
 [Here]  is an implementation of this distribution in R. Using their command  rsas, you can simulate from a distribution with specific parameter values and approximate its mean and variance using mean() and var().
sim <- rsas(1e5,3,2,-1,1.5)
mean(sim)
var(sim)

Alternatively, you can use their function dsas to approximate the mean using numerical integration:
mean_sas <- function(mu,sigma,epsilon,delta){
  tempf <- Vectorize(function(x) x*dsas(x,mu,sigma,epsilon,delta))
  val <- integrate(tempf,-Inf,Inf)$value
  return(val)
}

mean_sas(3,2,-1,1.5) 

Similarly for the variance:
var_sas<- function(mu,sigma,epsilon,delta){
  tempf1 <- Vectorize(function(x) x*dsas(x,mu,sigma,epsilon,delta))
  tempf2 <- Vectorize(function(x) x^2*dsas(x,mu,sigma,epsilon,delta))
  val <- integrate(tempf2,-Inf,Inf)$value - (integrate(tempf1,-Inf,Inf)$value)^2
  return(val)
}


var_sas(3,2,-1,1.5)

Finally, the same ideas can be applied to [their alternative version] of the SHASH distribution.
