Here's a pretty general construction. Take unit vectors $v_1,\dots,v_n$ in $\mathbb{R}^n$ and let $u$ be a random unit vector (chosen with the uniform probability measure on the unit sphere). Define $B_i$ to be 1 if the inner product of $u$ and $v_i$ is positive and -1 otherwise. Then the correlation between $B_i$ and $B_j$ is the inner product of $v_i$ and $v_j$. (I haven't checked that carefully but I think it's true.)
Added in the light of the comment below: OK, I should have checked. It's actually not the inner product but π minus twice the angle between the two vectors all over π. I.e., it depends linearly on the angle between the two vectors, is 1 when that angle is zero and -1 when it is π. The angle is the inverse cos of the inner product, which gives us a formula.
So it gives us a fairly big supply of matrices -- I can't quite decide whether it gives us all (or rather all for which the variables take two values, each with probability 1/2).
A more general-looking construction is this. Take any probability space and let $A_1,...,A_n$ be sets of measure 1/2. Pick a random point x and let $B_i$ be 1 if x is in $A_i$ and -1 otherwise. But that becomes trivial, because if you have any set of Bernoulli variables taking the values $\pm 1$ with probability 1/2, then you can set $A_i$ to be the set where $B_i=1$.
