I am trying to figure out the following problem. I am trying to evaluate the expectation: $E[\exp(-sX)\exp(-sY)]$, where $X$ and $Y$ are identically distributed, but correlated random variables, hence, they are NOT independent. I have access to the Laplace transform of the distribution (i.e., $E[\exp(-sX)]$) but not the exact PDF. I know that it is possible to take the inverse Laplace transform to get the PDF, but it might not be necessary due to the structure of the problem. I can also calculate the correlation coefficient $\rho(X,Y)$ from numerical data. So far I have done the following: $E[\exp(-sX)\exp(-sY)]=\\E_X[E_{Y|X}[\exp(-sX)\exp(-sY)]]=\\E_X[\exp(-sX)E_{Y|X}[\exp(-sY)]].$
Is there a simplification I am not seeing, or a way to incorporate the correlation coefficient and the known Laplace transforms $E[\exp(-sX)]$ and $E[\exp(-sY)]$.
Edit: What about the following approach: $Cov(A,B)=E[AB]-E[A]E[B]$, where $A=\exp(-sX)$, and $B=\exp(-sY)$. Therefore, $E[A]=\mathcal{L}_X(s)$ and $E[B]=\mathcal{L}_Y(s)$. Then, we calculate $Cov(A,B)$ from numerical data (simulations) and we have: $E[\exp(-sX)\exp(-sY)]=\mathcal{L}_X(s)\mathcal{L}_Y(s)+Cov(A,B)$
Suppose $(X,Y)$ is a 2-dimensional multivariate Gaussian with respective means $\mu_X,\mu_Y$, variances $\sigma^2_X,\sigma^2_Y$ and covariance $\rho_{XY}$. Note that I'm using the covariance rather than correlation just because it is simpler. Then $(X,Y) \sim \mathcal{N}(\mathbf{\mu},\mathbf{\Sigma})$, where $$\mathbf{\mu} =\begin{bmatrix} \mu_X\\ \mu_Y \end{bmatrix} \text{ and } \mathbf{\Sigma} = \begin{bmatrix} \sigma^2_X & \rho_{XY}\\ \rho_{XY} & \sigma^2_Y \end{bmatrix}.$$ Then the moment generating function $M_{XY}:\mathbb{R}^2\to\mathbb{R}$ of $(X,Y)$ is known to be given by: $$M_{XY}(\mathbf{t}) = \mathbb{E}\left[e^{\mathbf{t}\cdot (X,Y)^T}\right] = \exp\left(\mathbf{\mu}\cdot\mathbf{t} + \frac{1}{2}\mathbf{t}^T\mathbf{\Sigma}\mathbf{t}\right).$$ So, $$\mathbb{E}\left[\exp(-sX)\exp(-sY)\right] = M_{XY}(-s,-s) = \exp\left(-s\mu_X + \frac{s^2}{2}\sigma^2_X\right)\exp\left(-s\mu_Y + \frac{s^2}{2}\sigma^2_Y\right)\exp\left(s^2\rho_{XY}\right) = \mathbb{E}\left[\exp(-sX)\right]\mathbb{E}\left[\exp(-sY)\right]\exp(s^2\rho_{XY}).$$
Now let's assume that $X$ and $Y$ are only approximately Gaussian. Then when is this approximation appropriate? Let me first present a standard result for MGFs:
If the MGFs $M_{X'_n,Y'_n}(t)$ $(X'_n,Y'_n)$ are uniformly bounded for $|t| < |(-s,-s)|+\epsilon$ for some $\epsilon > 0$ and if $(X'_n,Y'_n)$ converge weakly to $(X',Y')$ then $$\mathbb{E}\left[\exp(-sX'_n-sY'_n)\right] = M_{X'_n,Y'_n}(-s) \to M_{X'Y'}(-s) = \mathbb{E}\left[\exp(-sX'-sY')\right].$$
How can we apply this? You said that you should be able to compute $\mathbb{E}[\exp(-sX)]$ for different values of $s$ and $X$ and $Y$ are identically distributed. Then by Holder's inequality, $$\mathbb{E}\left[\exp(-sX)\exp(-sY)\right]\leq \sqrt{\mathbb{E}[\exp(-2sX)]\mathbb{E}[\exp(-2sY)]} = \mathbb{E}[\exp(-2sX)].$$ So if $\mathbb{E}[\exp(-2s'X)]$ is not large compared to the Gaussian approximation of the MGF derived above when $s'$ is as close or closer to the origin than $(-s,-s)$, you can assume your MGF is well-behaved. Then the approximation will be about as good as your pair $(X,Y)$ is close to being Gaussian.
Hope that helps!