Calculating the MSE for assessment Let $X_1, \ldots, X_n \sim \mathcal{N}(\mu, \sigma^2)$ be the sample, when $\mu$, $\sigma$ are unknown.
We suggest assessment for $\sigma^2$: 
$$S^2 = \frac{\displaystyle\sum_{i=1}^n (X_i - \bar{X})^2}{n-1}$$
Now, I know this assesment is unbiased, so for calculating the MSE I have to calculate the expectation.
Any suggestions for this? I'm pretty stuck in this part.
 A: I assume that MSE means mean-squared error, so I interpret this question as a request to show how to calculate the variance of the estimator $S^2$.  I'd organise the algebra for this using the Kronecker delta $\delta _{ij}$, so that 
$$
\begin{eqnarray*}
X_{i}-\bar{X} &=&\sum_{j=1}^{n}X_{j}\left( \delta _{ij}-\frac{1}{n}\right) 
\text{  .} \\
E\left( X_{i}X_{j}\right)  &=&\sigma ^{2}\delta _{ij} \\
E\left( X_{i}X_{j}X_{k}X_{l}\right)  &=&\sigma ^{4}\left( \delta _{ij}\delta
_{kl}+\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}+3\delta _{ij}\delta
_{ik}\delta _{il}\right) \text{  .}
\end{eqnarray*}
$$
Then the variance of the estimator $S^2$ is
$$
\begin{eqnarray*}
E\left( \left( S^{2}-\sigma ^{2}\right) ^{2}\right)  &=&E\left( \left( \frac{%
\sum_{i=1}^{n}\left( \sum_{j=1}^{n}X_{j}\left( \delta _{ij}-\frac{1}{n}%
\right) \right) ^{2}}{n-1}-\sigma ^{2}\right) ^{2}\right)  \\
&=&E\left( \left( \frac{\sum_{i,j,k=1}^{n}X_{j}X_{k}\left( \delta _{ij}-%
\frac{1}{n}\right) \left( \delta _{ik}-\frac{1}{n}\right) }{n-1}-\sigma
^{2}\right) ^{2}\right)
\end{eqnarray*}
$$
The equations are a bit to long to write down, but the biggest term is
$$
E\left( \frac{%
\sum_{i_{1},j_{1},k_{1},i_{2},j_{2},k_{2}=1}^{n}X_{j_{1}}X_{k_{1}}X_{j_{2}}X_{k_{2}}\left( \delta _{i_{1}j_{1}}-%
\frac{1}{n}\right) \left( \delta _{i_{1}k_{1}}-\frac{1}{n}\right) \left(
\delta _{i_{2}j_{2}}-\frac{1}{n}\right) \left( \delta _{i_{2}k_{2}}-\frac{1}{%
n}\right) }{\left( n-1\right) ^{2}}\right)
$$
and it's just a question of applying the equation for $E\left( X_{i}X_{j}X_{k}X_{l}\right)$ and working through the algebra.
A: Here is the solution using the mathStatica add-on to Mathematica:
$S^2$ is also known as the  $2^{\text{nd}}$ h-statistic (unbiased estimator of the  $2^{\text{nd}}$ central moment). Since $S^2$ is an unbiased estimator of population variance, $\text{MSE}[S^2]$ = $\text{Var}(S^2)$. 
So we seek the $2^{\text{nd}}$ Central Moment (i.e. Variance) of the $2^{\text{nd}}$ h-statistic, so the answer is a one-liner:
MSE = CentralMomentToCentral[2, HStatistic[2][[2]] ] 


where $\mu_2$ and $\mu_4$ denote the  $2^{\text{nd}}$ and  $4^{\text{th}}$ central moments of the population. In the special case of the Normal distribution, just substitute in for $\mu_2$ and $\mu_4$:

