[Math] What’s the densitiy of the product of two independent Gaussian random variables

Question

Suppose that $X,Y$ are two scalar independent normal random variables, $X \sim N(\mu_X,\sigma_X^2)$, $Y \sim N(\mu_Y,\sigma_Y^2)$. I'm particularly interested about the case where we don't assume $\mu_X = \mu_Y = 0$.

I'm interested in the random variable $XY$. What can be said about its PDF?

There's an existing question where an answer explains that $XY$ is the difference of two chi-squared variables.

For the zero-mean case, we know that the PDF is the normal product distribution. Is there a non-zero-mean generalization of this?

I know that there's a 1970 SIAM paper by Springer and Thompson, but I don't have access to this. Is the part which is relevant for my question publicly available somewhere?

To add to my confusion, I found a note by Bromiley, where it is argued that the product of two normal independent random variables is a normal variable again – which I thought was not the case. The argument in the linked document goes like this: Products of gaussian PDFs are gaussian. The PDF of a product of two independent RVs is their convolution. The Fourier transform of a convolution is the product of the fourier transforms. Gaussians are mapped to gaussians under the (inverse) Fourier transform. Am I misunderstanding something? Is there something wrong with the proof?

Answer 1

$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ \begin{align} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\rm P}\pars{\xi}&= \int_{-\infty}^{\infty}\dd x\,{1 \over \root{2\pi}\sigma}\, \exp\pars{-\,{x^{2} \over 2\sigma^{2}}} \\ & \ \int_{-\infty}^{\infty}\dd y\,{1 \over \root{2\pi}\sigma}\, \exp\pars{-\,{y^{2} \over 2\sigma^{2}}}\ \overbrace{\delta\pars{\xi - xy}}^{\ds{\delta\pars{y - \xi/x} \over \verts{x}}} \\[3mm]&={1 \over 2\pi\sigma^{2}}\int_{-\infty}^{\infty} \exp\pars{-\,{1 \over 2\sigma^{2}}\bracks{x^{2} + {\xi^{2} \over x^{2}}}}\, {\dd x \over \verts{x}} \\[3mm]&={1 \over \pi\sigma^{2}}\int_{0}^{\infty} \exp\pars{-\,{1 \over 2\sigma^{2}}\bracks{x^{2} + {\xi^{2} \over x^{2}}}}\, {\dd x \over x}\tag{1} \end{align} $\ds{\delta\pars{x}}$ is the Dirac Delta Function.

With $\ds{x \equiv A\expo{\theta/2}\,,\quad A > 0\,,\quad\theta \in {\mathbb R}}$: \begin{align} {\rm P}\pars{\xi}&={1 \over \pi\sigma^{2}} \int_{-\infty}^{\infty} \exp\pars{-\,{1 \over 2\sigma^{2}} \bracks{A^{2}\expo{\theta} + {\xi^{2} \over A^{2}}\expo{-\theta}}}\, \pars{A\expo{\theta/2}\,\dd\theta/2 \over A\expo{\theta/2}} \end{align} We can choose $\ds{A}$ such that $\ds{A^{2} = {\xi^{2} \over A^{2}}\quad\imp\quad A = \verts{\xi}^{1/2}}$: \begin{align} {\rm P}\pars{\xi}&={1 \over 2\pi\sigma^{2}} \int_{-\infty}^{\infty} \exp\pars{-\,{\verts{\xi} \over \sigma^{2}}\cosh\pars{\theta}}\,\dd\theta ={1 \over \pi\sigma^{2}} \int_{0}^{\infty} \exp\pars{-\,{\verts{\xi} \over \sigma^{2}}\cosh\pars{\theta}}\,\dd\theta \end{align}

$$\color{#00f}{\large% {\rm P}\pars{\xi} = {1 \over \pi\sigma^{2}}\, {\rm K}_{0}\pars{\verts{\xi} \over \phantom{2}\sigma^{2}}} $$ where $\ds{{\rm K}_{\nu}\pars{z}}$ is a Second Kind Bessel Function.