Correlation Coefficient for Hypergeometric-type Distribution Problem Statement: A box contains $N_1$ white balls, $N_2$ black balls, and $N_3$ red balls
$(N_1+N_2+N_3=N).$ A random sample of $n$ balls is selected from the box (without replacement).
Let $Y_1,Y_2,$ and $Y_3$ denote the number of white, black, and red balls, respectively, observed
in the sample. Find the correlation coefficient for $Y_1$ and $Y_2.$ (Let $p_i=N_i/N$ for
$i=1,2,3.$)
My Work So Far: Because the sampling is without replacement, the distribution is much like the hypergeometric distribution. In fact, we have
the joint distribution function as the multivariate hypergeometric distribution:
$$
p(y_1,y_2,y_3)
=\frac{\displaystyle\binom{N_1}{y_1}\binom{N_2}{y_2}\binom{N_3}{y_3}}
{\displaystyle\binom{N}{n}}.
$$
Since we are not concerned about $Y_3,$ we can simplify a bit:
$$
p(y_1,y_2)
=\frac{\displaystyle\binom{N_1}{y_1}\binom{N_2}{y_2}\binom{N-(N_1+N_2)}{n-(y_1+y_2)}}
{\displaystyle\binom{N}{n}}.
$$
This is as far as I can get. The sums involved in calculating even something like $E(Y_1Y_2)$ are beyond my abilities (and even beyond the abilities of Mathematica). I know that the answer is
$$\operatorname{Cov}(Y_1,Y_2)=-n\,\frac{N-n}{N-1}\,\frac{N_1}{N}\,\frac{N_2}{N}.$$
My Question: How do I move on from here?
Thanks for your time!
 A: Let the white balls form a set $\mathcal W$ of size $N_1$ and the black balls form a set $\mathcal B$ of size $N_2,$ both of which are disjoint subsets of the set of all $N$ balls, $\mathcal U.$
For any ball $i\in \mathcal U,$ let $X_i$ be the indicator that ball $i$ is included in the sample $\mathcal S:$ that is, it equals $1$ when $i\in S$ and otherwise equals $0.$  The process of sampling makes these $X_i$ into random variables, because associated with any sample $\mathcal S$ is the vector of values of all the $X_i$ corresponding to $\mathcal S.$
Indicators enjoy nice connections with the underlying probabilities.  In particular, notice that since $X_i^2 = X_i,$
$$E[X_i^2] = E[X_i] = \Pr(X_i=1)1 + \Pr(X_i=0)0 = \Pr(X_i=1) = \Pr(i\in\mathcal{S}).$$
Let's pause to work out some important probabilities. To say the sample has size $n$ is equivalent to saying the sum of the indicators is $n:$
$$n = \sum_{i\in\mathcal U}X_i.$$
Since in simple random sampling all balls have equal chances $\pi_i = p$ (say) of being in the sample, by taking expectations we deduce
$$n = E\left[\sum_{i\in\mathcal U}X_i\right] = \sum_{i\in\mathcal{U}} E\left[X_i\right] = \sum_{i\in\mathcal{U}} \pi_i = \sum_{i\in\mathcal{U}} p = Np,$$
whence

$$\Pr(i\in\mathcal S) = \pi_i = p = \frac{n}{N}.$$

Similarly, given two distinct balls $i$ and $j\ne i,$ the chance that both $i$ and $j$ are in the sample (which I will write as $\pi_{ij}=q,$ which also does not depend on $i$ or $j$) can be found by taking expectations of the product $n\times n$ expressed as sums of indicators:
$$\begin{aligned}
n^2 &= E[n^2] = E\left[\sum_{i\in\mathcal{U}} X_i\, \sum_{j\in\mathcal{U}}X_j\right]\\
&= \sum_{i,\,j} E[X_iX_j]\\
&= \sum_i E[X_i^2] + \sum_{i\ne j}E[X_iX_j]\\
&= \sum_i E[X_i] + \sum_{i\ne j}\pi_{ij}\\
& = Np + N(N-1)q\\
&= n + N(N-1)q,
\end{aligned} $$
whence

$$\Pr(i\ne j\in\mathcal S) = \pi_{ij} = q = \frac{n^2-n}{N(N-1)} = \frac{\binom{n}{2}}{\binom{N}{2}},$$

the same result one would get from a combinatorial argument.
Returning to the question, we now have enough information to compute variances and covariances.  For example, with $i\ne j,$
$$\operatorname{Cov}(X_i,X_j) = E[X_iX_j] - E[X_i]E[X_j] = q - p^2 = -\frac{n(N-n)}{N^2(N-1)}$$
and, similarly,
$$\operatorname{Var}(X_i) = E[X_i^2] - E[X_i]E[X_j] = p - p^2 = \frac{n(N-n)}{N^2}.$$
Finally, we may compute variances and correlations of the $Y$'s by expressing them as suitable sums of indicators and applying the bilinear property of covariance.  For instance,

$$\begin{aligned}
\operatorname{Cov}(Y_1,Y_2) &= \operatorname{Cov}\left(\sum_{i\in\mathcal W}X_i,\sum_{j\in\mathcal{B}} X_j\right)\\&= \sum_{i\in\mathcal{W},\,j\in\mathcal{B}}\operatorname{Cov}(X_i,X_j)\\&= N_1\,N_2\left(-\frac{n(N-n)}{N^2(N-1)}\right)\end{aligned}$$

because $\mathcal W$ and $\mathcal B$ are disjoint.  This agrees with the formula quoted in the question.

You have now seen examples of all the techniques needed to complete the calculation of the correlation coefficient of $Y_1$ and $Y_2.$  (When computing $\operatorname{Var}(Y_k),$ you will need to split the sum into separate sums of variances $\operatorname{Cov}(X_i,X_i)=\operatorname{Var}(X_i)$  and covariances $\operatorname{Cov}(X_i,X_j),$ $j\ne i,$ much as I did in the calculation of $E[n^2],$ so make sure you're comfortable manipulating such double sums.)
In the spirit of self-study questions I will stop short of a complete answer.
