Find the MSE of a true response and its predicted value using OLS estimation From Theodoridis' Machine Learning, exercise 3.26.

Consider, once more, the same regression as that of Problem 3.8, but with $\boldsymbol\Sigma_{\boldsymbol\eta} = \mathbf{I}_N$.

For context, this is the regression model
$$y_n = \boldsymbol\theta^{T}\mathbf{x}_n + \eta_n\text{, } \qquad n = 1, 2, \dots, N$$
where $\boldsymbol\eta \sim \mathcal{N}(\mathbf{0}, \boldsymbol\Sigma_{\boldsymbol\eta})$ and the $\mathbf{x}_n$ are considered fixed.

Compute the MSE of the predictions $\mathbb{E}[(y-\hat{y})^2]$, where $y$ is the true response and $\hat{y}$ is the predicted value, given a test point $\mathbf{x}$ and using the LS estimator
$$\hat{\boldsymbol\theta}=(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}\text{.}$$
The LS estimator has been obtained via a set of $N$ measurements, collected in the (fixed) input matrix $\mathbf{X}$ and $\mathbf{y}$ [....] The expectation $\mathbb{E}[\cdot]$ is taken with respect to $y$, the training data $\mathcal{D}$, and the test points $\mathbf{x}$. Observe the dependence of the MSE on the dimensionality of the space.
Hint: Consider, first, the MSE, given the value of a test point $\mathbf{x}$, and then take the average over all the test points.

My attempt:
Theodoridis shows on p. 81 that the generalization error $\mathbb{E}_{y\mid\mathbf{x}}\mathbb{E}_{\mathcal{D}}\left[\left(y-f(\mathbf{x};\mathcal{D})\right)^2\right]$ at $\mathbf{x}$ is
$$\mathrm{MSE}(\mathbf{x}) = \sigma^2_{\eta}+\mathbb{E}_{\mathcal{D}}\left[\left(f(\mathbf{x};\mathcal{D})-\mathbb{E}_{\mathcal{D}}f(\mathbf{x};\mathcal{D})\right)^2\right]+\left(\mathbb{E}_{\mathcal{D}}f(\mathbf{x};\mathcal{D})-\mathbb{E}[y\mid\mathbf{x}]\right)^2$$
Setting $$\hat{y}=f(\mathbf{x};\mathcal{D})=\mathbf{x}^{T}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}$$
we obtain
$$\mathbb{E}_{\mathcal{D}}\hat{y}=\mathbf{x}^{T}\mathbb{E}_{\mathcal{D}}\left[(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}\right]=\mathbf{x}^{T}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{X}\boldsymbol\theta=\mathbf{x}^{T}\boldsymbol\theta$$
so $$\left(f(\mathbf{x};\mathcal{D})-\mathbb{E}_{\mathcal{D}}f(\mathbf{x};\mathcal{D})\right)^2 = \left\{\mathbf{x}^{T}\left[(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{X}\mathbf{y}-\boldsymbol\theta\right]\right\}^2$$
This looks disgusting (and not easily simplified), so I'm guessing I'm doing something wrong. We also have $\sigma^2_\eta = 1$ by assumption.
Edit: This appears to be solved in the 12th printing (Jan 2017) of Elements of Statistical Learning by Hastie et al. (see here), equation (2.47) on p. 37, but they seem to have skipped showing the details (not to mention, I find the notation confusing).
 A: Here are two facts that we can utilize: if $A$ is a constant matrix, then $$\operatorname{E}[AX]=A\operatorname{E}[X],\qquad\operatorname{Var}(AX)=A\operatorname{Var}(X)A'.$$
For convenience let's denote $\hat{y}=f(\mathbf{x};\mathcal{D})$. Notice that $\operatorname{E}_{\mathcal{D}}\left[\left(\hat{y}-\operatorname{E}[\hat{y}]\right)^2\right]$ is just a fancy way of writing $\operatorname{Var}_{\mathcal{D}}(\hat{y})=\operatorname{Var}_{\mathcal{D}}(x'(X'X)^{-1}X'y)$. Thus, $$\operatorname{E}_{\mathcal{D}}\left[\left(\hat{y}-\operatorname{E}[\hat{y}]\right)^2\right]=x'(X'X)^{-1}X'\operatorname{Var}_{\mathcal{D}}(y)X(X'X)^{-1}x=x'(X'X)^{-1}x.$$ On the other hand, you already know that $\operatorname{E}_{\mathcal{D}}[\hat{y}]=x'\theta$, since $\operatorname{E}[\eta]=0$, the bias part $$\operatorname{E}_{\mathcal{D}}[\hat{y}]-\operatorname{E}[y\mid x]=0$$ as well. Therefore, the MSE over one test point $x$ is $$\operatorname{MSE}(x)=1+x'(X'X)^{-1}x.$$ 
To average over the test set, note that $$x'(X'X)^{-1}x=\operatorname{tr}(x'(X'X)^{-1}x)=\operatorname{tr}((X'X)^{-1}xx').$$ So $$\operatorname{E}_x[x'(X'X)^{-1}x]=\operatorname{tr}((X'X)^{-1}\operatorname{E}_x(xx')).$$
