Could anybody show me how @Rob Hyndman calculates the variance of $\hat{y}$ in the following link
Obtaining a formula for prediction limits in a linear model :
EDIT: Basically I don't understand how come $X^*(X'X)^{-1}X'$ is not squared as well: $Var(y_{pred})=Var(\hat{y}+\epsilon)=(X^*\hat{\beta})^2\sigma^2+\sigma^2$
Best Answer
The link you provided has a small typo: "...and its variance by..." The variance of the fitted value is not what that expression is. It's the mean square prediction error (or MSE), which is strictly larger.
If $\operatorname{Var}(\mathbf{e}) = \sigma^2 I$, and $\mathbf{X}$ is of full rank, then \begin{align*} \operatorname{Var}(\hat{y}) &= \operatorname{Var}\left[\mathbf{X}^*\left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\mathbf{Y}\right] \tag{you}\\ &= \mathbf{X}^*\left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\operatorname{Var}\left[\mathbf{Y}\right]\left[\mathbf{X}^*\left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\right]'\tag{me} \\ &= \sigma^2 \mathbf{X}^*\left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\mathbf{X}\left(\mathbf{X}'\mathbf{X}\right)^{-1} \mathbf{X}^{*'} \tag{*}\\ &= \sigma^2 \mathbf{X}^*\left(\mathbf{X}'\mathbf{X}\right)^{-1} \mathbf{X}^{*'} . \end{align*} Also, $$ \operatorname{MSPE}(y_{\text{pred}}) = \operatorname{Var}(y_{\text{pred}} - \hat{y}) = \operatorname{Var}(y_{\text{pred}}) + \operatorname{Var}(\hat{y}) $$ because of independence (or uncorrelated-ness) of future data, and how the prediction errors have mean zero.