Distribution of sum of product-normal distributions. I know that the distribution of a product $Z=XY$ of two normally distributed variates $X$ and $Y$ with zero means is the product normal distribution [Mathworld].
What is the distribution of $Q=\sum Z$? By generating multiple samples of $X$ and $Y$, I can generate multiple variates $Z$ drawn from the product normal distribution.
I'm numerically experimenting with this, but would like to know if there can be an analytical expression.
[Mathworld]http://mathworld.wolfram.com/NormalProductDistribution.html
 A: Let $X \sim N(0,1)$ and $Y \sim N(0,1)$. Then, as noted, the pdf of $Z = X Y$ is $f(z)$:

The mgf of $Z$ is $E[e^{t Z}]$:

where I am using the Expect function from the mathStatica package for Mathematica to automate. Let $(Z_1, \dots, Z_n)$ denote a random sample of size $n$ drawn on $Z$, and let $Q = \sum_{i=1}^nZ_i$ denote the sample sum. Then, by independence, the characteristic function of $Q$ is:
$$\text{cf}(q) \quad = \quad \big(E\big[e^{i t Z}\big]\big)^n \quad = \quad \text{(mgf)}^{\ n} \big|_{t \rightarrow i t} \quad = \quad (1 + t^2)^{\ -\frac{n}{2}}$$
Inverting the characteristic function then yields a closed-form solution for the pdf of $Q$ as:
$$\text{pdf}(q) = \frac{\left| q\right| ^{\frac{n-1}{2}} K_{\frac{n-1}{2}}(\left| q\right| )}{2^{\frac{n-1}{2}} \sqrt{\pi} \,  \, \Gamma\left(\frac{n}{2}\right)}  $$
where $K_n(z)$ denotes the modified Bessel function of the second kind.
The following diagram plots the pdf of $Q$ when $n$ = 2, 3, 4 and 5.

Here is a Monte Carlo check when $n = 3$:

Notes:


*

*In the case of non-standardised variances, deriving a closed-form will be more difficult.

*As disclosure, I should perhaps add that I am one of the authors of the software used above.
