Thanks in advance.
Let $X$ a gamma-distributed random variable having scale $θ$ and shape $k$:
$$ X \sim \Gamma(k, \theta) \equiv \textrm{Gamma}(k, \theta) $$
with its probability density function is: $$ f(x;k,\theta) = \frac{x^{k-1}e^{-\frac{x}{\theta}}}{\theta^k\Gamma(k)} \quad \text{ for } x > 0 \text{ and } k, \theta > 0 $$
The sum of $N$ independent variables $X_i$ with Gamma distribution will be another Gamma distribution $$ \sum_{i=1}^N X_i \sim\mathrm{Gamma} \left( \sum_{i=1}^N k_i, \theta \right) $$
My question
What would be the distribution of $Z$?, being $Z$ $$ Z = \sum_{i=1}^N X_i $$ if $N$ is a gamma distributed random variable and $X$ is another gamma distribution random variable $$ X \sim \Gamma(k_1, \theta_1) \equiv \textrm{Gamma}(k_1, \theta_1) $$ $$ N \sim \Gamma(k_2, \theta_2) \equiv \textrm{Gamma}(k_2, \theta_2) $$
I have the numerical approach but I do not know how to get to an analytical solution