Assume $X, Y$ two independent exponential random variables with $\lambda = 1$, and $Z$ to be a gamma variate with $\alpha = \theta > 0, \beta = 1$, independent from $X$ and $Y$.
Next, assume $U = [1 + (X/Z)]^{-\theta}$, $V = [1 + (Y/Z)]^{-\theta}$. I'm trying to show that $$\mathbb{P}(U \leq u, V \leq v)=[\text{max}(u^{-\theta}+v^{-\theta}-1,0)]^{-1/\theta}$$
By replacing $U, V$ with their expressions we get (or not? see edit):
$$\mathbb{P}(U \leq u, V \leq v)=\mathbb{P}(X \leq Z[u^{-1/\theta} - 1], Y \leq Z[v^{-1/\theta} - 1])$$
Which in turn reduces to:
$$\int_{0}^{\infty}\mathbb{P}(X \leq Z[u^{-1/\theta} - 1], Y \leq Z[v^{-1/\theta} - 1]|Z=z)\mathbb{P}(Z=z)dz$$
Here by plugging-in the known distributions, we get
$$\frac{1}{\Gamma(\theta)}\int_{0}^{\infty}[1-e^{-z(u^{-1/\theta}-1)}][1-e^{-z(v^{-1/\theta}-1)}]z^{\theta-1}e^{-z}dz$$
Assuming I'm correct at evaluating the, i.e.,
$$ \int_{0}^{\infty}e^{-z(u^{-1/\theta}-1)}z^{\theta-1}e^{-z}dz = \int_{0}^{\infty}e^{-z(u^{-1/\theta})}z^{\theta-1}dz = [u^{-1/\theta}]^{-\theta} \Gamma(\theta) = u \Gamma(\theta), \tag{1}$$
By unbracing the integral I get the final result to be $ 1 - u - v + (u^{-1/\theta} + v^{-1/\theta} -1)^{-\theta} $, which is not exactly the expected result.
Where am I wrong? I'm not entirely sure on the $(1)$ result, also I might have messed something up in deriving the probability, so any hints would be appreciated!
EDIT: Or does raising a power in probability change the sign? That is, we get the probability at
$$\mathbb{P}(U \leq u, V \leq v)=\mathbb{P}(X > Z[u^{-1/\theta} - 1], Y > Z[v^{-1/\theta} - 1])$$In this case the final integral would result in $(u^{-1/\theta} + v^{-1/\theta} -1)^{-\theta} $, which, however, is still not the same result. Unless we're free here to set $\theta := 1/\theta$..