Distribution of $XY$ if $X \sim$ Beta$(1,K-1)$ and $Y \sim$ chi-squared with $2K$ degrees Suppose that $X$ has the beta distribution Beta$(1,K-1)$ and $Y$  follows a chi-squared with $2K$ degrees. In addition, we assume that $X$ and $Y$ are independent.
What is the distribution of the product $Z=XY$ .
Update
My attempt:
\begin{align}
f_Z &= \int_{y=-\infty}^{y=+\infty}\frac{1}{|y|}f_Y(y) f_X \left (\frac{z}{y} \right ) dy \\ &= \int_{0}^{+\infty}  \frac{1}{B(1,K-1)2^K \Gamma(K)}  \frac{1}{y} y^{K-1} e^{-y/2} (1-z/y)^{K-2} dy \\ &= \frac{1}{B(1,K-1)2^K \Gamma(K)}\int_{0}^{+\infty} e^{-y/2} (y-z)^{K-2} dy \\ &=\frac{1}{B(1,K-1)2^K \Gamma(K)} [-2^{K-1}e^{-z/2}\Gamma(K-1,\frac{y-z}{2})]_0^\infty \\ &= \frac{2^{K-1}}{B(1,K-1)2^K \Gamma(K)} e^{-z/2} \Gamma(K-1,-z/2)
\end{align}
Is it correct? if yes, how we call this distribution? 
 A: After some valuable remarks, I was able to find the solution:
We have $f_X(x)=\frac{1}{B(1,K-1)} (1-x)^{K-2}$ and $f_Y(y)=\frac{1}{2^K \Gamma(K)} y^{K-1} e^{-y/2}$.
Also, we have $0\le x\le 1$. Thus, if $x=\frac{z}{y}$, we get  $0 \le \frac{z}{y}  \le 1$ which implies that $z\le y \le \infty$.
Hence:
\begin{align}
f_Z &= \int_{y=-\infty}^{y=+\infty}\frac{1}{|y|}f_Y(y) f_X \left (\frac{z}{y} \right ) dy \\ &= \int_{z}^{+\infty}  \frac{1}{B(1,K-1)2^K \Gamma(K)}  \frac{1}{y} y^{K-1} e^{-y/2} (1-z/y)^{K-2} dy \\ &= \frac{1}{B(1,K-1)2^K \Gamma(K)}\int_{z}^{+\infty} e^{-y/2} (y-z)^{K-2} dy \\ &=\frac{1}{B(1,K-1)2^K \Gamma(K)} \left[-2^{K-1}e^{-z/2}\Gamma(K-1,\frac{y-z}{2})\right]_z^\infty \\ &= \frac{2^{K-1}}{B(1,K-1)2^K \Gamma(K)} e^{-z/2} \Gamma(K-1) \\ &= \frac{1}{2} e^{-z/2}
\end{align}
where the last equality holds since $B(1,K-1)=\frac{\Gamma(1)\Gamma(K-1)}{\Gamma(K)}$.
So $Z$ follows an exponential distribution of parameter $\frac{1}{2}$; or equivalently, $Z \sim\chi_2^2$.
A: There is a pleasant, natural statistical solution to this problem for integral values of $K$, showing that the product has a $\chi^2(2)$ distribution.  It relies only on well-known, easily established relationships among functions of standard normal variables.
When $K$ is integral, a Beta$(1,K-1)$ distribution arises as the ratio $$\frac{X}{X+Z}$$ where $X$ and $Z$ are independent, $X$ has a $\chi^2(2)$ distribution, and $Z$ has a $\chi^2(2K-2)$ distribution.  (See the Wikipedia article on the Beta distribution for instance.)  
Any $\chi^2(n)$ distribution is that of the sum of squares of $n$ independent standard Normal variates.  Consequently, $X+Z$ is distributed as the squared length of a $2 + 2K-2 = 2K$ vector with a standard multinormal distribution in $\mathbb{R}^{2K}$ and $X/(X+Z)$ is the squared length of the first two components when that vector is radially projected to the unit sphere $S^{2K-1}$.
The projection of a standard multinormal $n$-vector onto the unit sphere has a uniform distribution because the multinormal distribution is spherically symmetric. (That is, it is invariant under the orthogonal group, a result that follows immediately from two simple facts: (a), the orthogonal group fixes the origin and by definition does not change covariances; and (b) the mean and covariance completely determine the multivariate normal distribution. I illustrated this for the case $n=3$ at https://stats.stackexchange.com/a/7984).  In fact, the spherical symmetry immediately shows this distribution is uniform conditional on the length of the original vector.  The ratio $X/(X+Z)$ therefore is independent of the length.
What all this implies is that multiplying $X/(X+Z)$ by an independent $\chi^2(2K)$ variable $Y$ creates a variable with the same distribution as $X/(X+Z)$ multiplied by $X+Z$; to wit, the distribution of $X$, which has a $\chi^2(2)$ distribution.
A: I greatly deprecate the commonly used tactic of finding the density of $Z = g(X,Y)$ by computing first computing the joint density of $Z$ and
$X$ (or $Y$) because it is "easy" to use Jacobians, and then getting
$f_Z$ as a marginal density (cf. Rusty Statistician's answer). It is much easier to find the CDF of $Z$ directly and then differentiate
to find the pdf. This is the approach used below. 
$X$ and $Y$ are independent random variables with densities
$f_X(x) = (K-1)(1-x)^{K-2}\mathbf 1_{(0,1)}(x)$ and
$f_Y(y) = \frac{1}{2^K (K-1)!} y^{K-1} e^{-y/2}\mathbf 1_{(0,\infty)}(y)$.
Then, with $Z = XY$, we have for $z > 0$,
\begin{align}
P\{Z > z\} &= P\{XY > z\}\\
&= \int_{y=z}^\infty \frac{1}{2^K (K-1)!} y^{K-1} e^{-y/2}
\left[\int_{x=\frac{z}{y}}^1 (K-1)(1-x)^{K-2}\,\mathrm dx \right]
\,\mathrm dy\\
&= \int_{y=z}^\infty \frac{1}{2^K (K-1)!} y^{K-1} e^{-y/2}
\left(1-\frac{z}{y}\right)^{K-1}\,\mathrm dy\\
&= \int_{y=z}^\infty \frac{1}{2^K (K-1)!} (y-z)^{K-1} e^{-y/2}
\,\mathrm dy\\
&= e^{-z/2}\int_0^\infty \frac{1}{2^K (K-1)!}t^{K-1} e^{-t/2}
\,\mathrm dy~~~\scriptstyle{\text{on setting}~y-z = t}\\
&= e^{-z/2}\qquad\scriptstyle{\text{on noting that the integral is that
of a Gamma pdf}}\\
\end{align}
It is well-known that if $V \sim \mathsf{Exponential}(\lambda)$, then $P\{V > v\} = e^{-\lambda v}$. It follows that $Z = XY$ has an
exponential density with parameter $\lambda = \frac 12$,
which is also the $\chi^2(2)$ distribution.
