Given a known Gaussian distribution, $X \sim \mathcal N(\mu_x, \sigma_x^2)$, how does one determine the the distribution of $Y$ if $Y = X^2$?
Normal Distribution – PDF of the Square of a General Normal Random Variable
density functiondistributionsnormal distributionself-study
Related Solutions
simplify the term in the integral to
$T=e^{-\frac{1}{2}((\frac{\frac{z}{y}-\mu_x}{\sigma_x} )^2 -y)} y^{k/2-2} $
find the polynomial $p(y)$ such that
$[p(y)e^{-\frac{1}{2}((\frac{\frac{z}{y}-\mu_x}{\sigma_x} )^2 -y)}]'=p'(y)e^{-\frac{1}{2}((\frac{\frac{z}{y}-\mu_x}{\sigma_x} )^2 -y)} + p(y) [-\frac{1}{2}((\frac{\frac{z}{y}-\mu_x}{\sigma_x} )^2 -y)]' e^{-\frac{1}{2}((\frac{\frac{z}{y}-\mu_x}{\sigma_x} )^2 -y)} = T$
which reduces to finding $p(y)$ such that
$p'(y) + p(y) [-\frac{1}{2}((\frac{\frac{z}{y}-\mu_x}{\sigma_x} )^2 -y)]' = y^{k/2-2}$
or
$p'(y) -\frac{1}{2} p(y) (\frac{z \mu_x }{\sigma_x^2} y^{-2} \frac{z^2}{\sigma_x^2} y^{-3} -1)= y^{k/2-2}$
which can be done evaluating all powers of $y$ seperately
edit after comments
Above solution won't work as it diverges.
Yet, some others have worked on this type of product.
Using Fourrier transform:
Schoenecker, Steven, and Tod Luginbuhl. "Characteristic Functions of the Product of Two Gaussian Random Variables and the Product of a Gaussian and a Gamma Random Variable." IEEE Signal Processing Letters 23.5 (2016): 644-647. http://ieeexplore.ieee.org/document/7425177/#full-text-section
For the product $Z=XY$ with $X \sim \mathcal{N}(0,1)$ and $Y \sim \Gamma(\alpha,\beta)$ they obtained the characteristic function:
$\varphi_{Z} = \frac{1}{\beta^\alpha }\vert t \vert^{-\alpha} exp \left( \frac{1}{4\beta^2t^2} \right) D_{-\alpha} \left( \frac{1}{\beta \vert t \vert } \right)$
with $D_\alpha$ Whittaker's function ( http://people.math.sfu.ca/~cbm/aands/page_686.htm )
Using Mellin transform:
Springer and Thomson have described more generally the evaluation of products of beta, gamma and Gaussian distributed random variables.
Springer, M. D., and W. E. Thompson. "The distribution of products of beta, gamma and Gaussian random variables." SIAM Journal on Applied Mathematics 18.4 (1970): 721-737. http://epubs.siam.org/doi/10.1137/0118065
They use the Mellin integral transform. The Mellin transform of $Z$ is the product of the Mellin transforms of $X$ and $Y$ (see http://epubs.siam.org/doi/10.1137/0118065 or https://projecteuclid.org/euclid.aoms/1177730201). In the studied cases of products the reverse transform of this product can be expressed as a Meijer G-function for which they also provide and prove computational methods.
They did not analyze the product of a Gaussian and gamma distributed variable, although you might be able to use the same techniques. If I try to do this quickly then I believe it should be possible to obtain an H-function (https://en.wikipedia.org/wiki/Fox_H-function ) although I do not directly see the possibility to get a G-function or make other simplifications.
$M\lbrace f_Y(x) \vert s \rbrace = 2^{s-1} \Gamma(\tfrac{1}{2}k+s-1)/\Gamma(\tfrac{1}{2}k)$
and
$M\lbrace f_X(x) \vert s \rbrace = \frac{1}{\pi}2^{(s-1)/2} \sigma^{s-1} \Gamma(s/2) $
you get
$M\lbrace f_Z(x) \vert s \rbrace = \frac{1}{\pi}2^{\frac{3}{2}(s-1)} \sigma^{s-1} \Gamma(s/2) \Gamma(\tfrac{1}{2}k+s-1)/\Gamma(\tfrac{1}{2}k) $
and the distribution of $Z$ is:
$f_Z(y) = \frac{1}{2 \pi i} \int_{c-i \infty}^{c+i \infty} y^{-s} M\lbrace f_Z(x) \vert s \rbrace ds $
which looks to me (after a change of variables to eliminate the $2^{\frac{3}{2}(s-1)}$ term) as at least a H-function
what is still left is the puzzle to express this inverse Mellin transform as a G function. The occurrence of both $s$ and $s/2$ complicates this. In the separate case for a product of only Gaussian distributed variables the $s/2$ could be transformed into $s$ by substituting the variable $x=w^2$. But because of the terms of the chi-square distribution this does not work anymore. Maybe this is the reason why nobody has provided a solution for this case.
Following up on whuber's comment, $$\begin{align} \left.\left.\frac{1}{2(1-\rho^2)}\right(x^2+y^2-2\rho xy\right) &=\left.\left.\frac{1}{2(1-\rho^2)}\right(x^2-\rho^2x^2+(y^2-2\rho xy + \rho^2x^2)\right)\\ &= \frac{x^2}{2} + \frac{1}{2}\left(\frac{y-\rho x}{\sqrt{1-\rho^2}}\right)^2 \end{align}$$ and so $$\int_{-\infty}^\infty f_{X,Y}(x,y)\,\mathrm dy = \frac{e^{-x^2/2}}{\sqrt{2\pi}} \int_{-\infty}^\infty \frac{e^{-(y-\rho x)^2/2(1-\rho^2)}}{\sqrt{1-\rho^2}\sqrt{2\pi}}\,\mathrm dy.$$ I will leave it to the OP to complete the details and determine whether $\rho$ disappears or not when the integral is evaluated.
Best Answer
You have stumbled upon one of the most famous results of probability theory and statistics. I'll write an answer, although I am certain this question has been asked (and answered) before on this site.
First, note that the pdf of $Y = X^2$ cannot be the same as that of $X$ as $Y$ will be nonnegative. To derive the distribution of $Y$ we can use three methods, namely the mgf technique, the cdf technique and the density transformation technique. Let's begin.
Moment generating function technique.
Or characteristic function technique, whatever you like. We have to find the mgf of $Y=X^2$. So we need to compute the expectation
$$E\left[e^{tX^2}\right]$$
Using the Law of the Unconscious Statistician, all we have to do is compute this integral over the distribution of $X$. Thus we need to compute
$$\begin{align} E\left[e^{tX^2}\right] = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi}} e^{tx^2} e^{-\frac{x^2}{2}} dx &= \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} \exp\left\{ - \frac{x^2}{2} \left( 1- 2t \right) \right\} dt \\ & = \int_{-\infty}^{\infty} \frac{\left( 1-2t \right)^{1/2}}{\left( 1-2t \right)^{1/2}} \frac{1}{\sqrt{2\pi}} \exp\left\{ - \frac{x^2}{2} \left( 1- 2t \right) \right\} dt \\ & = \left(1-2t\right)^{-1/2} , \quad t<\frac{1}{2} \end{align}$$
where in the last line we have compared the integral with a Gaussian integral with mean zero and variance $\frac{1}{\left(1-2t\right)}$. Of course this integrates to one over the real line. What can you do with that result now? Well, you may apply a very complex inverse transformation and determine the pdf that corresponds to this MGF or you may simply recognise it as the MGF of a chi-squared distribution with one degree of freedom. (Recall that a chi-squared distribution is a special case of a gamma distribution with $\alpha = \frac{r}{2}$, $r$ being the degrees of freedom, and $\beta = 2$).
CDF technique
This is perhaps the easiest thing you can do and it is suggested by Glen_b in the comments. According to this technique, we compute
$$F_Y (y) = P(Y\leq y) = P(X^2 \leq y) = P(|X|\leq \sqrt{y})$$
and since distribution functions define the density functions, after we get a simplified expression we just differentiate with respect to $y$ to get our pdf. We have then
$$ \begin{align} F_Y (y) = P\left( |X|\leq \sqrt{y} \right) = P\left(-\sqrt{y} < X <\sqrt{y} \right) = \Phi\left(\sqrt{y}\right) - \Phi \left(- \sqrt{y}\right) \end{align} $$
where $\Phi(.)$ denotes the CDF of a standard normal variable. Differentiating with respect to $y$ we get,
$$f_Y(y) = F_Y^{\prime} (y) = \frac{1}{2 \sqrt{y}} \phi \left( \sqrt{y} \right) + \frac{1}{2 \sqrt{y}} \phi \left( -\sqrt{y} \right) = \frac{1}{\sqrt{y}} \phi(\sqrt{y}) $$
where $\phi(.)$ is now the pdf of a standard normal variable and we have used the fact that it is symmetric about zero. Hence
$$f_Y (y) = \frac{1}{\sqrt{y}} \frac{1}{\sqrt{2\pi}} e^{-\frac{y}{2}}, \quad 0<y<\infty$$
which we recognize as the pdf of a chi-squared distribution with one degree of freedom (You might be seeing a pattern by now).
Density transformation technique
At this point you might wonder, why we do not simply use the transformation technique you are familiar with, that is, for a function $ Y = g(X)$ we have that the density of $Y$ is given by
$$f_Y (y) = \left| \frac{d}{dy} g^{-1} (y) \right| f_X \left( g^{-1} (y) \right)$$
for $y$ in the range of $g$. Unfortunately this theorem requires the transformation to be one-to-one which is clearly not the case here. Indeed, we can see that two values of $X$ result in the same value of $Y$, $g$ being a quadratic transformation. Therefore, this theorem is not applicable.
What is applicable, nevertheless, is an extension of it. Under this extension, we may decompose the support of $X$ (support means the points where the density is non-zero), into disjoint sets such that $Y= g(X)$ defines a one-to-one transformation from these sets into the range of $g$. The density of $Y$ is then given by the sum over all of these inverse functions and the corresponding absolute Jacobians. In the above notation
$$f_Y (y) = \sum \left| \frac{d}{dy} g^{-1} (y) \right| f_X \left( g^{-1} (y) \right)$$
where the sum runs over all inverse functions. This example will make it clear.
For $y = x^2$, we have two inverse functions, namely $x = \pm \sqrt{y}$ with corresponding absolute Jacobian $\displaystyle{ \frac{1}{2\sqrt{y}} }$ and so the corresponding pdf is found to be
$$f_Y (y) = \frac{1}{2\sqrt{y}} \frac{1}{\sqrt{2\pi} } e^{-y/2} + \frac{1}{2\sqrt{y}} \frac{1}{\sqrt{2\pi} } e^{-y/2} = \frac{1}{\sqrt{y}} \frac{1}{\sqrt{2\pi}} e^{-y/2}, \quad 0<y<\infty$$
the pdf of a chi-squared distribution with one degree of freedom. On a side note, I find this technique particularly useful as you no longer have to derive the CDF of the transformation. But of course, these are personal tastes.
So you can go to bed tonight completely assured that the square of a standard normal random variable follows the chi-squared distribution with one degree of freedom.