Finding the covariance of two random sums I am trying to derive the covariance of two sample sums.
Some notation and details:
$x_i$ and $y_i$ are numeric values of two characteristics corresponding to member i of a finite population of N members.
$w_i$ is a random indicator variable taking a value of 1 with a probability of n/N, and is otherwise 0.  All of the N $w_i$ values sum to n.
So essentially, I’m taking a sample of size n from a population of size N, and I’m separately summing up the x and the y values for that sample.  I want to find the covariance of these two sums, and I think I’m close but could use some help.
$$
Cov\left(\sum_{i=1}^{N}{w_ix_i},\sum_{i=1}^{N}{w_iy_i}\right)
$$

*** Edit
Following is my original post (with incorrect derivation) but scroll down for the correct answer including some R code confirming the answer.  For a derivation of the answer, see the extremely helpful post by whuber below.

$$
=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-E\left(\sum_{i=1}^{N}{w_ix_i}\right)\right)\left(\sum_{i=1}^{N}{w_iy_i}-E\left(\sum_{i=1}^{N}{w_iy_i}\right)\right)\right]
$$
$$
=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\sum_{i=1}^{N}{x_iE\left(w_i\right)}\right)\left(\sum_{i=1}^{N}{w_iy_i}-\sum_{i=1}^{N}{y_iE\left(w_i\right)}\right)\right]
$$
$$
=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\frac{n}{N}\sum_{i=1}^{N}x_i\right)\left(\sum_{i=1}^{N}{w_iy_i}-\frac{n}{N}\sum_{i=1}^{N}y_i\right)\right]
$$$$
=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-n\mu_x\right)\left(\sum_{i=1}^{N}{w_iy_i}-n\mu_y\right)\right]$$$$
=\ E\left[\sum_{i=1}^{N}{w_ix_i}\sum_{i=1}^{N}{w_iy_i}-\ n\mu_y\sum_{i=1}^{N}{w_ix_i}-\ n\mu_x\sum_{i=1}^{N}{w_iy_i}+n^2\mu_x\mu_y\right]$$$$
=\ E\left[\sum_{i=1}^{N}{w_ix_iw_iy_i}\right]+E\left[\sum_{i=1}^{N}\sum_{j\neq i}^{N}{w_ix_iw_jy_j}\right]-\ n\mu_y\frac{n}{N}\sum_{i=1}^{N}x_i-\ n\mu_x\frac{n}{N}\sum_{i=1}^{N}y_i+n^2\mu_x\mu_y$$$$
=\ \left(\sum_{i=1}^{N}{x_iy_i}\right)\left(\frac{n}{N}\right)+\left(\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}\right)\left(\frac{n}{N}\right)\left(\frac{n-1}{N-1}\right)-\ n^2\mu_y\mu_x-\ n^2\mu_x\mu_y+n^2\mu_x\mu_y$$$$
=\ n\mu_{xy\left(i=j\right)}+n\left(n-1\right)\mu_{xy\left(i\neq j\right)}-\ n^2\mu_y\mu_x$$$$
=\ n^2\left(\frac{\mu_{xy\left(i=j\right)}+\left(n-1\right)\mu_{xy\left(i\neq j\right)}}{n}-\mu_x\mu_y\right)$$$$
=\ n^2(\mu_{xy}-\mu_x\mu_y)$$
$$=\ n^2Cov(x_i,y_i)$$
The above looks right to me, but I must be doing something wrong.  When I simulate the covariance in R, I find that it is not the equation above.  However, it seems (consistent with simulation) I can obtain the correct covariance by applying a finite population correction and by adding an n to the denominator of my result, as shown below.
$$\frac{N-n}{N}nCov(x_i,y_i)$$

*** End of original post:
Following is the answer and some R code that "confirms" it.

$$\frac{n(N-n)}{N-1}Cov(x_i,y_i)$$
(where Cov() represents the population covariance)
R Code with results that are consistent with this answer:
library(mvtnorm)
library(tidyverse)

N <- 5
n <- 3

val1 <- rmvnorm(n = N, mean = c(50, 100), sigma = matrix(c(15^2, 12^2, 12^2, 15^2), nrow = 2))
colnames(val1) <- c("var1","var2")

set.seed(83442)

