Derive $E(MS_A)$ in design of experiments 
Prove that
$$
E(MS_A) = E\left[\frac{SS_A}{a-1}\right] = \sigma^2 + \frac{bn \sum_{i=1}^{a} \tau_i^2}{a-1}.
$$

Context. In Design and Analysis of Experiments, suppose you're doing a factorial design with 2 factors (2-way ANOVA). Call factors A (with $a$ levels) and B (with $b$ levels), the sample size is $n$. The effects model is then given by
$$
Y_{ijk} = \mu + \tau_i + \beta_j + (\tau \beta)_{ij} + \varepsilon_{ijk}, \quad \varepsilon_{ijk} \sim N(0, \sigma^2)
$$
Attempt. I've been trying to prove this result for over a couple hours. The difficulty is dealing with $SS_A$. In particular,
$$
SS_A = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} \left( \bar{Y}_{i..} - \bar{Y}_{...} \right)^2
$$
 A: Suppose the fixed effects model is
$$y_{ijk}=\mu+\tau_i+\beta_j+\gamma_{ij}+\varepsilon_{ijk}\quad;\small \,i=1,\ldots,a,j=1,\ldots,b,k=1,\ldots,n \tag{$\star$}$$
Here $\mu$ is a general effect; $\tau_i,\beta_j,\gamma_{ij}$ are additional effects subject to $$\sum_i \tau_i=\sum_j \beta_j=\sum_i \gamma_{ij}=\sum_j \gamma_{ij}=0$$ These constraints are necessary for unique estimation of all parameters in $(\star)$. And $\varepsilon_{ijk}\stackrel{\text{i.i.d}}\sim N(0,\sigma^2)$ of course is the random error.
Then
$$\overline y_{i\cdot\cdot}=\mu+\tau_i+\overline\varepsilon_{i\cdot\cdot}$$
and $$\overline y_{\cdot\cdot\cdot}=\mu+\overline \varepsilon_{\cdot\cdot\cdot}$$
Therefore,
$$\sum_{i=1}^a (\overline y_{i\cdot\cdot}-\overline y_{\cdot\cdot\cdot})^2=\sum_{i=1}^a (u_i-\overline u)^2\,,$$
where $$u_i=\tau_i+\overline\varepsilon_{i\cdot\cdot}$$
Note that $u_i$'s are independent normal variables:
$$u_i \stackrel{\text{ind}}\sim N\left(\tau_i,\frac{\sigma^2}{bn}\right)$$
Hence,
\begin{align}
E\left[SSA\right]&=E\left[bn\sum_{i=1}^a (\overline y_{i\cdot\cdot}-\overline y_{\cdot\cdot\cdot})^2\right]
\\&=E\left[bn\sum_{i=1}^a (u_i-\overline u)^2\right]
\\&=bn\left\{\sum_{i=1}^a E\left[u_i^2\right]-a\,E\left[\overline u^2\right]\right\}
\end{align}
Simplifying this is straightforward using the distribution of the $u_i$'s.
Another way to do this is to use general theory of quadratic forms to conclude that $SSA/\sigma^2$ has a non-central $\chi^2$ distribution with degrees of freedom $a-1$ and non-centrality parameter $\frac{bn}{\sigma^2}\sum_i \tau_i^2$.
This would directly give $$E\left[\frac{SSA}{\sigma^2}\right]=a-1+\frac{bn}{\sigma^2}\sum_{i=1}^a \tau_i^2$$
