Let $\theta_i$ be a Bernoulli random variable of success probability $p$.
Given two fixed sets $A = \{a_1, a_2, ..., a_N\}$ and $B = \{b_1, b_2, ..., b_N\}$ where $a_i$ and $b_i$ are constants.
For any realisations $Z = \{\theta_1, \theta_2,..., \theta_N\}$, wee define two random variables $X$ and $Y$ such that:
\begin{equation} X = \sum_{i=1}^{N} a_i \theta_i \end{equation}
and,
\begin{equation} Y = \sum_{i=1}^{N} b_i \theta_i \end{equation}
Clearly, $X$ and $Y$ are correlated and in the large limit of $N$ each follows a normal distribution with means and variances depending on the fixed sets $A$ and $B$.
What is the joint distribution of $X$ and $Y$ ?
I know that if two variables, say $U$ and $V$, are jointly gaussian, then any linear combination $W = \alpha U + \beta V$ follows a gaussian distribution BUT the converse need not to be true.
However, what is funny with this question is that $X$ and $Y$ is marginally gaussian in the limit of large enough $N$ and the same holds for all linear combinations of the two. If I were to perform a simulation and draw the heatmap the joint PDF of $X$ and $Y$, I would get some shape that is not symmetrical with respect to the mean values of $X$ and $Y$ (see figure below where the red dotted lines represents the means of $X$ and $Y$)
Hence, the joint PDF of $X$ and $Y$ cannot be gaussian.
Any references or answer to this question would be gladly appreciated.
We will avoid the trivial case where $A$ and $B$ are identical
Note: This has been cross-posted on "cross-validated"