Variance of square of the mean of Bernoulli variable 
Let $X_1,\dots,X_n$ be independent Bernoulli variables, with probability of success $p_i,\ldots,p_n$ and let $Y_n =\dfrac1n\sum\limits^n_{i=1} (X_i - p_i )$. Find $\operatorname{Var}(Y_n^2)$.

What I know is that
$$E(Y_n) = 0, \quad   \operatorname{Var}(Y_n) = \frac1{n^2}\sum^n_{i=1}p_i(1-p_i)$$ and I know I can use the formula $$\operatorname{Var}(Y_n^2) = E(Y_n^4) - (E(Y_n^2))^2$$ but I have no idea how to calculate the $E(Y_n^4)$.
 A: Applying the binomial theorem gives \begin{align}n^4\Bbb E[Y^4]&=\Bbb E\left[\left(\sum_{i=1}^n(X_i-p_i)\right)^4\right]\\&=\Bbb E\chi^4-4\rho\Bbb E\chi^3+6\rho^2\Bbb E\chi^2-4\rho^3\Bbb E\chi+\rho^4\end{align} where $\chi^\bullet=\left(\sum\limits_{i=1}^nX_i\right)^\bullet$ and $\rho^\bullet=\left(\sum\limits_{i=1}^np_i\right)^\bullet$. Note that since $\{X_i\}$ are independent Bernoulli random variables, we have $X_i^t=X_i$ for any positive integer $t$ so that \begin{align}\Bbb E\chi&=\rho\\\Bbb E\chi^2&=\Bbb E\left[\sum_{i=1}^nX_i^2+2\sum_{\text{cyc}}X_iX_j\right]=\rho+2\sum_{\text{cyc}}p_ip_j\\\Bbb E\chi^3&=\Bbb E\left[\sum_{i=1}^nX_i^3+3\sum_{\text{cyc}}X_i^2X_j+6\sum_{\text{cyc}}X_iX_jX_k\right]=\rho+3\sum_{\text{cyc}}p_ip_j+6\sum_{\text{cyc}}p_ip_jp_k\\\Bbb E\chi^4&=\Bbb E\left[\sum_{i=1}^nX_i^4+4\sum_{\text{cyc}}X_i^3X_j+6\sum_{\text{cyc}}X_i^2X_j^2+12\sum_{\text{cyc}}X_i^2X_jX_k+24\sum_{\text{cyc}}X_iX_jX_kX_\ell\right]\\&=\rho+10\sum_{\text{cyc}}p_ip_j+12\sum_{\text{cyc}}p_ip_jp_k+24\sum_{\text{cyc}}p_ip_jp_kp_\ell\end{align} which can be substituted back into $\Bbb E[Y^4]$.
