# MGF of the absolute Value of a Skellam RV

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??

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)$$