Suppose we have pair $(X,Y)\sim Normal([\mu_x,\mu_y],{{\sigma_x^2\atop\rho \sigma_x\sigma_y } {\rho \sigma_x\sigma_y \atop \sigma_y^2} }] $ How is $U=X\cdot Y$ distributed? I've tried to compute this by substituting y=u/x in the bivariate normal pdf and taking integral(from $\infty$ to $\infty$ with respect to x. I find the pdf of U as the sum of two exponential distributions(one for U<0 and one for U>0) that are weighted unequally. Is my method valid, or do I have to deal with cumulative distribution functions instead?
Arkadiusz gives the answer in the case of two independent Gaussians. A simple technique to reduce the correlated case to the uncorrelated is to diagonalize the system. The intuition which I use is that for two random variables, we need two "independent streams of randomness," which we then mix to get the right correlation structure. Let $X \sim N(0,\sigma_x)$ and let $Z \sim N(0,1)$ be two independent normals. Define $Y = \tfrac{\rho \sigma_y}{\sigma_x}X + \sqrt{1\rho^2}\sigma_y Z$. Check that $\mathbb E Y^2 = \sigma_y^2$ and $\mathbb E XY = \rho \sigma_x \sigma_y$; this completely determines the bivariate Gaussian case you're interested in. Now, $XY = \tfrac{\rho \sigma_y}{\sigma_x} X^2 + \sqrt{1\rho^2}\sigma_y XZ$. The $X^2$ part has a $\chi^2$distribution, familiar to statistics students; the $XZ$ part is comprised of two independent Gaussians, hence Arkadiusz's answer gives the distribution of that random variable. Edit: As Robert Israel points out in the comments, I made a mistake in my final conclusion: the random variables $X^2$ and $XZ$ are uncorrelated, though certainly not independent. Nonetheless, the problem is essentially resolved at this point, since we have reduced the problem of understanding the product $XY$ to a sum of uncorrelated random variables $X^2$ and $XZ$ with known distributions. 


It seems that no closedform expression is known for the correlated case. For the uncorrelated case $\rho=0$ the distribution of XY is $\frac{1}{\pi \sigma_x \sigma_y}K_0(\frac{u}{\sigma_x \sigma_y})$, where $K_0(x)$ is the modified Bessel function of the second kind. This distribution differs from the distribution you gave (is it Laplace distribution?). A notable difference: it has kurtosis of 6 (sharper peak), comparing to 3 for Laplace distribution. On the Frequency Function of xy. C.Craig, 1936 Edit: answering the second question, the method is valid, but my guess is that the Jacobian determinant was skipped while doing the substitution. 


LaGatta's answer nails it, and may be useful for drawing simulations, etc. This is just a note to remind that if one is only interested in the mean of the product of normallydistributed (possibly correlated) random variables, then the answer is straightforward, using the identity E[x.y] = Cov(x,y) + E[x].E[y]. 

