Serial Correlation - Fixed Effects Estimation Proof Assume we have a standard panel model. $y_{it} = \alpha_{I} + \beta{x_{it}} + u_{it}$ Also assume $Cov(u_{is},u_{it}) = 0$. 
I would like to show the fixed effects estimation strategy entails $Cov(\ddot{u}_{it}, \ddot{u}_{is}|X) \neq 0$. In words, we need to cluster for serial correlation. 
Here's how far I can get into this proof. 
\begin{align*}
Cov(\ddot{u}_{i,s},\ddot{u}_{i,t}) &=  Cov(u_{is} - \bar{u}_{i}, u_{it} - \bar{u}_{i}|X) \\
&= E[(u_{is} - \bar{u}_{i})(u_{it} - \bar{u}_{i})|X] - \underbrace{E[u_{is} - \bar{u}_{i}|X]E[u_{it} - \bar{u}_{i}|X]}_\text{Zero because $E[\ddot{u_{i}}|X] =0$}\\
&= E[(u_{is} - \bar{u}_{i})(u_{it} - \bar{u}_{i})|X] \\
&= E[u_{is}u_{it} - u_{is}\bar{u}_{i} - u_{it}\bar{u}_{i} + \bar{u}^{2}_{i}|X]  \\
&= E[u_{is}u_{it}|X] - E[u_{is}\bar{u}_{i}|X] - E[u_{it}\bar{u}_{i}|X] + E[\bar{u}^{2}_{i}|X]
\end{align*}
The first term after the last equality is zero because we assume $Cov(u_{is},u_{it}) = 0$. 
Yet from here I can't solve it any further. I suspect the result is $-\frac{\sigma_{u}^{2}}{T}$. Let me know how I should see this. 
Edit: Yes, I am assuming $Var(u_{it}|X) = \sigma_{u}^{2} \; \forall \; i$.
 A: Let $u_i=(u_{it},...,u_{iT})^\top$ have zero mean then
$$Var(u_i) = \mathbb E[u_iu_i^\top] = \sigma^2I_T$$
Let $D$ be the first difference matrix and do first differences to get transformed model
$$Dy_i = DX_i\beta + Du_i$$
Then see that the variance is 
$$Var(Du_i) = D\mathbb E[u_iu_i^\top]D^\top=\sigma^2 DD^\top$$
Since $DD^\top$ is known you can use generalized least squares known to be optimal. Applying this to the transformed model result in the estimator
$$\hat \beta = \left(\sum_i X_i^\top D^\top(DD^\top)^{-1}DX_i^\top\right)^{-1}\sum_i X_i^\top D^\top(DD^\top)^{-1}Dy_i$$
but the matrix $ D\top(DD^\top)^{-1}D = Q$ where $Q$ is the demeaning matrix.
The covariance matrix $Var(\ddot u_i) = Var(Qu_i) = Q \mathbb E[u_iu_i^\top]Q^\top = \sigma^2_u QQ^\top = \sigma^2_u Q$ due to symmetry and idempotency of $Q$ matrix. 
It follows that diagonal terms are given as
$$Var(u_{it}) = \sigma_u^2 (1 - 1/T)$$ and that off-diagonal terms are given as
$$Cov(u_{it},u_{is}) = -\sigma_u^2/T$$
knowing the structure of the demeaning matrix here shown in a $3 \times 3$ version 
$$Q := \begin{pmatrix} 1-1/T & -1/T & -1/T \\ -1/T & 1 - 1/T & -1/T \\ -1/T & -1/T & 1 - 1/T\end{pmatrix}$$
Witout matrix algebra the off-diagonal terms $s\not=t$ can be found by seeing that
$$Cov(\ddot u_{it},\ddot u_{is}) = \mathbb E[\ddot u_{it}\ddot u_{is}] = \mathbb E[(u_{it} - \bar u_i)(u_{is} - \bar u_i)]$$
first identity follows by definition of covariance and use of the fact that $\mathbb E[\ddot u_{it}]=0$ second from defintion of $\ddot u_{it}$. Continuing 
$$= \mathbb E[u_{it}u_{is}] - \mathbb E[u_{it}\bar u_{i}] - \mathbb E[\bar u_{i}u_{is}] + \mathbb E[\bar u_i \bar u_i]$$
First term is $0$ by assumption, second and third term are $-\sigma_u^2/T$ because $u_{it}$ is uncorrelated with all terms in the average $\bar u_i$ except one which is $u_{it}$ hence $\mathbb E[u_{it}\bar u_i] = \frac{1}{T} \mathbb E[u_{it}u_{it}] = \frac{\sigma_u^2}{T}$. The final term $\mathbb E[\bar u_i \bar u_i] = \frac{1}{T^2} \mathbb E[(u_{i1} + ...+ u_{iT})(u_{i1} + u_{iT})] = \frac{T}{T^2} \sigma_u^2 = \frac{1}{T}\sigma^2_u.$ Insert these and get the result as stated above. 
However this serial correlation is of little consequence for standard errors because
$$\sqrt N (\hat \beta - \beta) = \left(\frac{1}{N} \sum_i \ddot X_i^\top \ddot X_i\right)^{-1} \left(\frac{1}{\sqrt N}\ddot X_i^\top \ddot u_i\right),$$
and $\ddot X_i^\top \ddot u_i = X_i^\top Q^\top Q u_i = X_i^\top Q^\top u_i = \ddot X u_i$ showing that demeaning of errors can be removed such that
$$\sqrt N (\hat \beta - \beta) = \left(\frac{1}{N} \sum_i \ddot X_i^\top \ddot X_i\right)^{-1} \left(\frac{1}{\sqrt N}\ddot X_i^\top  u\right),$$
implying that the asymptotic variance is $\sigma^2_u \mathbb E[\ddot X_i^\top \ddot X_i ]$ under the assumption that $\mathbb E[u_iu_i^\top\lvert \ddot X_i] = \sigma_u^2 I_T$ as was assumed.
The standard errors are not clustered because the above removal of the $\ddot u$ in the expression for the asymptotic covariance implies that there is no need for estimation of $\hat \Omega = \frac{1}{N} \sum_i \hat{\ddot u}\hat{\ddot u}^\top$. As they say in the reference you post "In the setting without clustering, the key assumption is that $\Omega$ is diagonal." And although $\Omega := Var(Qu)$ is not diagonal this has no consequence for the asymptotic variance where the variance $Var(u)=\sigma^2I_T$ which is diagonal  is used. Offcourse if you want to calculate robust standard errors the story is different, then there is clustering by individual.
