[notafraudster]

You can arbitrarily (rank-)correlate any variables of any distribution using copulas as an intermediary.

So basically, draw a multivariate correlated standard normal with the intended correlation. Transform the standard normal draws into quantiles of the standard normal (this works the marginal distributions are also normal). Now you have X and Y quantiles for your target distributions and you can draw from them.

The distributional transformations will slightly attenuate the correlation, and the choppier the distributions the more attenuation you'll get. Additionally, while you can do this for 3+ variables at once, there are constraints on the possible sigma matrix describing n-dimensional correlations.

If the variables have no fixed distribution you can use the eCDF of the actual data, so it's even possible to import e.g. population and income data and produce a permutation that gives you the correlation desired.

I agree that it is fairly difficult to do this if you have arbitrary DGPs with complex interdependencies in them--if you correlate on variables A and B, then it is difficult to guarantee the observed correlation between f(A) and g(B) and vice versa--but still you can provide a lot of utility with the copula method.