PDF, mean and variance of product of two dependent random variables Consider two dependent random variables 
$$X \sim \text{Ga}(\alpha,\beta)
\quad \quad \quad \quad
Y|X \sim \text{N} \Big( 0,\frac{1}{X} \Big).$$
How do I compute the PDF, mean and variance of the random variable $Z = X \times Y$?
 A: The moments are obtained using iterated moment formulae.  Using the law of iterated expectation you have:
$$\begin{equation} \begin{aligned}
\mathbb{E}(Z) 
&= \mathbb{E}(\mathbb{E}(Z|X)) \\[6pt]
&= \mathbb{E}(\mathbb{E}(X \cdot Y|X)) \\[6pt]
&= \mathbb{E}(X \cdot \mathbb{E}(Y|X)) \\[6pt]
&= \mathbb{E}(X \cdot 0) \\[6pt]
&= \mathbb{E}(0) = 0. \\[6pt]
\end{aligned} \end{equation}$$
Using the law of iterated variance you have:
$$\begin{equation} \begin{aligned}
\mathbb{V}(Z) 
&= \mathbb{V}(\mathbb{E}(Z|X)) + \mathbb{E}(\mathbb{V}(Z|X)) \\[6pt]
&= \mathbb{V}(\mathbb{E}(X \cdot Y|X)) + \mathbb{E}(\mathbb{V}(X \cdot Y|X)) \\[6pt]
&= \mathbb{V}(X \cdot \mathbb{E}(Y|X)) + \mathbb{E}(X^2 \cdot \mathbb{V}(Y|X)) \\[6pt]
&= \mathbb{V}(X \cdot 0) + \mathbb{E} \Big( X^2 \cdot \frac{1}{X} \Big) \\[6pt]
&= \mathbb{V}(0) + \mathbb{E}(X) \\[6pt]
&= \frac{\alpha}{\beta}. \\[6pt]
\end{aligned} \end{equation}$$
To find the PDF of $Z$ we can start by observing that this distribution is symmetric around the point $z=0$, so we can proceed for the positive case and then use symmetry for the rest.  To find the distribution function you can use the law of total probability for any $z \geqslant 0$ to get:
$$\begin{equation} \begin{aligned}
F_Z(z) 
= \mathbb{P}(Z \leqslant z) 
&= \mathbb{P}(X \cdot Y \leqslant z) \\[6pt]
&= \int \limits_0^\infty \mathbb{P}(X \cdot Y \leqslant z|X=x) \cdot f_X(x) \ dx \\[6pt]
&= \int \limits_0^\infty \mathbb{P} \Big( Y \leqslant \frac{z}{x} \Big| X=x \Big) \cdot f_X(x) \ dx \\[6pt]
&= \int \limits_0^\infty \Phi \Big( \frac{z}{\sqrt{x}} \Big) \cdot \text{Ga}(x|\alpha, \beta) \ dx. \\[6pt]
\end{aligned} \end{equation}$$
For all $z \geqslant 0$ the density corresponding to this distribution function is:
$$\begin{equation} \begin{aligned}
f_Z(z) 
= \frac{dF_Z}{dz}(z) 
&= \frac{d}{dz} \int \limits_0^\infty \Phi \Big( \frac{z}{\sqrt{x}} \Big) \cdot \text{Ga}(x|\alpha, \beta) \ dx \\[6pt]
&= \int \limits_0^\infty \frac{\partial}{\partial z}  \Phi \Big( \frac{z}{\sqrt{x}} \Big) \cdot \text{Ga}(x|\alpha, \beta) \ dx \\[6pt]
&= \int \limits_0^\infty \frac{1}{\sqrt{x}} \cdot \phi \Big( \frac{z}{\sqrt{x}} \Big) \cdot \text{Ga}(x|\alpha, \beta) \ dx \\[6pt]
&= \int \limits_0^\infty x^{-1/2} \cdot \frac{1}{\sqrt{2 \pi}} \exp \Big( -\frac{1}{2} \frac{z^2}{x} \Big) \cdot \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} \exp(- \beta x) \ dx \\[6pt]
&= \frac{1}{\sqrt{2 \pi}} \frac{\beta^\alpha}{\Gamma(\alpha)} \int \limits_0^\infty x^{\alpha -3/2} \exp \Big( - \beta x -\frac{1}{2} \frac{z^2}{x} \Big) \ dx. \\[6pt]
\end{aligned} \end{equation}$$
Using symmetry, we then have:
$$f_Z(z) 
= \frac{1}{\sqrt{2 \pi}} \frac{\beta^\alpha}{\Gamma(\alpha)} \int \limits_0^\infty x^{\alpha -3/2} \exp \Big( - \beta x -\frac{1}{2} \frac{z^2}{x} \Big) \ dx
\quad \quad \quad \text{for } z \in \mathbb{R}.$$
There is no closed-form expression for this integral, so the density cannot be simplified any further.  The density can be computed numerically using standard methods of numerical computation of integrals.
A: Use the multivariate change of variable formula.
Let $Z = XY$ and $V = Y$ so that $Z/V = X$.
Then we have the multivariate change of variable formula, as shown here:
$f_{Z,V}(z,v) = f_{X,Y}(x = \frac{Z}{V}, y = V) * |J|$ where $|J|$ is the Jacobian.
Now, before we do that, we need to figure out what $f_{X,Y}(x,y)$ is. 
First, note that the density of Y is actually a conditional density, that is $Y|X \sim N(0, \frac{1}{X})$
Since we know $f_{Y|X}(y) = \frac{f_{X,Y}(x,y)}{f_{X}(x)}$ we can simply multiply $f_{Y|X}(y)$ (which is given to us) by $f_{X}(x)$ (which is also given as a Gamma($\alpha, \beta$)) to obtain the joint density. Then, using the above steps (which I am assuming you probably know how to do) you can obtain the joint density of Z and V. 
Then, integrate the resulting joint density  over all values of V to get the marginal of Z, $f_Z(z)$, from which you can then obtain the desired moments needed for the mean and variance of Z.
Hope this helps. 
