Is it possible to simulate pairs of random variables with a given marginal distribution and population correlation where one random variable is larger than the other?
More formally, I need to simulate pairs of random variables $(X_1, Y_1), \dots(X_n, Y_n),$ where $X_1, \dots X_n \sim f(\cdot;\Theta)$, $Y_1, \dots Y_n \sim g(; \Psi)$, $f$ and $g$ are continuous probability distributions, and $Cor(X_i, Y_i) = \rho$. Is it possible to simulate these pairs of random variables such that $X_i \leq Y_i$ for all $1 \leq i \leq n$?
I can weaken the conditions slightly: a strict inequality, $X_i < Y_i$, would be fine, the population or sample correlation can close to $\rho$, $X$ and $Y$ are from the same family of distributions but with different parameters.
Without the inequality constraint, this is easy to simulate using a copula: Generate $(W_{i1} ,W_{i2}) \sim \mathcal{N}_2(0, \Sigma)$, use probability-integral-transform to get correlated $U[0,1]$ random variables, then plug them into the marginals. I don't know how to do this with the constraint. My instinct is to keep working with a copula of some kind, but I am not tied to this approach. Is this kind of random variable generation possible?