Distribution of a product of normal distributions : why am I wrong?

Question

Let $X$ and $Y$ be two independent normal distributions with mean $0$ and variance $1$ for simplicity.

I want to find the distribution of $XY$.

Attempt :

$P(XY=w)=\int_{-\infty}^{+\infty}P(X=s)P(Y=\frac{w}{s})ds=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-\frac12s^2-\frac12w^2/s^2}ds$

Using $\int_{-\infty}^{+\infty}e^{-ax^2-b/x^2}dx=\sqrt{\frac{\pi}{a}}e^{-2\sqrt{ab}}$ (Glasser's theorem or variable change)

with $a=\frac12$ and $b=\frac 12w^2$

I obtain :

$P(XY=w)=\frac{1}{\sqrt{2\pi}}e^{-|w|}$

but integrating over $w$ yields $\sqrt{\frac2\pi}\approx0.8<1$

So, which cases am I missing ? Is it because $s=0$ is included and it is illicit ? Even then, I would expect to get something strictly superior to one and not inferior since the integrand is always positive... unless I am missing special cases ?

I know I can find the distribution of $XY$ by a google search, but I'd still like to know where I am making a mistake, so that it doesn't happen again. Sorry if this is trivial, sometimes I just can't see it.

Answer 1

Firstly, your notation is confused.   You should not use a probability mass function when you mean a probability density function.

Secondly, the convolution is based on the chain rule of differentiation and the law of total probability.

$$\begin{align}f_{XY}(w) ~&=~ \int_\Bbb R \underset{\text{Jacobian Determinant}}{\underbrace{\left\lVert\frac{\partial(s,w/s)}{\partial (s,w)}\right\rVert} f_X(s)}f_Y(w/s)\operatorname d s \\[1ex] &=~ \frac 1{2\pi} \int_\Bbb R\lvert s^{-1}\rvert \exp(-s^2/2)\exp(-w^2/2s^2)\operatorname d s \\[1ex] &=~\frac 1\pi \int_0^\infty s^{-1}\,\mathsf e^{-(s^2+w^2/s^2)/2}\operatorname d s\end{align}$$

Thirdly, that's not going to resolve into elementary functions.

(Hint Topic: Modified Bessel Function of the Second Kind.)

Answer 2

I've redone the calculations an elementary way, which, in my opinion is better for beginners. (It is lengthier though.)

The pdf of $XY$ can be interpreted the following way

$$f_{XY}(w)=\lim_{\Delta w\to 0}\frac{F_{XY}(w+\Delta w)-F_{XY}(w)}{\Delta w}=\lim_{\Delta w\to 0}\frac{P(w\le XY<w+\Delta w)}{\Delta w}.$$

Let's deal with the probability in the numerator

$$P(w\le XY<w+\Delta w)=\frac1{\sqrt{2\pi}}\int_{-\infty}^{\infty}P(w\le XY<w+\Delta w\mid X=x)e^{-\frac{x^2}2} \ dx=$$ $$=\frac1{\sqrt{2\pi}}\int_{-\infty}^{\infty}P\left(w\le xY<w+\Delta w \right)e^{-\frac{x^2}2} \ dx.$$

So,

$$\lim_{\Delta w\to 0}\frac{P(w\le xY<w+\Delta w)}{\Delta w}=$$ $$=\frac1{\sqrt{2\pi}}\int_{-\infty}^{\infty}\lim_{\Delta w\to 0}\frac{P\left(w\le xY<w+\Delta w\right)}{\Delta w}e^{-\frac{x^2}2} \ dx.\tag 1$$

Now, if $x>0$

$$\lim_{\Delta w\to 0}\frac{P\left(\frac wx\le Y<\frac wx+\frac {\Delta w}x\right)}{\Delta w}=\lim_{\Delta w\to 0}\frac{F_Y\left(\frac wx+\frac{\Delta w}x\right)-F_Y\left(\frac wx\right)}{\Delta w}=$$ $$=\frac{d F_Y\left(\frac wx\right)}{dw}=\frac1xf_Y\left(\frac wx\right)=\frac1{\mid x \mid}\frac1{\sqrt{2\pi}}e^{-(\frac wx)^2\frac12}$$ If $x<0$ $$\lim_{\Delta w\to 0}\frac{P\left(\frac wx\ge Y>\frac wx+\frac {\Delta w}x\right)}{\Delta w}=\lim_{\Delta w\to 0}\frac{F_Y\left(\frac wx\right)-F_Y\left(\frac wx+\frac {\Delta w}x\right)}{\Delta w}=$$ $$=-\lim_{\Delta w\to 0}\frac{F_Y\left(\frac wx+\frac {\Delta w}x\right)-F_Y\left(\frac wx\right)}{\Delta w}=$$ $$=\frac1{\mid x \mid }f_Y\left(\frac wx\right)=\frac1{\mid x\mid}\frac1{\sqrt{2\pi}}e^{-(\frac wx)^2\frac12}$$

The case $x=0$ can be considered as if it was omitted from the domain of the integral $(1)$ which looks like this now

$$f_{XY}(w)=\frac1{2\pi}\int_{-\infty}^{\infty}\frac{e^{-\frac12\left((\frac wx)^2+{x^2}\right)}}{\mid x \mid } \ dx=\frac1{\pi}\int_0^{\infty}\ \frac{e^{-\frac12\left((\frac wx)^2+{x^2}\right)}}{x}dx.$$