numSims <- 100000
covList1 <- as.list(1:numSims)
for(i in 1L:numSims) {
  covList1[[i]] <- val1[sample.int(N, n),] %>%
    as_tibble() %>%
    summarize(var1 = sum(var1), var2 = sum(var2))
}

covDF1 <- do.call("bind_rows", args = covList1)

(cov1 <- cov(covDF1$var1, covDF1$var2)) # Covariance from simulations: 195.7401
((N-n)/(N-1))*n*cov(val1[, 1], val1[, 2])*(N-1)/N # Covariance with N-1 denominator: 196.2536
# Note that the (N-1)/N at the end adjusts for the fact that we need a population (not sample) covariance. 

The results of this comparison improve with more simulations.
 A: To help make the ideas clear, I will use capital letters for random variables.
Everything follows from the restriction $\sum W_i=n,$ because that implies this sum has zero variance. Since each $W_i$ is a Bernoulli variable,
$$\operatorname{Var}(W_i) = \frac{n(N-n)}{N^2}.$$
Computing the variance of the sum and assuming, as is the case with simple random sampling, that for $i\ne j$ $\operatorname{Cov}(W_i,W_j)$ does not depend on $i$ or $j,$ we find
$$\begin{aligned}
0 &= \operatorname{Var}\left(\sum_{i=1}^N W_i\right) \\
&=\sum_{i=1}^N \operatorname{Var}\left(W_i\right) + \sum_{i\ne j}^N \operatorname{Cov}(W_i,W_j) \\
&= N\frac{n(N-n)}{N^2} + N(N-1)\operatorname{Cov}\left(W_1,W_2\right),
\end{aligned}$$
enabling us to solve for the covariance as

$$\operatorname{Cov}\left(W_i,W_j\right) = \operatorname{Cov}\left(W_1,W_2\right) = -\frac{n(N-n)}{N^2(N-1)}.$$

