S(x,y)=exp[−(αx+βy+θαβxy)] and f(x,y)=[(1−θ)αβ+θα2βx+θαβ2y+θ2α2β2xy]S(x,y) respectively, where x,y>0, α,β≥0 and 0<θ<1.

I want to generate random samples from this distribution using mathematica. It is important to point out that the marginal distributions of X and Y are exponential with parameters α and β.