Consider a pair of dependent zero mean unit variance Gaussians, $$X,Y \sim \mathcal{MVN}\left(\vec{0},\begin{pmatrix}1 & \rho \\ \rho & 1\end{pmatrix}\right).$$
Their product $Z:=X\cdot Y$ is known to follow the variance gamma distribuiton with density
$$p(z)=\frac{1}{\pi\sqrt{1-\rho^2}} \exp\left[\frac{\rho z}{1-\rho^2}\right] K\_{0}\left[\frac{\vert z\vert}{1-\rho^2}\right],$$ where $K\_0$ is the order zero modified Bessel function of the second kind (see [Nadaraja and Pogány, 2016](https://doi.org/10.1016/j.crma.2015.10.019)).


**Is the joint distribution of $X,Y,Z$ known?**


The distribution of $W:=aX +bY+cZ$ is also of interest.


**Note:** This is cross-posted from [math.stackexchange.com](https://math.stackexchange.com/questions/4051637/joint-distribution-of-dependent-gaussians-and-their-product). Given the recency of the results of which I'm looking for an extension, I think this might be more appropriate here. Let me know in the comments if I'm mistaken!



Joint distribution of dependent Gaussians and their product
The joint distribution of $X,Y,Z$ is of course known: it is the joint distribution of $X,Y,XY$, where $(X,Y)\sim N\_2(0,0,1,1,\rho)$.


However, the joint distribution of $X,Y,Z$ is not absolutely continuous (because $Z=XY$ is a function of $X$ and $Y$); that is, the joint distribution of $X,Y,Z$ does not a joint pdf.


However, it is easy to find the joint pdf of $Z$ and $X$ (say), and then just use the formula $Y=Z/X$ to compute any reasonable probabilities involving $X,Y,Z$. Indeed, using transformation technique (that is, change of variables in double integrals), we see that the joint pdf (say $g$) of $Z$ and $X$ is given by
$$g(z,x)=\frac1{|x|}\,f(x,z/x)$$
for real $z$ and $x\ne0$, where $f$ is the joint pdf of $(X,Y)\sim N\_2(0,0,1,1,\rho)$.


