Computes the probability density function of the 9-parameter JSBB distibution given by
$$
f_{Y_{1},Y_{2}}\bigl(y_1,y_2\big \vert\Theta\bigr) = f_{Y_1, Y_2}(y_1, y_2) =\frac{\delta_1\delta_2\lambda_1\lambda_2\exp\Bigl\{\frac{-z^{2}_{1}-z^{2}_{2} +2\rho z_{1}z_{2}}{2(1-\rho^2)}\Bigr\}}{2\pi \sqrt{1-\rho^2}\bigl(y_1-\xi_1\bigr)\bigl(y_2-\xi_2\bigr)\bigl(\lambda_1+\xi_1-y_1\bigr)\bigl(\lambda_2+\xi_2-y_2\bigr)},
$$
where
$$
z_{i}=\delta_i \log \Bigl(\frac{y_{i}-{\xi}_i}{{\xi}_i+{\lambda}_i-y_{i}}\Bigr)+\gamma_{i},
$$
for \(i=1,2\). The parameter space of SBB distribution is \(\Theta=({\bf{\delta}},{\bf{\gamma}},{\bf{\lambda}},{\bf{\xi}}, \rho)^{\top}\) in which \({\bf{\delta}}=(\delta_1,\delta_2)^{\top}\), \({\bf{\gamma}}=(\gamma_1,\gamma_2, \rho)^{\top}\), \({\bf{\lambda}}=(\lambda_1,\lambda_2)^{\top}\), and \({\bf{\xi}}=(\xi_1,\xi_2)^{\top}\). The supports of marginals are \(\xi_1<y_1<\lambda_1+\xi_1\) and \(\xi_2<y_2<\lambda_2+\xi_2\).
The support of the parameter space is \(\delta_1>0,\delta_2>0,-\infty<\gamma_1<+\infty,-\infty<\gamma_2<+\infty, \lambda_1>0,\lambda_2>0, -\infty<\xi_1<+\infty, -\infty<\xi_2<+\infty\) and \(-1<\rho<+1\).