Please Help to solve the following

The earliest and the simplest known bivariate exponential distribution, introduced by Gumbel (1960), has joint survivor function and joint probability density function given by:
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 β.