Skewness of Tweedie distribution Tweedie distributions are a family of distributions from the exponential dispersion family that have  power-law mean-variance relationship:
\begin{align}
\mathbb E[X] &= \mu \\
\operatorname{Var}[X]&=\phi \mu^p
\end{align}
What is the formula for skewness?

For small integer values of $p$, these are well-known distributions (Gaussian, Poisson, gamma, inverse Gaussian). The most interesting case is $1<p<2$ which corresponds to the compound Poisson-Gamma distribution that has point mass at zero and is continuous for $x>0$. The formulas for the density as well as for the parameters of Poisson and Gamma in terms of $\mu$ and $\phi$ can be found e.g. in  Dunn & @GordonSmyth (2005) Series evaluation of Tweedie exponential
dispersion model densities.
I couldn't find the skewness formula anywhere so I derived it myself, and am posting this Q&A to share the result.
 A: Exponential dispersion family is a broad family of distributions allowed in GLMs. The general form of the PDF can be written as follows:
$$f(x;\theta,\phi)=a(x,\phi)\exp\Big[\frac{1}{\phi}\big(x\theta-\kappa(\theta)\big)\Big].$$
The term $\kappa(\theta)$ is denoted with kappa because it is intimately related to the cumulants. Specifically, cumulant generating function (CGF) is given by
$$K(t;\theta,\lambda)=\frac{1}{\phi}\big(\kappa(\theta + t\phi)-\kappa(\theta)\big)$$
(see Wikipedia or Eq 2.6 in Jørgensen 1987, or Jørgensen's The Theory of Dispersion Models, 1997. Note that with $\phi=1$ the family reduces to the natural exponential family, see Wikipedia for its CGF.)
It follows that the first three cumulants are given by:
\begin{align}
\kappa_1 &= \kappa'(\theta)\\
\kappa_2 &= \phi\kappa''(\theta)\\
\kappa_3 &= \phi^2\kappa'''(\theta)
\end{align}
(Again note that for the natural exponential family cumulants are simply derivatives of $\kappa(\theta)$.)
For Tweedie distribution it must hold that
\begin{align}
\kappa_1 &= \kappa'(\theta) = \mu\\
\kappa_2 &= \phi\kappa''(\theta) = \phi\mu^p
\end{align}
so it follows that $$\kappa_3=\phi^2\kappa'''(\theta)=\phi^2(\kappa''(\theta))'=\phi^2(\mu^p)'=\phi^2p\mu^{p-1}\mu'=\phi^2p\mu^{p-1}\mu^p=\phi^2p\mu^{2p-1}.$$
Now we can compute skewness:
$$\operatorname{Skewness}[X]=\frac{\kappa_3}{\kappa_2^{3/2}}=\frac{\phi^2p\mu^{2p-1}}{(\phi\mu^p)^{3/2}}=\phi^{1/2}p\mu^{p/2-1}.$$
As a sanity check, this formula yields correct values for $p=0$, $p=1$, and $p=2$; these are skewness formulas for the Gaussian, Poisson, and gamma.
Let's verify that it works correctly for $1<p<2$:
# Tweedie random generation, using compound Poisson-Gamma representation
def tweediernd(n=1, p=1.5, phi=10, mu=1):
    # See Dunn & Smyth paper linked above for these formulas
    lambd = mu**(2-p)/(2-p)/phi   # Poisson rate
    alpha = -(2-p)/(1-p)          # gamma shape
    beta = phi*(p-1)*mu**(p-1)    # gamma scale

    x = np.zeros(n)
    for i in range(n):
        x[i] = np.sum(np.random.gamma(alpha, scale=beta, 
                      size=np.random.poisson(lambd)))
    return x

np.random.seed(42)
x = tweediernd(n=10000)
print('Mean:    ', np.mean(x))            # 1
print('Variance:', np.var(x))             # 10
print('Skewness:', scipy.stats.skew(x))   # sqrt(10)*1.5 = 4.74

This yields:
Mean:     0.996421833721
Variance: 9.86859188577
Skewness: 4.763172234662853

