Finding joint distribution from pdf of sum Let $X$ and $Y$ be dependent random variables, where $X$ follow a Poisson distribution with rate $\lambda_x$. Let their sum $Z = X + Y$ follows a Poisson distribution with rate $\lambda_z$. 
Is this enough information to figure out the joint distribution of $(X, Y)$?
 A: The answer is: it is impossible in general to find joint distribution of $X$ and $Y$ given $X$ follows a Poisson distribution with rate $λ_x$ and $Z=X+Y$ follows a Poisson distribution with rate $λ_z$. Other words, one can construct different examples of joint distributions of $X,Y$ such that all the above conditions are satisfied. 
I assume here $\lambda_z>\lambda_x$.
1st example. Let $X\sim\Pi_{\lambda_x}$ and $Z\sim\Pi_{\lambda_z}$ be independent r.v.'s. Set $Y=Z-X$. Then $Z=X+(Z-X)=X+Y$.
In this case joint distribution of $(X,Y)=(X,Z-X)$ looks as follows:
$$
P(X=k, Z-X = m) = P(X=k, Z=m+k) = \dfrac{\lambda_x^k}{k!}e^{-\lambda_x}\times\dfrac{\lambda_z^{m+k}}{(m+k)!}e^{-\lambda_z},\ k\geq 0, \ m\geq -k.
$$
2nd example. For simplicity, I construct joint distribution for $\lambda_x=2$. Let $X\sim \Pi_2$ and construct $Y_0$ depending on $X$ s.t. $X+Y_0\sim\Pi_2$. After that we can choose another $Y_1\sim \Pi_{\lambda_z-\lambda_x}$ independent on $X$ and finally get $Y=Y_0+Y_1$. 
Let the joint pmf of $(X,Y_0)$ looks as follows:
$$\begin{array}{r|c|c|c|c|c}
X & 0 & 1 & 2 & 3 & \ldots \\
Y_0\quad & & & & & \\ 
\hline
0 & \frac12e^{-2} & \frac32e^{-2} & \frac{2^2}{2}e^{-2} & \frac{2^3}{3!}e^{-2} & \ldots \\
1 & \frac12e^{-2} & 0 & 0 & 0 & \ldots \\
-1 & 0 & \frac12e^{-2} & 0 & 0 & \ldots \\
\hline
\end{array}$$
Here 
$$
P(X=0,Y_0=0)= P(X=1,Y_0=-1)=P(X=0,Y_0=1) = \frac12e^{-2},
$$
$$
P(X=1,Y_0=0) = \frac32e^{-2}.
$$
Then 
$$
P(X+Y_0=0)=P(X=0,Y_0=0)+ P(X=1,Y_0=-1)=e^{-2}=P(X=0),$$
$$
P(X+Y_0=1)=P(X=0,Y_0=1)+ P(X=1,Y_0=0)=2e^{-2}=P(X=1).$$
Therefore, $X+Y_0$ is distributed as $X$. Next, choose another $Y_1\sim \Pi_{\lambda_z-2}$ independent on $X$ and $Y_0$ and finally get $Y=Y_0+Y_1$. 
Then $Z=X+Y=X+Y_0+Y_1\sim\Pi_{\lambda_z}$ is distributed like a sum of two independent Poisson r.v. $X+Y_0\sim\Pi_{2}$ and $Y_1\sim \Pi_{\lambda_z-2}$. 
Joint distribution of $(X,Y)=(X,Y_0+Y_1)$ looks as follows:
$$
P(X=k, Y=m) = P(X=k, Y_0=0)P(Y_1=m)+ P(X=k, Y_0=1)P(Y_1=m-1)+P(X=k, Y_0=-1) P(Y_1=m+1)$$
For example, when $k=0$, $m=1$
$$P(X=0, Y=1) = P(X=0, Y_0=0)P(Y_1=1)+ P(X=0, Y_0=1)P(Y_1=0)+P(X=0, Y_0=-1) P(Y_1=2)= \dfrac12(\lambda_z-2)e^{-\lambda_z}+\dfrac12e^{-\lambda_z}=\dfrac12(\lambda_z-1)e^{-\lambda_z}.
$$
This distribution differ from the distribution in 1st example.
