Uncertainty/Standard Deviation of Monte Carlo methods I am using a Monte Carlo method to estimate the expected value of the results of certain simulations.
Consider this simplified case: $X, Y$ are independent random variables and $g(X,Y)$ is a nonlinear function of $X,Y$. I would like to use Monte Carlo to estimate $\mathbb{E}_{X,Y} g(X,Y)$.
In my code I do the following:

*

*Draw $N_Y$ random samples ($y_1, \dots, y_{N_y}$) of $Y$


*For each $y_j$, draw $N_X$ fresh new random samples ($x_1^{(j)}, \dots, x_{N_x}^{(j)}$) of $X$ and compute
$$ g_j = \frac{1}{N_x} \sum_{i=1}^{N_x} g(x_i^{(j)}, y_j)$$


*I then estimate
$$\mathbb{E}g(X,Y) \approx  \frac{1}{N_y} \sum_{j=1}^{N_y} g_j$$
The reason why I am using this method (rather than just take $N_x \times N_y$ random samples $(x_i, y_i)$) is because it is computationally heavy to draw random samples of $Y$ (large Erdos-Renyi graphs), and in this way, I can take $N_y << N_x$.
I know that because of the linearity of the expected value, the above estimate is indeed $\mathbb{E}g(X,Y)$. My question is regarding its uncertainties.
My question is:

*

*How do you compute the uncertainties for it, e.g. std?  How do they compare with the standard Monte Carlo method?

Edit after @Xi'an's answer
For easiness let me introduce some notation. There are three possible estimators:

*

*Draw $N_y$ random samples  of $Y$ and $N_x$ random samples of $X$, then estimate  $\mathbb{E}g(X,Y)$ by $S_1$
$$ S_1 =  \frac{1}{N_y N_x} \sum_{j=1}^{N_y} \sum_{i=1}^{N_x}  g(x_i, y_j)$$


*Draw $N_y \times N_x$ random pairs of $Y$ and $X$, then estimate  $\mathbb{E}g(X,Y)$ by $S_2$
$$ S_2 =  \frac{1}{N_y N_x} \sum_{j=1}^{N_y N_x} g(x_j, y_j)$$


*Draw $N_y$ random samples of $Y$ and for each $y_j$, draw $N_X$ fresh new random samples ($x_1^{(j)}, \dots, x_{N_x}^{(j)}$) of $X$ , then estimate  $\mathbb{E}g(X,Y)$ by $S_3$
$$ S_3 =  \frac{1}{N_y N_x} \sum_{j=1}^{N_y} \sum_{i=1}^{N_x}  g(x_i^{(j)}, y_j)$$
My question is: Can we say anything, a priori, about the variance of these estimators or when one outperforms the others?
Notice that these estimators have very different amount of "randomness".

Since $X$ and $Y$ are independent, there is not need to simulate new samples $x^{(j)}$ for different $j$'s. That is,
$$\frac{1}{N_xN_y}\sum_{i=1}^{N_x}\sum_{j=1}^{N_y}g(x_i,y_j)\tag{1}$$
is an unbiased estimator of $\mathbb E_{X,Y}[g(X,Y)]$. Indeed, for all pairs $(i,j)$,
$$\mathbb E_{X,Y}[g(X_i,Y_j)]=\int_\mathfrak{X}\int_\mathfrak{Y}g(x,y)f_X(x)f_Y(y)\,\text dx\text dy$$
The variance of (1) may be larger than the variance of
$$\frac{1}{N_xN_y}\sum_{i=1}^{N_x}\sum_{j=1}^{N_y}g(x_i^{(j)},y_j)\tag{2}$$
since $g(X_i,Y_j)$ is likely to be positively correlated with $g(X_i,Y_k)$ for $j\ne k$, but the cost of producing the $N_x\times N_y$ $X_i^{(j)}$'s may offset the gain in variance provided by (2).
