Variance of the random sum of a Poisson? We have that $N$ and $X_1, X_2, \dots$ are all independent. We also have $\operatorname{E} [X_j] = \mu$ and $\operatorname{Var}[X_j] = σ^2$.
Then, we introduce an integer–valued random variable, $N$, which is the random sum such that:
$$Z = \sum_{j=1}^{N+1}X_j.$$
Assuming that $N$ is distributed $\sim\mathrm{Poisson}(\lambda)$, what is the first moment and what is the variance of $Z$?
For a normal Poisson distribution, I know the variance is just $\lambda$, as is the mean. I'm having trouble understanding the implication of having the bounds be poisson distributed. Normally, I would just say "variance of the sum is the sum of the variance," but I don't think that's how it works with random sums. Any hints/guidance appreciated.
 A: \begin{align}
\operatorname{var}(Z) & = \operatorname{var}(\operatorname{E}(Z\mid N)) + \operatorname{E}(\operatorname{var}(Z\mid N)) & & (\text{This is the law of total variance.}) \\[10pt]
& = \operatorname{var}(N\mu) + \operatorname{E}(N\sigma^2) \\[10pt]
& = \mu^2 \operatorname{var}(N) + \sigma^2 \operatorname{E}(N) \\[10pt]
& = \mu^2 \lambda + \sigma^2 \lambda = \lambda \operatorname{E}(X_1^2).
\end{align}
Indeed, generally the $n$th cumulant of a compound Poisson distribution is the mean of the simple Poisson distribution times the $n$th raw moment of the distribution that gets compounded.
PS: The above applies if the sum is from $0$ to $N$; I'll leave it as an exercise to figure out whether something needs to change if it's from $1$ to $N+1$.
A: Hint: rewrite
$$
Z = \sum_{n=1}^\infty X_n\mathbb{1}_{\{N+1 \geq n\}}
$$ and apply your "expectation of the sum is the sum of the expectations" idea.

Following the comment below, more detail.
Write $Y_n = \mathbb{1}_{\{N+1 \geq n\}}$, which is independent of $X_n$. 


*

*Then for the expectation, you have $$\mathbb{E}[Z] = \sum_{n=1}^\infty \mathbb{E}[X_nY_n] = \sum_{n=1}^\infty \mathbb{E}[X_n] \mathbb{E}[Y_n]= \sum_{n=1}^\infty \mu \mathbb{P}\{N \geq n-1\} = \mu\sum_{n=0}^\infty \mathbb{P}\{N \geq n\} $$ using independence of $X_n$ and $Y_n$, which yields $$\mathbb{E}[Z] = \mu\lambda$$ since $\sum_{n=0}^\infty\mathbb{P}\{N \geq n\} = \mathbb{E}[N].$ 


*

*The variance will be a bit (a lot?) less straightforward, though, since the $Y_n$'s are not independent. But you have $\mathbb{E}[Z]^2$ already, so it only remains to compute $\mathbb{E}[Z^2]$ — it's not very enjoyable, but you can do so by expanding the sum:
$$\begin{align}
\mathbb{E}[Z^2] &= \mathbb{E}\left[\sum_{n=1}^\infty\sum_{m=1}^\infty X_nY_nX_mY_m \right ] =
\sum_{n=1}^\infty\sum_{m=1}^\infty \mathbb{E}\left[X_nY_nX_mY_m \right ]
\\
&=
\sum_{n=1}^\infty \mathbb{E}\left[X^2_n\right]\mathbb{E}\left[Y_n^2\right ]
+ 2\sum_{n=1}^\infty\sum_{m=n+1}^\infty \mathbb{E}\left[X_n\right]\mathbb{E}\left[X_m\right]\mathbb{E}\left[Y_n Y_m \right ]
\\ &=
\sum_{n=1}^\infty \mathbb{E}\left[X^2_n\right]\mathbb{E}\left[Y_n\right ]
+ 2\mu^2\sum_{n=1}^\infty\sum_{m=n+1}^\infty \mathbb{E}\left[Y_nY_m \right ]
\end{align}$$
using that $Y_n^2=Y_n$ (it's a random variable being either $0$ or $1$). Now, observe that for $m \geq n$, $Y_nY_m = Y_m$, so you get
$$\begin{align}
\mathbb{E}[Z^2] 
&=
(\sigma^2+\mu^2)\sum_{n=1}^\infty \mathbb{E}\left[Y_n\right ]
+ 2\mu^2\sum_{n=1}^\infty\sum_{m=n+1}^\infty \mathbb{E}\left[Y_m \right ] \\
&=
(\sigma^2+\mu^2)\sum_{n=1}^\infty \mathbb{E}\left[Y_n\right ]
+ 2\mu^2\sum_{m=1}^\infty (m-1)\mathbb{E}\left[Y_m \right ]\\
&=
(\sigma^2+\mu^2)\sum_{n=0}^\infty \mathbb{P}\{N \geq n\}
+ 2\mu^2\sum_{m=0}^\infty m\mathbb{P}\{N \geq m\} \\
&=
(\sigma^2+\mu^2)\lambda
+ 2\mu^2\sum_{m=0}^\infty m\mathbb{P}\{N \geq m\} \\
\end{align}$$
and you can continue by manipulating the last sum (it's quite tedious, but it works).


A: If $T = X_1 + X_2 + \cdots +X_N$, where $X_i$ are iid and $N$ is
independent of $X_i$, then a standard result, usually derived
via a conditioning argument, is that $E(T) = E(N)E(X)$ and 
$V(T) = E(N)V(X) + V(N)[E(X)]^2.$ (Search the Internet for
proofs under 'Random sum of random variables'; the UNL page
seems complete.) A small adjustment will
deal with your question; see comments in the simulation code below. (Question edited: What used to be $S$ is
now $Z,$ original $S$ used below.)
In addition to its mean and variance, the $distribution$ of $S$
can be approximated through simulation. Suppose that $X_i \sim Exp(rate=1/3)$ and $N \sim Pois(\lambda = 5).$ Then a simulation
of 100,000 sums $S$ can be performed as follows. For example,
one can approximate $P(S > 20),$ which is not available from
the formulas above.
 m = 10^5;  s = numeric(m);  rate=1/3;  lam=5
 for (i in 1:m) {
    n = rpois(1, lam)
    s[i] = sum(rexp(n+1, rate))  }  # note the 'n+1'
 mean(s);  var(s)
 ## 17.93764  # approximates E(S) = (5 + 1)(3)
 ## 97.76707  # approximates V(S) = 6(9) + 5(9)
 mean(s > 20)
 ## 0.36613   # approximates P(S > 20)

You can $roughly$ verify your formula for $V(S)$ by comparing
it with the simulated value above. (Variances deal with squared
units, so numerical approximations can be less precise.) 
The histogram below shows the simulated distribution of $S$. The dotted density
curve of $Gamma(6, rate=1/3)$ would apply to the $deterministic$
sum of $n = 6$ random variables $X_i \sim Exp(rate=1/3)$. Using a
$random$ Poisson sum the simulated distribution of the histogram has a larger variance.
The simulated distribution is markedly right-skewed, so that
$P(S > 20)$ cannot be well approximated by the false assumption
that $S$ must be nearly normally distributed as $Norm(E(S)=18,
SD(S)=\sqrt{99}),$ solid green density curve.

Note: See a similar discussion at  http://math.stackexchange.com/questions/1483316/. 
