What is the expectation of exponential of the product of two random variables? I am looking for examples of probability distributions that would allow me to characterize the distribution (at least approximately) and to compute the first two moments exactly of:
$$
e^{aXY}
$$
where $a$ is a positive real and $X$ and $Y$ are scalar random variables defined on the positive real line and jointly distributed with mean/variance $(\mu_{X},\sigma^2_{X})$ and $(\mu_{Y},\sigma^2_{Y})$ and covariance $\sigma_{XY}$. 
Any suggestions very welcome. I'm looking for examples, so hopefully there should be more than one approach to this.
 A: I'm entering this as an "answer" to my own question even though it's probably not much of an answer. But hopefully this will help in showing what I'm after. I will edit it or delete it if asked. Edit: This is still work in progress. I should probably mention that Stat-R, cardinal, StasK, and Macro have guided me every step of the way to ... where I have stopped now. Please do feel free to edit/correct this draft.
I'm looking for examples of a joint distributions of the dependent rvs $X$ and $Y$ such that (1) the distribution can be characterized, at least approximately, (2) the following can be computed exactly:
$$ 
E[e^{aXY}] 
$$
Example 1. Let:
$$ 
X = x + \epsilon 
$$
$$ 
Y = y + \epsilon 
$$
$$ 
\epsilon \sim \mathcal{N}(0,\sigma^{2}) 
$$
where $x$ and $y$ are fixed real numbers (the known means of $X$ and $Y$). 
Remark: 
$$ 
e^{\epsilon} \sim \mathcal{LN}(0,\sigma^{2})
$$
$$
E[e^{\epsilon}]=e^{\sigma^{2}/2}
$$
It follows from a substitution of $x+\epsilon$ and $y+\epsilon$ into $XY$ that:
$$ 
E[e^{aXY}] = e^{axy} \times E[e^{a(x+y)\epsilon} e^{a\epsilon^{2}}] 
$$
Edit: The $\epsilon^{2}$ is not negligeable, as explained in the comments. Thus,
$$ E[e^{aXY}] = e^{axy} \times E[e^{a(x+y)\epsilon}] \times E[e^{a\epsilon^{2}}] + e^{2axy} \times COV[e^{a(x+y)\epsilon},e^{a\epsilon^{2}}]
$$
The first expectation on the rhs:
$$
E[e^{a(x+y)\epsilon}] = e^{a^{2}(x+y)^{2}\sigma^{2}/2} 
$$
The second expectation on the rhs features the square of a Normal, which is a Chi-squared. Edit: I have been shown, in the comments, how to compute the expectation by exploiting the fact that it's an evaluation of the MGF of a chi-squared, since $(\epsilon/\sigma)^{2}\sim\chi_{1}^{2}$. Therefore, 
$$
E[e^{a\epsilon^{2}}] = (1-2a/\sigma^{2})^{-k/2} = 1/\sqrt{1-2a/\sigma^{2}}, \ a < \sigma^{2}/2
$$
where $k$ denotes the degrees of freedom, which in this case is $k=1$, as the chi-squared is obtained from $1$ squared Normal. If $a>\sigma^{2}/2$, then the expectation diverges.
Lastly, the covariance term. I think here I need to compute this from the integral, hopefully some known results can be used along the way. If someone knows how to finish this it would be great, but I'm stuck.
