Estimate of variance of sum of Bernoulli variables Let $X_1,\ldots,X_n$ be Bernoulli distributed with parameter $p$. We define $\overline X = \frac{1}{n}\sum^n_{i=1} X_i$ and consider a point estimate $T=\overline X(1- \overline X)$. I want to compute if the estimator is biased or equivalent, if $\mathbb{E}[T]=\sigma^2$. I computed this two different ways, with two different results:
Firstly: $\mathbb{E}[T]=\mathbb{E}[\overline X - \overline X^2]=\mathbb{E}[\overline X]-\mathbb{E}[\overline X^2]$.
For $\mathbb{E}[\overline X]$ is used the linearity of the mean value operator and that the expectation of a single Bernoulli variable is $p$. I then get $\mathbb{E}[\overline X]=\frac{1}{n}np=p$. For $\mathbb{E}[\overline X^2]$ I used that $\overline X \approx \frac{1}{n} \operatorname{Bin}(n,p)$ and looked up the second moment which is $n^2p^2+np(1-p)$. The first way would then yield
$\mathbb{E}[T]=\mathbb{E}[\overline X]-\mathbb{E}[\overline X^2]=p - \frac{1}{n^2}(n^2p^2+np(1-p))=p[(1-p)(1-\frac{1}{n})]$. Thus, the estimator is biased.
Secondly: I tried to calculate $\mathbb{E}[\overline X^2]$ in another way. I used that for each $X_i\approx \operatorname{Bernoulli}(p)$ the second moment is given by $\mathbb{E}[X_i^2]=p\cdot 1 + (1-p)\cdot 0=p$ and then used that 
$\mathbb{E}[\overline X^2]=\frac{1}{n^2}\mathbb{E}[\sum_{i,j=1}^nX_iX_j]=\frac{1}{n^2}\sum_{i=1}^{n^2} \mathbb{E}[X_i^2]=\frac{1}{n^2}n^2p=p$ which is obviously different from my frist method.
I am thinking it should have worked both ways, but I can't find the error.

The issue is when you write, in your second approach,
$$
\frac{1}{n^2}\mathbb{E}\left[\sum_{i,j=1}^nX_iX_j\right]=\frac{1}{n^2}\sum_{i=1}^{n^2} \mathbb{E}[X_i^2]
$$
as the equality is not true. Indeed, what you do have is
$$
\mathbb{E}[X_iX_j]= \begin{cases}\mathbb{E}[X_i^2] = \mathbb{E}[X_i] =p & i=j\\ \mathbb{E}[X_i]\mathbb{E}[X_j] = p^2& i\neq j \end{cases}
$$
here the case $i\neq j$ follows from independence.
