Suppose $X_i, Y_i, Z_i$, $i = 1, ..., n$ are independent Bernoulli random variables, with $X_i \sim Bern(\theta)$, $Y_i \sim Bern(\eta)$, $Z_i \sim Bern(\pi)$, where $\theta, \eta, \pi \in (0, 1)$. Define $S = \sum_{i=1}^n X_i$, $T = \sum_{i=1}^n X_iY_i$, $U = \sum_{i=1}^n X_iZ_i$
Find $Cov(S, T)$ and $Cov(T, U)$.
$$Cov(S, T) = E(ST) - E(S)E(T)$$
$$Cov(S, T) = E(ST) - E(\sum_{i=1}^n X_i)E(\sum_{i=1}^n X_iY_i)$$
$$Cov(S, T) = E(ST) - \sum_{i=1}^n E(X_i)\sum_{i=1}^n E(X_iY_i)$$
$$Cov(S, T) = E(ST) - \sum_{i=1}^n \theta\sum_{i=1}^n E(X_i)E(Y_i)$$
$$Cov(S, T) = E(ST) - \sum_{i=1}^n \theta\sum_{i=1}^n \theta\eta$$ $$Cov(S, T) = E(ST) - n\theta n\theta\eta$$ $$Cov(S, T) = E(ST) - n^2\theta^2\eta$$
I'm struggling to calculate E(ST) and I'm not sure if my thought process for calculating E(S) and E(T) is correct. Any advice? I'm assuming the same approach can be taken to find Cov(T, U).
Hint: Write out $ST$ as a double sum: $$ST = \left(\sum_i X_i \right)\left(\sum_j X_jY_j\right) = \sum_i \sum_j X_iX_jY_j$$ where both $i$ and $j$ range from $1$ to $n$. You will need to consider the case $i=j$ separately from the case $i\ne j$, because the variables $X_i$ and $X_j$ are independent only when $i\ne j$. When $i=j$ you have $X_i^2$, whose expectation you can deduce from $\operatorname{Var}(X_i)$ and $E(X_i)$. (Meanwhile, $Y_j$ is independent of any combination of $X$'s.)
You can handle $TU$ similarly.