What is the distribution of a random variable whose pmf is a function of other pmfs?

I'm interested to find out what is the resulting distribution for a discrete random variable $$Y$$ whose probability mass function is defined as:

$$Pr(Y = k) = C \cdot \frac{Pr(X_1 = k) \cdot Pr(X_2 = k)}{Pr(X_3 = k)}$$ for $$k = 0, 1, 2,\ldots$$, where each of the (independent) $$X_i$$'s follows a Negative Binomial distribution with $$X_i \sim NB(r_i, p_i)$$ for $$i = 1, 2, 3$$, while $$C$$ is the normalizing constant.

I've managed to show that when $$X_i \sim Poisson(λ_i)$$ then $$Y \sim Poisson(\frac{λ_1 λ_2}{λ_3})$$ but have so far failed to get somewhere with the negative binomial case.

Starting from:

$$Pr(X_i = k) = {k + r_i - 1 \choose k} {p_i}^{r_i} (1 - p_i)^k$$

I've managed to find (I think!) that the pmf of $$Y$$ would include this quantity: $$\frac{{p_1}^{r_1} {p_2}^{r_2}}{{p_3}^{r_3}} \cdot \bigg[\frac{(1-p_1)(1-p_2)}{(1-p_3)}\bigg]^k$$ and other factorials / Gamma functions involving $$r_1, r_2, r_3$$ and $$k$$, but couldn't reach the end. My gut feeling (backed by some simulations) is that $$Y$$ will also follow a negative binomial distribution but can't prove if that's the case and, if so, what the corresponding expressions for $$r_y$$ and $$p_y$$ (in terms of the $$r_i$$'s and the $$p_i$$'s) will be.

Any help will be much appreciated. Thanks!

On the face of it it seems your suggestion could work for the negative binomial if $$r_1=r_2=r_3$$ even if $$p_1,p_2,p_3$$ differ so long as $$p_3$$ is small enough, but might not work if $$r_1,r_2,r_3$$ differ.

Suppose $$r_1=r_2=r_3=r$$. Then

$$\Pr(Y = k) = C \cdot \frac{\Pr(X_1 = k) \cdot \Pr(X_2 = k)}{\Pr(X_3 = k)} \\ = C \cdot \frac{{k + r - 1 \choose k} {p_1}^{r} (1 - p_1)^k \cdot {k + r - 1 \choose k} {p_2}^{r} (1 - p_2)^k}{{k + r - 1 \choose k} {p_3}^{r} (1 - p_3)^k} \\ = {k + r - 1 \choose k} C \cdot \frac{{p_1}^{r} \cdot {p_2}^{r} }{ {p_3}^{r} }\left(\frac{(1 - p_1) \cdot (1 - p_2)}{ (1 - p_3)}\right)^k$$

and now let $$p=1-\frac{(1 - p_1) \cdot (1 - p_2)}{ (1 - p_3)}$$ and $$C=\frac{ {p}^{r} \cdot{p_3}^{r} }{{p_1}^{r} \cdot {p_2}^{r} }$$ so $$\Pr(Y = k) = {k + r - 1 \choose k} p^r \left(1-p\right)^k$$ which is a negative binomial, at least when $$0 < p < 1$$ which requires $$1-p_3 > (1-p_1)\cdot (1-p_2)$$; you need that condition anyway as otherwise $$\frac{\Pr(X_1 = k) \cdot \Pr(X_2 = k)}{\Pr(X_3 = k)}$$ would be an increasing function of $$k$$ and so have an infinite sum making $$C$$ zero.

Now for a counterexample where $$r_1,r_2,r_3$$ differ. Suppose $$r_1=r_2=1$$ but $$r_3=2$$, and $$p_1=p_2=p_3=\frac12$$. Then

$$\Pr(Y = k) = C \cdot \frac{{k \choose k} \frac12 \frac1{2^{k}}\cdot {k \choose k} \frac12 \frac1{2^{k}}}{{k+1 \choose k} \frac1{2^2} \frac1{2^{k}}} = C \cdot \frac{ 1}{(k+1){2^{k}} }$$ which will require $$C=\frac{1}{2\log_e(2)}$$. This $$Y$$ will have mean $$\frac{1}{\log_e(2)}-1\approx 0.4427$$ and variance $$\frac{2}{\log_e(2)}- \frac{1}{\log_e(2)^2}\approx 0.8040$$ and $$\Pr(Y = 0) = C \approx 0.7213$$.

If this $$Y$$ were negative binomial then this mean and variance would suggest $$p = \frac{\log_e(2)\left(1-\log_e(2)\right)}{2\log_e(2)-1}\approx 0.5506$$ and $$r=\frac{(1-\log_e(2))^2}{\log_e(2)^2 +\log_e(2)-1} \approx 0.5424$$. That $$r$$ is not an integer, but even if we allowed that, those $$p$$ and $$r$$ would suggest $$\Pr(Y = 0)=p^r \approx 0.7235$$ which is wrong, so we can conclude that this $$Y$$ does not have a negative binomial distribution