Show that the co-variance between $X_j$ and $X_k$ is $\frac {-pq}{N-1}, j \ne k$ An urn containing $pN$ white and $qN$ black balls, the total number of balls being $N$. Balls are drawn one by one without being returned to the urn until a certain number $n$ of balls is reached. 
Let $ X_i=
\begin{cases}
 1&\text{if the i$_{th}$ drawn ball is white }\ \\
 0&\text{if the i$_{th}$ drawn ball is black}\
\end{cases}
$ 
Show that the co-variance between  $X_j$ and $X_k$ is $\dfrac {-pq}{N-1}, j \ne k$
Attempt: The joint probability density function of $X_j X_k$ can be computed as :
$\begin{array}{|c|c|c|c|c|}
\hline
Y=X_jX_k& (1\times1=1) & (1\times0=0) & (0\times1=0) & (0\times0=0) \\ \hline
 P(Y)& p^2& pq&pq &q^2\\ \hline
\end{array}$
The co-varirance between $X_j \text{ and}\ X_k = E [ ~\{X_j - E(X_j) \} \{ X_k-E(X_k)\}~] = E(X_j X_k)-E(X_j)E(X_k)\\ = 1 \times p^2 - [p \times p ] = 0
$
Where could I be going wrong? Is there a conceptual error somewhere?
Thanks a lot for the help
 A: Your calculations assume draws with replacement and therefore yields uncorrelated $X_j$ and $X_k$. First of all, $E[X_j]$ is trivial and equal to $p$ as you've also noted. We need to find $E[X_jX_k]=P(X_j=1,X_k=1)$. 
Here, the key property is 
$$P(X_j=1,X_k=1)=P(X_1=1,X_2=1)=\frac{pN}{N}\frac{pN-1}{N-1}$$
(to see it, start with $P(X_1=1,X_3=1)=\sum_{x_2\in\{0,1\}}P(X_1=1,X_2=x_2,X_3=1)$)
Which yields $\operatorname{cov}(X_j,X_k)=p\frac{pN-1}{N-1}-p^2=\frac{-pq}{N-1}$
A: Visualization:
Suppose that there are $N=4$ balls in the urn---one white and three black.
We take $n=3$ balls from the urn without replacement. We are interested
in random variables $X_1,$ which takes the value 1 when the first ball
is white and 0 otherwise and $X_3,$ which takes the value 1 when the third
ball is white and 0 otherwise. Intuitively, it is clear that $X_1$ and $X_3$ are negatively associated because choosing a white ball on the first draw makes it impossible to get a white ball on the third.
Then in @gunes' Answer (+1), we have $E(X_1) = E(X_3) = p = .25,$
$E(X_1X_3) = (.25)(0/4) = 0,$ and $Cov(X_1,X_3) = 0 - (.25)^2 =
-0.0625.$
A simulation in R of a million such 3-ball draws gives results that agree
with the above to about 3 decimal places.
set.seed(1025)
urn = c(0,0,0,1)
N=10; n=3; m=10^6; x1 = x3 = numeric(m)
for(i in 1:m) {
  x = sample(urn, n)
  x1[i] = x[1];  x3[i] = x[3] }

e.1 = mean(x1);  e.3 = mean(x3);  e.13 = mean(x1*x3)
e.1;  e.3;  e.13
[1] 0.249513     # aprx p = 0.25
[1] 0.249795     # aprx p = 0.25
[1] 0

e.13 - e.1*e.3
[1] -0.0623271   # aprx cov = -0.0625 from above
cov(x1, x3)
[1] -0.06232716  # aprx cov = -0.0625 using R fcn

A key point is that $0 = P(X_1 = 1, X_3 = 1) \ne P(X_1 = 1)P(X_3 = 1) = 0.0625.$
mean(x1==1 & x3==1)
[1] 0
mean(x1==1)* mean(x3==1)
[1] 0.0623271     # aprx 0.0625

One way to show this in a plot is to choose 10,000 of the 1 million
points simulated, and jitter them (with random uniform noise) to avoid overplotting. Then the probability at a point is represented by the density of points in a square with the point (red) at its center.

The left-hand panel shows results from the simulation above. At right, is a similar plot from balls sampled with replacement, so that draws are independent.
