Distribution of $XY$ when $(X,Y) \sim BVN(0,0,1,1,\rho)$ The question is pretty much in the title,
I need to find an approximate distribution of $XY$ when $(X,Y)$ follow a Bivariate Normal Distribution where $X$ and $Y$ are each $N(0,1)$ distributed and $Cov(X,Y) = \rho.$
I proceeded to find the required using transformation of variables, but I am stuck at the following integration,
$\displaystyle \int_{-\infty}^{\infty} \left| \frac{1}{u}\right| \text{exp} \left\{ \frac{-1}{2(1-\rho^2)}\left(u^2 + \frac{v^2}{u^2}\right)\right\}\,du$
How do I proceed to solve this integral? Or is there any other way obtaining the required distribution?
This problem is however an intermediate step of another question that I am trying to solve,

Suppose ${(X_1, Y_1), . . . ,(X_n, Y_n)}$ is a random sample from a bivariate normal distribution with $E(X_i) = E(Y_i) = 0, \, Var(X_i) = Var(Y_i) = 1$ and unknown $Corr(X_i, Y_i) = \rho \in (−1, 1),$ for all
  $i = 1, . . . , n.$ Define $W_n = \frac{1}{n}\sum_{i = 1}^n X_iY_i$.
  For large $n$, obtain an approximate level $(1 − \alpha)$ two-sided
  confidence interval for $\rho$, where $0 < \alpha < 1.$

Though, I know about the Fisher's z transform, and its asymptotic distribution. But I think the question requires me to give an answer in terms of $W_n$ only. How should I solve this? 
Can I use CLT in this? Because I ran a simulation for 1000 times and the distribution of $W_n$ came out to be convincingly normal. 
 A: Your original question does not require you to find the distribution of $XY$ when $(X,Y)$ is jointly normal. Here is a hint for that question:
Let $Z_i=X_iY_i$, so that $Z_1,Z_2,\ldots,Z_n$ are i.i.d variables. Hence by classical CLT we have
$$\frac{\sqrt n(W_n-\operatorname E(Z_1))}{\sqrt{\operatorname{Var}(Z_1)}}\stackrel{L}\longrightarrow N(0,1)\tag{*}$$
, where $W_n=\frac{1}{n} \sum\limits_{i=1}^n Z_i$.
Since we know the conditional distributions of a bivariate normal distribution, use law of total expectation to find $\operatorname E(Z_1)$:
$$\operatorname E(Z_1)=\operatorname E\left[\operatorname E(X_1Y_1\mid Y_1)\right]=\operatorname E\left[Y_1 \operatorname E(X_1\mid Y_1)\right]$$
Do this similarly for $\operatorname E(Z_1^2)$ and hence find $\operatorname{Var}(Z_1)=\operatorname{E}(Z_1^2)-[\operatorname{E}(Z_1)]^2$
Now $(*)$ actually gives you a pivot for constructing an asymptotic confidence interval for $\rho$. 
A: An exact answer is given in Probability Distributions Involving Gaussian Random Variables: A Handbook for Engineers, Scientists and Mathematicians (along with many other results.)
If $(X,Y)$ is bivariate normal distributed with zero means and correlation $\rho$, then the density function of the product $XY$ is given by
$$
   f_{XY}(z)= \frac1{\pi \sigma_1 \sigma_2} \exp\left\{\frac{\rho z}{\sigma_1 \sigma_2 (1-\rho^2)}\right\} K_0\left( \frac{|z|}{\sigma_1\sigma_2 (1-\rho^2)} \right)
$$ where $K_0$ is the modified Bessel function of the second kind. 
EDIT   answer for additional question in comment

The reference book given gives the characteristic function. For the moment generating function I get $M_{XY}(t)= \left(1 - 2 t \rho - t^2(1-\rho^2)\right)^{-\frac{1}{2}}$ for $-1<t-1$ (so you missed on $t$.) 
