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.

2$\begingroup$ No, the density of the sum of two independent random variables is the convolution of their density functions. Unfortunately, $X^2$ and $XZ$ are certainly not independent (although they are uncorrelated). $\endgroup$ – Robert Israel Jun 16 '11 at 19:07

$\begingroup$ @Robert Israel, thank you for your comment. Of course, you are absolutely right. I have edited my response to fix this error. $\endgroup$ – Tom LaGatta Jun 16 '11 at 20:20

2$\begingroup$ How does it solve the problem to represent XY as a sum of uncorrelated random variables with known distributions? $\endgroup$ – Douglas Zare Oct 1 '15 at 10:10
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 $\operatorname{E}XY= \operatorname{Cov}(X,Y) + (\operatorname{E}X)(\operatorname{E}Y)$.

$\begingroup$ That is a nice observation, pedro. If one is furthermore interested in any higher moments $\mathbb E (XY)^k$, then I suggest using the decomposition $Y = \tfrac{\rho \sigma_y}{\sigma_x}X + \sqrt{1\rho^2}\sigma_y Z$ and using the independence of $X$ and $Z$. As Arkadiusz says, there is likely no known closedform expression for the product, though in practice we can estimate probabilities using Chebyshev's inequality and arbitrary moments. $\endgroup$ – Tom LaGatta Jun 16 '11 at 20:26