Square root of normal distribution

Question

Let $X$ and $Y$ be independent random variates with the same probability distribution, $P(x)$. Assuming that the product $Z=XY$ is a random variate with normal distribution, say $$f_Z(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2} x^2}$$ what is the form of $P(x)$? Does it exist in closed form?

EDIT: I just realized that a product distribution is necessarily divergent at zero, which means that the normal distribution cannot be a product distribution (at least when $P(x)$ is interpreted as a normal function; maybe there is a generalized function that solves the problem?). Do you agree?

Answer 1

First: No helvio, the density of the product of two iid random variables (r.v.'s) is not necessarily divergent at 0. E.g., take the product of two iid r.v.'s each having a Gamma distribution with shape parameter 2.

Second: Let $U:=\ln|Z|$ and $W:=\ln|X|$. Then the characteristic function (c.f.) $f$ of $U$ is given by the formula $f(t)=2^{i t/2} \Gamma((1+i t)/2)/\sqrt{\pi}$ for real $t$. So, the c.f. of $W$ (say $g$) would be a square root of $f$, and then the density (say $q$) of $W$ would be given by the formula $q(x)=\frac1{2\pi}\,\int_{-\infty }^{\infty } g(t) e^{-i t x} \, dt$ for real $x$. I have tried to evaluate the latter integral numerically, and that seems to suggest that $q(3)<0$; however, I have not done a rigorous estimate of the approximation error.

On the other hand, numerical experiments using Bochner's theorem on the positive definite functions suggest that $g$ is a c.f.! If so, then the corresponding r.v. $W$ exists, and then it suffices to let $X$ have the distribution of $\varepsilon e^W$, where $\varepsilon$ is a Rademacher r.v. independent of $W$.

Addendum: The calculation of $q(3)$ mentioned previously was likely mistaken, as I forgot to make sure that $g(t)$ be taken to be the value of the square root of $z:=f(t)$ on the Riemann surface for $\sqrt z$, so as for $g$ to be a continuous function -- as it must be.

Another approach to the problem: Note that $\mu_j:=E|Z|^j=2^{j/2} \Gamma \left(\frac{j+1}{2}\right)/\sqrt{\pi }$ and then one would have $\nu_j:=E|X|^j=\sqrt{\mu_j}$ for $j>0$. To prove that such a desired r.v. $X$ exists, it is enough to show that for all nonnegative integers $n$ the determinants, say $d_n$ and $e_n$, of the matrices $(\nu_{i+j})_{i,j=0}^n$ and $(\nu_{i+j+1})_{i,j=0}^n$ are nonnegative (see e.g. Ch. V in Kreǐn, M. G. and Nudel'man, A. A., The Markov moment problem and extremal problems). These determinants are indeed positive for $n\le9$. One can hopefully find a recursion for $d_n$ and $e_n$ to show that they are indeed nonnegative for all $n$. Note here that $\mu_j$ is an integer or an integral multiple of $\sqrt{2/\pi}$.