Possible methods for parameter estimation of a compound Poisson Let $X_i$ be iid and each $X_i$ takes one value among $(0,1,-1)$ with probs $(p_1,p_2,1-p_1-p_2)$ respectively. Let $N$ be a Poisson RV with mean $\lambda$, and $$Z=\sum_{i=0}^N X_i$$
be a compound Poisson accordingly. 
Question: Given $n$ independent observations $(Z_1,\cdots,Z_n)$ of $Z$, how can we make inference about the parameters $(p_1,p_2,\lambda)$? Is there any common technique for the parameter estimation of compound Possion process? (It seems that MLE works out poorly becasue the ML function entails infinite summation.)

Example plot of the histogram of many observations of $Z$

 A: Note first that those $X_i = 0$ contribute nothing to the sum.  We can exclude them by remembering that if we generate a Poisson variate $x$ with mean $\lambda$ and then a Binomial $z$ with probability parameter $p$ and size parameter equal to $x$, $z \sim \text{Poisson}(p\lambda)$; consequently the number of nonzero elements of our sum, label it $M$, $\sim \text{Poisson}((1-p_1)\lambda)$.
Now we can see that we have a problem; we can't distinguish between different values of $(p_1, \lambda)$ pairs, because we don't know how many of the $X_i = 0$ in each sample, and we have no way of finding out.   $\lambda = 1,000$ and $p = 0.99$ are indistinguishable in terms of the final distribution from $\lambda = 10$ and $p = 0$; the distribution of $M$, and therefore the distribution of the sum itself, is the same in either case.  Thus, the parameters are not identifiable on the basis of the collected sample statistics (although if we also collected the number of zero observations in each sample, they would be.)
If we alter the problem to the tractable one of $X_i \in \{-1,1\}$ with probabilities $(p, 1-p)$ instead of the tri-valued $X_i$ in the original problem, then we can apply a method of moments estimation based upon the first two moments, or maximum likelihood, or something else.
The MOM estimator is based on equating the sample mean and variance to the population moments; in this case, $\mathbb{E}(x) = \lambda(2p-1)$, and $\sigma^2(x)$ can be found by noting that the variance of a single observation of $X$ equals $4p(1-p)$, therefore the variance of the sum of $M$ observations is $4p(1-p)M$, and a little more math gets us to $\sigma^2(x) = 4p(1-p)\lambda + (2p-1)^2\lambda = \lambda$.  This latter result is surprising, at least it surprised me, so after triple-checking the math ($p(\text{correct})$ still $< 1$, of course), I simulated it:
x <- rep(0, 1000000)
p <- 0.3
lambda <- 5
for (j in 1:length(x)) {
  M <- rpois(1,lambda)
  z <- rbinom(1,M,p)
  x[j]  <- z - (M-z)
}
> mean(x)
[1] -1.997802
> var(x)
[1] 5.008178

At this point how to calculate the MOM estimate should be clear.
Note that $\lambda$ in this problem is equivalent to $(1-p_1)\lambda$ in the original problem; we can estimate the original problem's $(1-p_1)\lambda$, but not the individual components thereof.  
A: I have no experience on this field, but what about method of moments?
We could work it out:
\begin{align*}
\mathbb{E}\left(Z\right)&=\mathbb{E}\left[\mathbb{E}\left(Z\mid N\right)\right]=\mathbb{E}\left[N\cdot\left(p_2-1+p_1+p_2\right)\right]=\lambda \left(2p_2+p_1-1\right)\\
\mathbb{E}\left(Z^2\right)&=\mathbb{E}\left[\mathbb{E}\left(Z^2\mid N\right)\right]=\mathbb{E}\left[\sum_{i=1}^N\sum_{j=1}^N \mathbb{E}\left(X_i X_j\right)\right]=\\
&=\mathbb{E}\left[N\left(1-p_1\right)+N\left(N-1\right)\left(p_2^2+\left(1-p_1-p_2\right)^2-2p_2\left(1-p_1-p_2\right)\right)\right]=\\
&=\lambda\left(1-p_1\right)+\lambda^2\left(2p_2+p_1-1\right)^2\\
\mathbb{E}\left(Z^3\right)&=\mathbb{E}\left[\mathbb{E}\left(Z^3\mid N\right)\right]=\mathbb{E}\left[\sum_{i=1}^N\sum_{j=1}^N \sum_{k=1}^N\mathbb{E}\left(X_i X_jX_k\right)\right]=\\
&=\mathbb{E}\left[N\left(p_2-1+p_1+p_2\right)+3N\left(N-1\right)\left(\left(1-p_1\right)p_2-\left(1-p_1\right)\left(1-p_1-p_2\right)\right)\right.+\\
&+\left.N\left(N-1\right)\left(N-2\right)\left(p_2^3+3\left(1-p_1-p_2\right)^2p_2-\left(1-p_1-p_2\right)^3\right.\right.-\\
&\left.\left.-3\left(1-p_1-p_2\right)p_2^2\right)\right]=\\
&=\lambda\left(2p_2+p_1-1\right)+3\lambda^2\left(1-p_1\right)\left(2p_2+p_1-1\right)+\lambda^3\left(2p_2+p_1-1\right)^3
\end{align*}
This system has no solution, as the third equation is just the linear combination of the first two. Perhaps try higher moments?
Manual calculation of higher moments gets very tedious in the above way, but moment generating function comes to the rescue. The MGF of $X$ is $p_1+p_2e^t+\left(1-p_1-p_2\right)e^{-t}$, so the MGF of $Z$ is $$M\left(t\right)=e^{\lambda\left(p_1+p_2e^t+\left(1-p_1-p_2\right)e^{-t}-1\right)}.$$
Now, the problem becomes apparent: derivatives of $M\left(t\right)$ can be calculated recursively from the first two terms, i.e. every moment will be a linear combination of the first two moments! (As moments are the respective derivatives of the MGF at $t=0$; you can easily check that it works for this three as well.) But we have three parameters!
That is: it is impossible to solve this problem with the method of moments.
...and this way, I think we have a very-very-very-very complicated and cumbersome proof for what @jbowman solved in a single paragraph.
