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
