If $X$ takes on a Gamma Distribution, how can I find $X^2$, $X^3$, etc? I am trying to take consecutive powers of a Gamma Distribution. For example, if 
$X \sim  \textrm{Gamma}(k, \theta)$, 
I would like to find $X^2$, $X^3$, and in general $X^m$ for $m>0$. 
The pdf using the shape-scale parametrization 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$.
I have started to do this by working by the following:
$P(X^m \leq x) = P(-\sqrt[m]{x} \leq X \leq \sqrt[m]{x}) = \int_{-\sqrt[m]{x}}^{\sqrt[m]{x}} \frac{t^{k-1}e^{-\frac{t}{\theta}}}{\theta^k\Gamma(k)} \,dt$.
Then, $\quad f_{X^m}(x) = \frac{d}{dx}\int_{-\sqrt[m]{x}}^{\sqrt[m]{x}} \frac{t^{k-1}e^{-\frac{t}{\theta}}}{\theta^k\Gamma(k)} \,dt$ $\quad$
which by the Fundamental Theorem of Calculus equates to:
$\dfrac{(\sqrt[m]{x})^{k-1}e^{-\dfrac{\sqrt[m]{x}}{\theta}}}{\theta^k\Gamma(k)}\cdot(\sqrt[m]{x})'-\dfrac{(-\sqrt[m]{x})^{k-1}e^{-\dfrac{-\sqrt[m]{x}}{\theta}}}{\theta^k\Gamma(k)}\cdot(-\sqrt[m]{x})'$.
After this, I am stuck because there are many different cases to consider, such as when $m$ is an integer, etc. 
Could anyone tell me if my strategy is correct? Thanks!
 A: Your strategy will fail.  Powers of $X$ follow Generalized Gamma distributions.  Their densities are given by
$$f(x; k, \gamma, 1) = \frac{1}{\Gamma (k)} x^{\gamma  k}  e^{-x^{\gamma }} \left(\gamma \frac{dx}{x} \right).$$
(The third parameter, here set to a unit value $\sigma=1$, is a scale parameter.)  For the power $m$ of $X \sim \Gamma(k)$, the shape parameter is $\gamma = 1/m$. This all becomes obvious when you think of $x^\gamma$ as being another variable, say $u$, and note that 
$$\frac{du}{u} = \frac{d(x^\gamma)}{x^\gamma} = \gamma \frac{d x}{x}.$$
The remaining factors, $u^k e^{-u}/\Gamma(k)$, give the PDF of a Gamma distribution (with respect to the measure $du/u$).
The issue before us is whether for any given $k$ it is possible to find other parameters $k^\prime$ and $\sigma^\prime$ such that $f(x,k, \gamma, \sigma)=f(x;k^\prime, 1, \sigma^\prime)$, where $\gamma = 1/m$ for some integral $m \gt 1$: that is what it would mean for  $X^m$ to have an (ordinary) Gamma distribution.
We can settle this by looking at moments.  Integration shows that the (non-central) moments of the generalized Gamma distribution are
$$\mu_i = \sigma^i\frac{\Gamma(k+m i)}{\Gamma(k)}.$$ 
For integral values of $m$ their successive ratios are
$$\frac{\mu_{i+1}}{\mu_i} = \frac{\sigma^{i+1}\Gamma(k + mi + m)}{\sigma^i\Gamma(k + m i)} = \sigma(k + mi)(k+mi+1)\cdots (k+mi+m-1).$$
Trying to solve equations (or rather, prove they have no solutions) involving these formulas looks messy.  Instead, consider the limiting values 
$$\lim_{i\to \infty}\frac{\mu_{i+1}}{i\mu_i}.$$
For $m=1$ (a true Gamma distribution) this is the limiting value of $\sigma(k+i)/i$, obviously equal to $\sigma$, whereas for $m\gt 1$ this is the limiting value of $\sigma(k+mi)\cdots(k+mi+m-1)/i$, which (being of order $i^{m-1}$) diverges.  Therefore the distribution of $X^m$ for $m\gt 1$ (and integral) cannot have the same moments as any possible Gamma function, proving it is not a Gamma function.

It should be clear that the generalized Gamma distribution family is indeed closed under taking (positive) powers.
