Let $X$ and $Y$ be two independent continuous random variables with pdfs $f_X$ and $f_Y$, respectively. Let $\varphi_1$ and $\varphi_2$ be two continuous functions from ${\mathbb R}$ to ${\mathbb R}$. I want to calculate $E[\varphi_1(X)\varphi_2(Y)]< \infty$ numerically. This is,
$$I = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \varphi_1(x)\varphi_2(y) f_X(x)f_Y(y)dxdy\\ =\int_{-\infty}^{\infty} \varphi_1(x)f_X(x)dx\int_{-\infty}^{\infty}\varphi_2(y)f_Y(y) dy< \infty.$$
Using Monte Carlo integration, I can either approximate $I$ using $$I \approx \frac{1}{n}\sum_{j=1}^{n} \varphi_1(x_j)\varphi_2(y_j),$$ where $(x_j,y_j)$ is an independent sample from the joint distribution of $(X,Y)$. Alternatively, I can use the approximation: $$I \approx \left[\frac{1}{n}\sum_{j=1}^{n} \varphi_1(x_j)\right]\left[\frac{1}{n}\sum_{j=1}^{n} \varphi_2(y_j)\right] .$$ In my case, and minding potential implementation errors, I am obtaining different results.
Question. Is there any reasons to prefer one approximation over the other?
Simulating from the pdfs is straightforward in my case.