Solution (source code)

= Solution

Write
$$
Y=\mu+\varepsilon,
\qquad
Y^*=\mu+\varepsilon^*,
$$
where $\varepsilon$ and $\varepsilon^*$ are independent, centered, and have covariance matrix $\sigma^2I_n$. Then
$$
HY-Y^*=-(I-H)\mu+H\varepsilon-\varepsilon^*.
$$
The deterministic term is orthogonal in expectation to the two centered random terms, and the random terms are independent. Since the <hat matrix> is a symmetric idempotent projection of rank $p$,
$$
\mathbb E\|H\varepsilon\|^2
=\sigma^2\operatorname{tr}(H^TH)
=\sigma^2\operatorname{tr}H=p\sigma^2,
$$
while $\mathbb E\|\varepsilon^*\|^2=n\sigma^2$. Therefore
$$
\boxed{
\mathbb E\|HY-Y^*\|^2
=\|(I-H)\mu\|^2+(n+p)\sigma^2
}.
$$
In this <bias-variance decomposition for linear prediction>, $\|(I-H)\mu\|^2$ is squared model bias, $p\sigma^2$ is variance from fitting $p$ coefficients, and $n\sigma^2$ is irreducible noise in the future response. Enlarging the model tends to reduce the first term while increasing the fitted-model variance, which is the bias-variance tradeoff.