Let $X \sim \text{N}(0, \sigma_X^2)$ and $Y \sim \text{N}(0, \sigma_Y^2)$ be independent normal random variables with zero mean, but (possibly) different variances. Given some constants $a$ and $b$, I would like to obtain an expression for:
$$\mathbb{E}(\exp( ab X Y)).$$
I am able to solve the univariate case where I only have one of these random variables in the exponential, but I cannot solve this bivariate case.
Without loss of generality, the factor $ab$ can be omitted as including this factor has the same effect as changing $\sigma_X^2\sigma_Y^2$ by a factor of $a^2b^2$.
Using the law of total expectation, the definition of the moment generating function, the mgfs of the normal and chi-square distributions, \begin{align} Ee^{XY} &=E(Ee^{XY}|Y) \\&=EM_X(Y) \\&=Ee^{\frac12\sigma_X^2Y^2} \\&=Ee^{\frac12\sigma_X^2\sigma_Y^2\chi_1^2} \\&=M_{\chi_1^2}(\frac12\sigma_X^2\sigma_Y^2) \\&=(1 - \sigma_X^2\sigma_Y^2)^{-1/2}. \end{align} This only works for $\sigma_X\sigma_Y<1$ which suggest that the expectation is perhaps not finite when $\sigma_X\sigma_Y\ge 1$.
Verifying the result against simulations:
> X <- rnorm(1e+6,sd=.9)
> Y <- rnorm(1e+6,sd=.7)
> mean(exp(X*Y))
[1] 1.284138
> (1-.9^2*.7^2)^(-.5)
[1] 1.287672
