Compound Binomial-Exponential: Closed form for the PDF? Setup:
Consider the random variable $Y_N$ derived from 
$$Y_N = \sum_{i=0}^N X_i$$
where $N$ is a random variable with distribution $p_n = {M \choose n} p^n q^{M-n}$ (binomial i.e. $M$ Bernoulli trials) and the $X_i$ are random variables with distribution $f(x) = \lambda e^{-\lambda x}$ (exponential). 
I'd like to get the pdf of $y$. 
Previously I asked this question which led me to expect the distribution is the sum 
$$p(y) =  \sum_{i=1}^M {M \choose i} p^i q^{M-i} \frac{\lambda^iy^{i-1}e^{-\lambda y}}{(i-1)!} $$
or the integral (inverse Fourier transform of the characteristic function)
$$p(y) = \int_{-\infty}^{\infty} \frac{dk}{2\pi} e^{-iky}(\frac{p\lambda}{\lambda-ik}+q)^M = \int_{-\infty}^{\infty} \frac{dk}{2\pi}e^{-iky}(\frac{i\lambda + kq}{i\lambda + k})^M$$
The second equality follows from $p+q=1$. 
I'd like to solve either one of these problems in closed form if possible. The moment generating function these derive from is
$$ p(s) = (\frac{p\lambda}{\lambda-s}+q)^M.$$
The characteristic function is $p(s=ik)$. 
Question:
What is the distribution $p(y)$ of $Y_N$? 
Attempts:
I tried binomial expanding $p(s)$ in the fourier transform form, which leads to a sum over integrals like $\int_{-\infty}^{\infty} du u^{-m}e^{-a u}$ where $m$ is a positive integer -- almost a sum of gamma functions but with negative parameters and over the wrong limits. 
I tried manipulaing the sum form into something I could use the binomial formula on, but I can't obtain binomial coefficients from the factorial prefactors in the sum. 
Any thoughts? I had thought a compound binomial-exponential process would have been well studied with a closed form solution but I couldn't find anything immediately relevant. 
Other:
This paper explores a very similar problem but I'll have to read carefully to see how it relates (link).
 A: I was able to compute the integral you presented using contour integration in the complex plane.
