# Is there a simple way to simulate $N$ binary random variables with marginal probability $p$ such that they all have correlation $\rho$?

I am wondering if there a simple way to simulate $N$ binary random variables $X_1, \ldots, X_N$ each with marginal probability $p$ such that they all have the condition of:

$$Corr(X_i, X_j) = \rho \qquad \forall \ \ i \neq j$$

? I found a case for $N=2$, but am wondering if this can be extended for all $N$. Any thought would be greatly appreciated.

If $\rho\in[0,1]$ the following method works: with probability $\alpha$ pick the $X_i$ iid, with probability $1-\alpha$ pick the $X_i$ all equal. If $\alpha=0$ the $X_i$ are perfectly correlated, if $\alpha=1$ they are perfectly uncorrelated and if $\alpha=1-\rho$ the $X_i$ have the desired correlation of $\rho$.
The case $\rho<0$ is trickier. In any case, since the covariance matrix of the $X_i$ must be positive semidefinite, we have $\rho\ge -1/(N-1)$, which for $N>2$ restricts the $\rho$ for any collection of rvs, binary or not. Even in the $N=2$ case the restriction to binary rvs makes things delicate. The pair $X_1=\mathbb 1[U\le p],$ $X_2=\mathbb1[1-U\le p],$ where $U$ is uniformly distributed on $[0,1]$ achieves the most negative $\rho$, which is $\rho_p= -(1-p)/p$ in the case $p\ge1/2$ and $\rho_p=-p/(1-p)$ if $p\le 1/2$. If the target $\rho$ does not lie in the interval $[\rho_p,0]$, no such binary random variables exist; if it does, stochastically mixing this recipe with the iid recipe will do the job.
I do not know how to generalize to the negative $\rho$ and $N>2$ case.