Problem with the density of the compound distribution My problem is to calculate $E[\max(S-5000, 0)]$ where 
$$S = \sum_{i=1}^{N} X_i,$$ 
$N$ is a random variable with geometric distribution, parametrized as follows:
$$P(N=n) = \frac{\beta^n}{(1+\beta)^{n+1}}, ~~~~~ n=\{0,1,2, ...\}$$
and ${X_1, X_2...}$ are all independent and exponentially distibuted with density:
$$f_X(x) = \lambda e^{-\lambda x}.$$
I came up with an idea that I can find the density function of $S$ which is (according to my knowledge):
$$f_S(s) = \sum_{n=0}^{\infty} f_{X_1 + X_2 + ... +X_n}(s) \cdot P(N=n).$$
However, I encountered a problem - $N$ may have values $\{0,1,2,...\}$ (natural numbers WITH zero) but I can't calculate the density of the sum of zero random variables (the number of random variables to be sumed is zero when $n=0$). 
Is it possible to calculate such density? It is there any other way to calculate $E[\max(S-5000, 0)]$?
 A: First note that the sum of $n$ i.i.d. exponential random variables has an Erlang distribution, that is, for any $n\geqslant1$, the sum $\sum_{i=1}^n X_i$ has probability density 
$$\frac{(\lambda t)^{n-1}}{(n-1)!}\lambda e^{-\lambda t}\mathsf 1_{(0,\infty)}(t).  $$
(How sum of exponential variables is a gamma variable) It follows then that the distribution of $S$ is given by $\nu_S = \nu_S^d + \nu_S^c$, where for each Borel set $B$, $$\nu_S^d(B)  =\mu(B\cap\{0\})\mathbb P(N=0)   $$ where $\mu$ is counting measure, and $$\nu_S^c(B)=\int_B f_S\ \mathsf d\lambda$$
where $\lambda$ is Lebesgue measure. We compute $$\nu_S^d(B) =\mu(B\cap\{0\})(1+\beta)^{-1}  $$ and
\begin{align}
f_S(t) &= \sum_{n=0}^\infty f_{S\mid N=n}(t)\mathbb P(N=n)\\
&= \sum_{n=1}^\infty \frac{(\lambda t)^{n-1}}{(n-1)!}\lambda e^{-\lambda t}\frac{\beta^n}{(1+\beta)^{n+1}}\\
&= \frac{\beta \lambda e^{-\lambda t}}{(1+\beta)^2}\sum_{n=0}^\infty\frac{\left(\frac{\beta\lambda t}{1+\beta}\right)^n}{n!}\\
&=\frac{\beta\lambda}{(1+\beta)^2} \exp\left(-\frac{\lambda t}{1+\beta}\right)\mathsf 1_{(0,\infty)}(t).
\end{align}
It follows that 
\begin{align}
\mathbb E[(S-5000)^+] &= \int_{5000}^\infty (t-5000)\frac{\beta\lambda}{(1+\beta)^2} e^{-\frac{\lambda t}{1+\beta}}\ \mathsf dt\\
&= \frac\beta\lambda \exp\left(-\frac{5000\beta\lambda}{1+\beta}\right).
\end{align}
