MGF of the absolute Value of a Skellam RV

Question

I am trying to derive the moment generating function for the absolute value of a Skellam random variable $Skellam(\lambda_1, \lambda_2)$

Suppose $X_1 \sim Pois(\lambda_1)$ and $X_2 \sim Pois(\lambda_2)$ are independent. Deriving the MGF of $W = X_1 - X_2 \sim Skellam(\lambda_1, \lambda_2)$ is straight forward enough, since:

$M_{X_i}(t) = e^{\lambda_i (e^t - 1)}$

therefore

$M_{W}(t) = E(e^{tW})=E(e^{t(X_1-X_2)})=E(e^{tX_1-tX_2})=E(e^{tX_1})E(e^{-tX_2})=e^{\lambda_1 (e^t - 1)}e^{\lambda_2 (e^{-t} - 1)} = e^{-(\lambda_1 + \lambda_2) + \lambda_1 e^t + \lambda_2 e^{-t}}$

I'd like to find $M_{|W|}(t)$. My first thought was to use conditional expectation

$$ \begin{align} M_{|W|}(t) &= E(e^{t|W|}) \\ &= E(e^{t|X_1 - X_2|}) \\ &= E\left(e^{t(X_1 - X_2)}\right)P(X_1 > X_2) + E\left(e^{t(X_2 - X_1)}\right)P(X_1 \le X_2) \\ &= e^{-(\lambda_1 + \lambda_2) + \lambda_1 e^t + \lambda_2 e^{-t}} P(X_1 > X_2) + e^{-(\lambda_1 + \lambda_2) + \lambda_2 e^t + \lambda_1 e^{-t}} P(X_1 \le X_2) \end{align} $$

Which is where I get stuck. Can anyone help me get one or two steps beyond this?? Is there another, more elegant approach to this problem??

Answer 1

Let $Z$ follow a Skellam distribution with rates $\lambda_1$ and $\lambda_2$, and $K = |Z|$. The pmf for $K$ is: $$p(k; \lambda_1, \lambda_2) = \begin{cases} e^{-\lambda_1 - \lambda_2} \left( \left(\frac{\lambda_1}{\lambda_2}\right)^{\frac{k}{2}} I_k(2\sqrt{\lambda_1 \lambda_2}) + \left(\frac{\lambda_2}{\lambda_1}\right)^{\frac{k}{2}} I_{-k}(2\sqrt{\lambda_1 \lambda_2}) \right) &\text{if } k > 0\\ e^{-\lambda_1 - \lambda_2}I_0 (2\sqrt{\lambda_1 \lambda_2})& \text{if } k = 0\end{cases}$$

Here $I_k(a)$ is a modified Bessel function of the first kind and has the symmetry property $I_{k}(a) = I_{-k}(a)$, so the moment generating function of $k$ is

$$\begin{aligned} \mathcal{M}(t; \lambda_1, \lambda_2) = e^{-\lambda_1 - \lambda_2} \left(\sum_{k=0}^{\infty} e^{tk} I_k(2\sqrt{\lambda_1 \lambda_2}) \big[\big(\frac{\lambda_1}{\lambda_2}\big)^{\frac{k}{2}} + \big(\frac{\lambda_2}{\lambda_1}\big)^{\frac{k}{2}} \big] - I_0 (2\sqrt{\lambda_1 \lambda_2}) \right) \end{aligned} $$

Written in this form, recognize that the sum can be written in terms of a special function known as Marcum's $Q$ (used, for example, in the cdf of the noncentral $\chi^2$ distribution). A definition of $Q$ is:

$$ Q(\sqrt{2b},\sqrt{2a}) = e^{-a - b} \sum_{k=0}^\infty \left(\frac{a}{b}\right)^{\frac{k}{2}} I_k(2\sqrt{a b}) $$

So that the moment-generating function becomes:

$$\begin{aligned} \mathcal{M}(t;\lambda_1, \lambda_2) = e^{-\lambda_1 - \lambda_2} \big(&Q(\sqrt{2\lambda_2e^{-t}},\sqrt{2\lambda_1e^t}) e^{\lambda_1e^t + \lambda_2e^{-t}} + \\ &Q(\sqrt{2\lambda_1e^{-t}},\sqrt{2\lambda_2e^t}) e^{\lambda_2e^t + \lambda_1e^{-t}} - \\ &I_0 (2\sqrt{\lambda_1 \lambda_2})\big) \end{aligned}$$

---

Differentiating $Q(\sqrt{2\lambda_1e^{-t}}, \sqrt{2\lambda_2e^t})$ w.r.t. $t$ gives:

$$Q'(\sqrt{2\lambda_1e^{-t}}, \sqrt{2\lambda_2e^t}) = e^{ -\lambda_1 e^t - \lambda_2 e^{-t}} (\lambda_2e^{-t} I_0(2\sqrt{\lambda_1 \lambda_2 }) + \sqrt{\lambda_2 \lambda_1} I_1(2\sqrt{\lambda_1 \lambda_2 }) )$$

So, differentiating the MGF about $t=0$ yields the expected value of $k$:

$$ \begin{aligned} \mathbb{E}(K; \lambda_1, \lambda_2) = 2 &e^{-\lambda_1 - \lambda_2} \big( \lambda_2 I_0(2\sqrt{\lambda_1 \lambda_2 }) + \sqrt{\lambda_1 \lambda_2} I_1(2\sqrt{\lambda_1 \lambda_2 }) \big) + \\ &(\lambda_2 - \lambda_1)\left(1 - 2 Q(\sqrt{2\lambda_1}, \sqrt{2\lambda_2}) \right) \end{aligned} $$

If we take $\lambda_1 = \lambda_2 = \lambda$, the expectation simplifies to:

$$ \mathbb{E}_{\lambda}(K) = 2\lambda e^{-2\lambda} \left( I_0(2\lambda) + I_1(2\lambda) \right) $$