For the case of $\lambda > 0$, $y <0$, it is easy to show via contour integration that
$$\int_{-\infty}^\infty \dfrac{dk}{2\pi}e^{-iyk}\left(\dfrac{i\lambda p}{k - (-i\lambda)}+q\right)^M = 0$$
For the more interesting case of $\lambda > 0$, $y > 0$, one can show via contour integration, and algebraic manipulations to generate a Laurent Series, that
$$\begin{align*}\int_{-\infty}^\infty \dfrac{dk}{2\pi}e^{-iyk}\left(\dfrac{i\lambda p}{k - (-i\lambda)}+q\right)^M &= -2\pi i \space \mathrm{Res}_{z=-i\lambda}\left[\dfrac{1}{2\pi}e^{-iyz}\left(\dfrac{i\lambda p}{z - (-i\lambda)}+q\right)^M\right]\\
\\
&= -2\pi i \space \mathrm{Res}_{z=-i\lambda}\left[\dfrac{1}{2\pi}e^{-iyz}\sum_{j=0}^{M}\binom{M}{j}\left(\dfrac{i\lambda p}{z - (-i\lambda)}\right)^jq^{M-j}\right]\\
\\
&= -2\pi i \space \mathrm{Res}_{z=-i\lambda}\left[\dfrac{1}{2\pi}\left(\sum_{n=0}^\infty \dfrac{\left(-iy\right)^n e^{-\lambda y}}{n!}\left(z-\left(-i\lambda\right)\right)^n\right)\left(\sum_{j=0}^{M}\binom{M}{j}\left(\dfrac{i\lambda p}{z - (-i\lambda)}\right)^jq^{M-j}\right)\right]\\
\\
&= -2\pi i \dfrac{1}{2\pi}\sum_{j=1}^M \dfrac{\left(-iy\right)^{j-1}e^{-\lambda y}}{\left(j-1\right)!}\binom{M}{j}\left(i\lambda p\right)^jq^{M-j} \\
\\
&= \sum_{j=1}^M \binom{M}{j} p^j q^{M-j}\dfrac{\lambda^j y^{j-1} e^{-\lambda y}}{\left(j-1\right)!}\\
\\
&\mbox{(confirmation of exactly what you already thought to be true)}\\
\\
&= \dfrac{q^M e^{-\lambda y}}{y}\sum_{j=1}^M \binom{M}{j} \left(\dfrac{\lambda y p}{q}\right)^j \dfrac{1}{\left(j-1\right)!}\\
\end{align*}$$
For the case of $\lambda > 0$, $y =0$, the integral diverges due to the "$+q$".  The integral diverges at this single value of $y$ with a weight of $q^M$, given some hindsight with the convolution theorem and the definition of a CDF.
The integral really does not provide any new information regarding the series.  You need to manipulate the series form into something you find palatable.
The above development of the Laurent Series, to find the residue, does provide some insights into how portions of each term in the final sum arise.  The binomial coefficients come from the binomial expansion, as you would expect.  The factorial factors come from the Taylor Series expansion of the exponential function.  The final summation comes from multiplying series terms together and adding to find the coefficient of the $(z-(-i\lambda))^{-1}$ term of the Laurent Series.
Consolidating the results for all real values of $y$, the expression should really be written as:
$$\begin{align*}p(y) &= q^M\delta(y)+H(y)\dfrac{q^M e^{-\lambda y}}{y}\sum_{j=1}^M \binom{M}{j} \left(\dfrac{\lambda y p}{q}\right)^j \dfrac{1}{\left(j-1\right)!}\\
\\
&= q^M\delta(y)+ H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} \left(\dfrac{p}{q}\lambda y\right)^{j} \dfrac{1}{j!}\\
\end{align*}$$
Where $H(y)$ is the Heaviside unit step and $\delta(y)$ is the Dirac Delta function.
Using this expression, convolving the $M=1$ case with itself $M-1$ times, yields the correct expression for the $M^{th}$ case, as expected given the convolution theorem of the Fourier Transform.
Using the above expression for $p(y)$, we find
$$\int_{-\infty}^\infty p(y) dy = (q+p)^M =1$$
which is what one would expect. 
To replace the series, Wolfram Alpha comes up with this "closed" form in terms of a hypergeometric function:
$$\begin{align*}p(y) &= q^M\delta(y)+H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} \left(\dfrac{p}{q}\lambda y \right)^{j} \dfrac{1}{j!}\\
\\
&= q^M\delta(y)+ H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y} M \space {}_1F_1 \left(1-M;2;-\dfrac{p}{q}\lambda y\right)
\end{align*}$$
Update to derive the CDF from the PDF
To derive the CDF, let's rewrite the PDF a little bit.  I'll be using the following expression involving the Incomplete Upper Gamma Function in the process
$$\dfrac{\partial}{\partial y}\Gamma(j,\lambda y) =\dfrac{\partial}{\partial (\lambda y)}\Gamma(j,\lambda y) \cdot \dfrac{\partial}{\partial y}\lambda y = -(\lambda y)^{j-1}e^{-\lambda y}\lambda$$
So working the PDF a bit
$$\begin{align*}p(y) &= q^M\delta(y)+H(y) \dfrac{p}{q} q^{M} \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} \left(\dfrac{p}{q}\lambda y \right)^{j} \dfrac{1}{j!}\\
\\
&= \delta(y)q^M+H(y) \lambda e^{-\lambda y}\sum_{j=0}^{M-1} \binom{M}{j+1} q^{M-(j+1)}p^{j+1}\left(\lambda y \right)^{j} \dfrac{1}{j!}\\
\\
&= \delta(y)q^M+H(y) \lambda e^{-\lambda y}\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{1}{(j-1)!}\left(\lambda y \right)^{j-1} \\
\\
&= \delta(y)q^M-H(y) \sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{1}{(j-1)!}(-1)\left(\lambda y \right)^{j-1} e^{-\lambda y}\lambda \\
\\
&= \delta(y)q^M-H(y) \sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\partial}{\partial y}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)} \\
\end{align*}$$
To integrate the PDF to get the CDF, I'll be using the following relations:
$$\int_{-\infty}^{x} H(t)f'(t) dt = H(x)\int_0^x f'(t) dt = H(x)\left[f(x) - f(0)\right]$$
$$\int_{-\infty}^x \delta(t) f(t)dt = H(x)f(0)$$
So, integrating the PDF to find the CDF:
$$\begin{align*}P(y) &= \int_{-\infty}^y p(t) dt\\
\\
&= \int_{-\infty}^y \delta(t)q^M dt - \sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\int_{-\infty}^y H(t)\dfrac{\partial}{\partial t}\dfrac{\Gamma(j,\lambda t)}{\Gamma(j)} dt\\
\\
&= H(y)q^M -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j} H(y)\left[\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}-\dfrac{\Gamma(j,0)}{\Gamma(j)}\right]\\
\\
&= H(y)\left[q^M -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j} \left(\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}-1\right)\right]\\
\\
&= H(y)\left[q^M +\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j} -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}\right]\\
\\
&= H(y)\left[(q+p)^M -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}\right]\\
\\
P(y) &= H(y)\left[1 -\sum_{j=1}^{M} \binom{M}{j} q^{M-j}p^{j}\dfrac{\Gamma(j,\lambda y)}{\Gamma(j)}\right]\\
\end{align*}$$
