How to generate random samples from Gumbel’s bivariate exponential distribution?

Question

The earliest and the simplest known bivariate exponential distribution, introduced by Gumbel (1960), has joint survivor function and joint probability density function given by:

\begin{equation}\label{eq:sf} S(x, y) = \exp\left[- \left(\alpha\,x + \beta\,y + \theta\,\alpha\,\beta\,x\,y\right)\right] \end{equation} and \begin{equation}\label{eq:pdf} f(x, y) = \left[(1 - \theta)\,\alpha\,\beta + \theta\,\alpha^2\,\beta\,x + \theta\,\alpha\,\beta^2\,y + \theta^2\,\alpha^2\,\beta^2\,x\,y\right]\,S(x, y) \end{equation} respectively, where $ x, y >0 $, $ \alpha, \beta \geq 0 $ and $ 0 < \theta < 1 $.

I want to generate random samples from this distribution. It is important to point out that the marginal distributions of X and Y are exponential with parameters $\alpha$ and $\beta$.

Answer 1

Given the joint density\begin{equation} f(x, y) = \left[(1 - \theta)\,\alpha\,\beta + \theta\,\alpha^2\,\beta\,x + \theta\,\alpha\,\beta^2\,y + \theta^2\,\alpha^2\,\beta^2\,x\,y\right]\,S(x, y) \end{equation} the conditional density of $Y$ given $X=x$ is such that \begin{align} f(y|x) &\propto \left[(1 - \theta)\,\alpha\,\beta + \theta\,\alpha^2\,\beta\,x\right]\, S(x, y)+ \left[\theta\,\alpha\,\beta^2 + \theta^2\,\alpha^2\,\beta^2\,x\right]\,y\,S(x, y)\\ &\propto \left[1 - \theta+ \theta\,\alpha x\right]\, S(x, y)+ \left[\theta + \theta^2\,\alpha\,\beta\,x\right]\,y\,S(x, y)\\ &\propto \left[1 - \theta+ \theta\,\alpha x\right]\, \exp\left[- \left(\beta\,y + \theta\,\alpha\,\beta\,x\,y\right)\right]+ \left[\theta + \theta^2\,\alpha\,\beta\,x\right]\,y\,\exp\left[- \left( \beta\,y + \theta\,\alpha\,\beta\,x\,y\right)\right]\\ &\propto \left[1 - \theta+ \theta\,\alpha x\right]\, \exp\left[- \beta\,y\left( 1+ \theta\,\alpha\,x\right)\right]+ \left[\theta + \theta^2\,\alpha\,\beta\,x\right]\,y\,\exp\left[- \beta\,y\left( 1 + \theta\,\alpha x \right)\right]\\\end{align} which means that the conditional distribution of $Y$ given $X=x$ is a mixture of two distributions,

1. the exponential $\mathcal{E}(\beta\,\left\{ 1+ \theta\,\alpha\,x\right\})$ distribution
2. the Gamma $\mathcal{G}a(2,\beta\,\left\{ 1+ \theta\,\alpha\,x\right\})$ distribution

with respective weights $1$ and $$ \beta^{-1}\left[\theta + \theta^2\,\alpha\,\beta\,x\right]\times\left( 1 + \theta\,\alpha x \right)^{-2}$$up to normalisation. It is therefore straightforward to simulate $X$ marginally and then $Y$ given this realisation of $X$.