I've recently encountered the bivariate Poisson distribution, but I'm a little confused as to how it can be derived.
The distribution is given by:
$P(X = x, Y = y) = e^{-(\theta_{1}+\theta_{2}+\theta_{0})} \displaystyle\frac{\theta_{1}^{x}}{x!}\frac{\theta_{2}^{y}}{y!} \sum_{i=0}^{min(x,y)}\binom{x}{i}\binom{y}{i}i!\left(\frac{\theta_{0}}{\theta_{1}\theta_{2}}\right)^{i}$
From what I can gather, the $\theta_{0}$ term is a measure of correlation between $X$ and $Y$; hence, when $X$ and $Y$ are independent, $\theta_{0} = 0$ and the distribution simply becomes the product of two univariate Poisson distributions.
Bearing this in mind, my confusion is predicated on the summation term - I'm assuming this term explains the the correlation between $X$ and $Y$.
It seems to me that the summand constitutes some sort of product of binomial cumulative distribution functions where the probability of "success" is given by $\left(\frac{\theta_{0}}{\theta_{1}\theta_{2}}\right)$ and the probability of "failure" is given by $i!^{\frac{1}{min(x,y)-i}}$, because $\left(i!^{\frac{1}{min(x,y)-i!}}\right)^{(min(x,y)-i)} = i!$, but I could be way off with this.
Could somebody provide some assistance on how this distribution can be derived? Also, if it could be included in any answer how this model might be extended to a multivariate scenario (say three or more random variables), that would be great!
(Finally, I have noted that there was a similar question posted before (Understanding the bivariate Poisson distribution), but the derivation wasn't actually explored.)