Consequently, assuming $N\gt 1,$ for fixed coefficients $(x_i)$ and $(y_i)$ and writing $\bar x = \sum x_i/N,$ $\bar y = \sum y_i/N,$ and $\overline{xy}=\sum_{i}x_iy_i/N,$ we find
$$\begin{aligned}
\operatorname{Cov}\left(\sum_{i=1}^N x_iW_i, \sum_{j=1}^N y_jW_j\right) &= \sum_{i=1}^N x_iy_i \operatorname{Var}\left(W_i\right) + \sum_{i\ne j}^N x_iy_j\operatorname{Cov}\left(W_i,W_j\right) \\
&= \frac{n(N-n)}{N^2}\sum_{i=1}^Nx_iy_i - \frac{n(N-n)}{N^2(N-1)}\sum_{i\ne j}^N x_iy_j \\
&= \frac{n(N-n)}{N}\overline{xy} - \frac{n(N-n)}{N-1} \bar{x}\bar{y} + \frac{n(N-n)}{N(N-1)}\overline{xy}\\
&= \frac{n(N-n)}{N-1}\left(\overline{xy} - \bar{x}\bar{y} \right).
\end{aligned}$$
(When $N=1$ the double sum doesn't appear and the result easily reduces to $0.$)
If we draw one of the $(x_i,y_i)$ randomly and equiprobably from all $N$ of these paired values, we have a bivariate random variable $(X,Y),$ enabling the result to be written

$$\operatorname{Cov}\left(\sum_{i=1}^N x_iW_i, \sum_{j=1}^N y_jW_j\right) = \frac{n(N-n)}{N-1} \operatorname{Cov}(X,Y).$$


I was tempted to check this result with simulation, but elected to use an exhaustive enumeration of all the possible samples instead, of which there are $\binom{N}{n}.$  For small $N$ this is feasible and gives precise results.  The output computes the covariance of the weighted sums in three ways: using the formula in terms of $\overline{xy}-\bar{x}\bar{y},$ the formula in terms of $\operatorname{Cov}(X,Y),$ and--this is the verification--the population covariance of all possible sample sums.
An example of its output for $N=20,$ $n=15$ is
Direct Formula Covariance formula         Exhaustive 
     -3.035239          -3.035239          -3.035239 

demonstrating agreement in this case.
Here's the R code.
#
# Create *any* bivariate population you like.
#
N <- 20
# set.seed(17)
x <- rnorm(N)
y <- rnorm(N) - x
#
# Specify the sample size.
#
n <- 15
if(choose(N, n) > 1e6) stop("Are you sure you want to do this?", call.=FALSE)
# 
# Compute the distribution of the sample sum.
#
W <- combn(1:N, n)
wx <- apply(W, 2, function(w) sum(x[w]))
wy <- apply(W, 2, function(w) sum(y[w]))
#
# Compare various formulae.
#
c(`Direct Formula` = n * (N-n) / (N-1) * (mean(x*y) - mean(x)*mean(y)),
`Covariance formula`=n * (N-n) / N * cov(x, y),
Exhaustive = cov(wx, wy)*(1 - 1/length(wx)))

# plot(wx, wy) # Can be interesting...
```

A: Based on some of the methods whuber used in his answer, I decided to derive this covariance again, but this time in the way that I originally attempted the derivation (starting by showing the covariance as an expectation and going from there).  I get the right answer now:
$$Cov\left(\sum_{i=1}^{N}{w_ix_i},\sum_{i=1}^{N}{w_iy_i}\right)$$
$$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-E\left(\sum_{i=1}^{N}{w_ix_i}\right)\right)\left(\sum_{i=1}^{N}{w_iy_i}-E\left(\sum_{i=1}^{N}{w_iy_i}\right)\right)\right]$$
$$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\sum_{i=1}^{N}{x_iE\left(w_i\right)}\right)\left(\sum_{i=1}^{N}{w_iy_i}-\sum_{i=1}^{N}{y_iE\left(w_i\right)}\right)\right]$$
$$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\frac{n}{N}\sum_{i=1}^{N}x_i\right)\left(\sum_{i=1}^{N}{w_iy_i}-\frac{n}{N}\sum_{i=1}^{N}y_i\right)\right]$$
$$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-n\mu_x\right)\left(\sum_{i=1}^{N}{w_iy_i}-n\mu_y\right)\right]$$
$$=\ E\left[\sum_{i=1}^{N}{w_ix_i}\sum_{i=1}^{N}{w_iy_i}-\ n\mu_y\sum_{i=1}^{N}{w_ix_i}-\ n\mu_x\sum_{i=1}^{N}{w_iy_i}+n^2\mu_x\mu_y\right]$$
$$=\ E\left[\sum_{i=1}^{N}{w_ix_iw_iy_i}\right]+E\left[\sum_{i=1}^{N}\sum_{j\neq i}^{N}{w_ix_iw_jy_j}\right]-\ n\mu_y\frac{n}{N}\sum_{i=1}^{N}x_i-\ n\mu_x\frac{n}{N}\sum_{i=1}^{N}y_i+n^2\mu_x\mu_y$$
$$=\ \left(\sum_{i=1}^{N}{x_iy_i}\right)\left(\frac{n}{N}\right)+\left(\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}\right)\left(\frac{n}{N}\right)\left(\frac{n-1}{N-1}\right)-\ n^2\mu_y\mu_x-\ n^2\mu_x\mu_y+n^2\mu_x\mu_y$$
$$=\left(\frac{n\left(N-1\right)}{N\left(N-1\right)}\right)\left(\sum_{i=1}^{N}{x_iy_i}\right)+\left(\frac{n\left(n-1\right)}{N\left(N-1\right)}\right)\left(\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}\right)-n^2\mu_x\mu_y$$
$$=\frac{n\left(N-n\right)}{N\left(N-1\right)}\sum_{i=1}^{N}{x_iy_i}+\frac{n\left(n-1\right)}{N\left(N-1\right)}\sum_{i=1}^{N}{x_iy_i}+\frac{n\left(n-1\right)}{N\left(N-1\right)}\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}-n^2\mu_x\mu_y$$
$$=\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_{xy}+\frac{Nn\left(n-1\right)}{N^2\left(N-1\right)}\sum_{i=1}^{N}x_i\sum_{j=i}^{N}y_j-n^2\mu_x\mu_y$$
$$=\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_{xy}+\frac{Nn\left(n-1\right)}{\left(N-1\right)}\mu_x\mu_y-\frac{n^2\left(N-1\right)}{\left(N-1\right)}\mu_x\mu_y$$
$$=\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_{xy}-\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_x\mu_y$$
$$=\frac{n\left(N-n\right)}{\left(N-1\right)}\left(\mu_{xy}-\mu_x\mu_y\right)$$
$$=\frac{n\left(N-n\right)}{\left(N-1\right)}Cov(x_i,y_i)$$
