A few more steps of the Bias - Variance decomposition
Indeed, the full derivation is rarely given in textbooks as it involves a lot of uninspiring algebra. Here is a more complete derivation using notation from the book "Elements of Statistical Learning" on page 223
If we assume that $Y = f(X) + \epsilon$ and $E[\epsilon] = 0$ and $Var(\epsilon) = \sigma^2_\epsilon$ then we can derive the expression for the expected prediction error of a regression fit $\hat f(X)$ at an input $X = x_0$ using squared error loss
$$Err(x_0) = E[ (Y - \hat f(x_0) )^2 | X = x_0]$$
For notational simplicity let $\hat f(x_0) = \hat f$, $f(x_0) = f$ and recall that $E[f] = f$ and $E[Y] = f$
\begin{aligned}
E[ (Y - \hat f)^2 ] &= E[(Y - f + f - \hat f )^2]
\\
& = E[(y - f)^2] + E[(f - \hat f)^2] + 2 E[(f - \hat f)(y - f)]
\\
& = E[(f + \epsilon - f)^2] + E[(f - \hat f)^2] + 2E[fY - f^2 - \hat f Y + \hat f f]
\\
& = E[\epsilon^2] + E[(f - \hat f)^2] + 2( f^2 - f^2 - f E[\hat f] + f E[\hat f] )
\\
& = \sigma^2_\epsilon + E[(f - \hat f)^2] + 0
\end{aligned}
For the term $E[(f - \hat f)^2]$ we can use a similar trick as above, adding and subtracting $E[\hat f]$ to get
\begin{aligned}
E[(f - \hat f)^2] & = E[(f + E[\hat f] - E[\hat f] - \hat f)^2]
\\
& = E \left[ f - E[\hat f] \right]^2 + E\left[ \hat f - E[ \hat f] \right]^2
\\
& = \left[ f - E[\hat f] \right]^2 + E\left[ \hat f - E[ \hat f] \right]^2
\\
& = Bias^2[\hat f] + Var[\hat f]
\end{aligned}
Putting it together
$$E[ (Y - \hat f)^2 ] = \sigma^2_\epsilon + Bias^2[\hat f] + Var[\hat f] $$
Some comments on why $E[\hat f Y] = f E[\hat f]$
Taken from Alecos Papadopoulos here
Recall that $\hat f$ is the predictor we have constructed based on the $m$ data points $\{(x^{(1)},y^{(1)}),...,(x^{(m)},y^{(m)}) \}$ so we can write $\hat f = \hat f_m$ to remember that.
On the other hand $Y$ is the prediction we are making on a new data point $(x^{(m+1)},y^{(m+1)})$ by using the model constructed on the $m$ data points above. So the Mean Squared Error can be written as
$$ E[\hat f_m(x^{(m+1)}) - y^{(m+1)}]^2$$
Expanding the equation from the previous section
$$E[\hat f_m Y]=E[\hat f_m (f+ \epsilon)]=E[\hat f_m f+\hat f_m \epsilon]=E[\hat f_m f]+E[\hat f_m \epsilon]$$
The last part of the equation can be viewed as
$$ E[\hat f_m(x^{(m+1)}) \cdot \epsilon^{(m+1)}] = 0$$
Since we make the following assumptions about the point $x^{(m+1)}$:
- It was not used when constructing $\hat f_m$
- It is independent of all other observations $\{(x^{(1)},y^{(1)}),...,(x^{(m)},y^{(m)}) \}$
- It is independent of $\epsilon^{(m+1)}$
Other sources with full derivations