Showing unbiasedness of the variance estimator: $E(\hat \sigma^2)=\sigma^2$ There is an infinite population characterized by a probability model. A simple random sample of size $n$ is $X_1,X_2,...,X_n$. Moreoever, $X_i$ are i.i.d. according to the model. 
The population mean is $\mu = E(X_i)$. The sample mean is $\hat \mu=\frac {1}{n}\sum^n _{i=1}X_i$. $E(\hat \mu)=\mu$.
Variance of the estimator is $V(\hat \mu)=\frac {\sigma^2}{n}$, where ${\sigma^2}=V(X_i)$. Also, $\hat \sigma^2=\frac {1}{n-1}\sum^n_{i=1}(X_i-\hat \mu)^2$. 
I want to show that $E(\hat \sigma^2)=\sigma^2$ .
My confusion here is that if $X_i$ is random variable, does this mean the sample mean $\hat \mu$ is also a random variable? I also tried the fact that $\sum^n_{i=1}(X_i-\hat \mu)^2 = \sum^n_{i=1}[(X_i- \mu)+(\mu-\hat \mu)]^2$, but I'm stuck on how to proceed forward. 
 A: Expanding, we have
\begin{equation}
 \hat \sigma^2=\frac {1}{n-1}\sum^n_{i=1}(X_i-\hat \mu)^2
    =\frac {1}{n-1}\sum^n_{i=1} X_i^2 - \frac {2}{n-1}\hat{\mu} \sum^n_{i=1}X_i + \frac {n}{n-1}\hat{\mu}^2
\end{equation}
Taking expectations
\begin{equation}
 E \hat \sigma^2
 =\frac {1}{n-1}\sum^n_{i=1} E X_i^2 - \frac {2}{n-1} E (\hat{\mu} \sum^n_{i=1}X_i) + \frac {n}{n-1}E\hat{\mu}^2
\end{equation}
Notice that
\begin{align}
 E X_i^2 &= \sigma^2 +\mu^2\\
 E (\hat{\mu} \sum^n_{i=1}X_i) &= nE (\hat{\mu} \frac{1}{n} \sum^n_{i=1}X_i)=nE \hat{\mu}^2
\end{align}
So
\begin{equation}
 E \hat \sigma^2
 =\frac {n}{n-1}(\sigma^2 + \mu^2)- \frac {2n}{n-1} E \hat{\mu}^2 + \frac {n}{n-1}E\hat{\mu}^2
  =\frac {n}{n-1}(\sigma^2 + \mu^2- E \hat{\mu}^2 ) \tag{1}
\end{equation}
Notice that 
\begin{equation}
 E \hat{\mu}^2 
 =
 \frac{1}{n^2}
 E (\sum_{i=1}^n X_i)^2
 =
 \frac{1}{n^2}
 \sum_{i=1}^n \underbrace{E X_i^2}_{\mu^2 + \sigma^2} + \sum_{i \neq j} \underbrace{E(X_iX_j)}_{\mu^2}
\end{equation}
we get
\begin{equation}
 E \hat{\mu}^2
 =
 \frac{1}{n^2}
 (n (\mu^2 + \sigma^2) + (n^2-n)(\mu^2)) 
 =
 \frac{1}{n}\sigma^2 + \mu^2
\end{equation}
Replacing in $(1)$, we get
\begin{equation}
 E \hat \sigma^2
  =\frac {n}{n-1}(\sigma^2- \frac{1}{n}\sigma^2 )
  =
  \frac{n}{n-1}\frac{n-1}{n}\sigma^2 
  =
  \sigma^2
\end{equation}
