simulating dependent bernoulli variates

I have two binary random variables, say $y_A\sim Bernoulli\left(p_A\right)$ and $y_B\sim Bernoulli\left(p_B\right)$ that represent probability of default. I need to superimpose a sort of dependency on both and I was thinking to set $p_A=k_A*\tilde G$ and $p_B=k_B*\tilde G$ being $G$ a gamma variate with $E\left[G\right]=1$ and a certain variance. Nevertheless I am not sure if it is a sound approach since outcomes greater than one are theoretically possible.... Is there another approach to get the desired result?

• Macke et al. (2009) use dichotomized Gaussians to model dependent Bernoulli variables. Perhaps they also work for your use case? – Lucas Feb 22 '18 at 20:24
• One solution appears at stats.stackexchange.com/questions/279706. The idea is that after specifying $p_A$ and $p_B$, there's only one free parameter left, which you can use to determine the correlation. You can look at this as a categorical distribution on the four possible outcomes, which immediately gives a simple way to simulate from it. – whuber Feb 22 '18 at 20:55