What is the MSE of a model when the true data generating process is known? Say I fit an OLS model to get an estimate of the coefficients $\hat{\beta}$.
For some reason, I also happen to know with complete certainty:

*

*That the predictors come from a multivariate normal with mean at the origin and some covariance matrix

$$ \mathbf{x}_i \sim \operatorname{MVN}(\mathbf{0}, \Sigma)$$

*

*That the coefficients of the true data generating process are $\tilde{\beta}$.


*And that the error term $e_i \sim \mathcal{N}(0, \sigma)$ and are iid.
Because I know the true generating process, I should be able to compute the MSE of my model exactly.  When the model only has a single predictor which has standard normal distribution, the calculation is straight forward.  I'm wondering if someone can help guide me through the general case I've made here.
How can I calculate the MSE of my model exactly given I know the data generating process?
 A: I took this question to be asking about the MSE of the regression estimator, rather than the MSE of a prediction.  The comments by the OP have subsequently clarified that he wanted the latter, so this answer is looking at something different to what he wants.  I am going to leave this answer here anyway, since it is useful for analysis of the former.

In a regression context, usually the explanatory variables would be treated as fixed, so the MSE calculation would condition on them.  I will give you the conditional result first, followed by looking at the unconditional result.  Under OLS estimation you have $\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} = (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{\varepsilon}$, which means that the squared-norm of the estimation error can be written as the quadratic form:
$$\begin{align}
||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}||^2
&= (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\text{T} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \\[6pt]
&= ((\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{\varepsilon})^\text{T} (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{\varepsilon} \\[6pt]
&= \boldsymbol{\varepsilon}^\text{T} \mathbf{x} (\mathbf{x}^\text{T} \mathbf{x})^{-2} \mathbf{x}^\text{T} \boldsymbol{\varepsilon} \\[6pt]
&= \boldsymbol{\varepsilon}^\text{T} \mathbf{A}(\mathbf{x}) \boldsymbol{\varepsilon}, \\[6pt]
\end{align}$$
where $\mathbf{A}(\mathbf{x}) \equiv \mathbf{x} (\mathbf{x}^\text{T} \mathbf{x})^{-2} \mathbf{x}^\text{T}$ is an $n \times n$ matrix that is fully determined by the explanatory variables.  Using a well-known rule for the expected value of a quadratic form, the MSE conditional on the explanatory variables is:
$$\begin{align}
\text{MSE}(\hat{\boldsymbol{\beta}}, \boldsymbol{\beta} | \mathbf{x})
&\equiv \mathbb{E}(||\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}||^2 | \mathbf{x}) \\[6pt]
&= \mathbb{E}(\boldsymbol{\varepsilon}^\text{T} \mathbf{A}(\mathbf{x}) \boldsymbol{\varepsilon}) \\[6pt]
&= \text{tr}(\sigma^2 \mathbf{A}(\mathbf{x}) \mathbf{I}) \\[6pt]
&= \sigma^2 \text{tr}(\mathbf{A}(\mathbf{x})). \\[6pt]
\end{align}$$
Consequently, given any distribution for the (random) design matrix $\mathbf{X}$ we get the marginal MSE:
$$\begin{align}
\text{MSE}(\hat{\boldsymbol{\beta}}, \boldsymbol{\beta})
&= \sigma^2 \mathbb{E}( \text{tr} (\mathbf{A}(\mathbf{X}))) \\[6pt]
&= \sigma^2 \int \text{tr}(\mathbf{A}(\mathbf{x})) \cdot p (\mathbf{x}) \ d \mathbf{x}. \\[6pt]
\end{align}$$
Given your stipulated distribution for the design matrix of your regression, it should be possible to compute the expected value above, giving you the marginal MSE for the analysis.  I am not aware of a closed form solution in this case, so I would suggest simulating the expectation of interest through importance sampling or some MCMC method.
A: As mentioned in the comments, the question is about the expected predictive performance of a given, fitted model on new test data drawn from the data generating process (DGP). In this context, the mean squared error (MSE) is the expected squared prediction error, where the expectation is taken over the distribution of the test data, and the  fitted parameters are held fixed.
According to the DGP defined in the question, let $x \sim \mathcal{N}(\vec{0}, \Sigma)$ be a new test point with corresponding response $y \sim \mathcal{N}(x \beta + \beta_0, \sigma^2)$. The predicted response is $\hat{y} = x \hat{\beta} + \hat{\beta}_0$. Under these conditions, the MSE is:
$$MSE(y, \hat{y}) = \sigma^2
+ (\beta-\hat{\beta})^T \Sigma (\beta-\hat{\beta})
+ (\beta_0-\hat{\beta}_0)^2$$
This holds regardless of how the parameters were fit (e.g. using OLS or something entirely different). The derivation is below.
MSE for a fixed test point
The MSE for a given/fixed test point $x$ is the expected squared prediction error, where the expectation is taken w.r.t. the conditional distribution of the  response, given the test point:
$$MSE(y, \hat{y} \mid x) \ = \
E_{y|x} \big[ (y - \hat{y})^2 \big]$$
$$= E_{y|x} \big[ y^2 - 2 y \hat{y} + \hat{y}^2 \big]$$
By linearity of expectation:
$$= E_{y|x}[y^2] - 2 \hat{y} E_{y|x}[y] + \hat{y}^2$$
Recall that $\operatorname{Var}(y|x) = E_{y|x}[y^2] - E_{y|x}[y]^2$:
$$MSE(y, \hat{y} \mid x) \ = \
\operatorname{Var}(y|x)
+ E_{y|x}[y]^2
- 2 \hat{y} E_{y|x}[y]
+ \hat{y}^2$$
Factorize:
$$= \operatorname{Var}(y|x)
+ \Big( \hat{y} - E_{y|x}[y] \Big)^2$$
According to the DGP, the conditional mean of of the response is $E_{y|x}[y] = \beta^T x + \beta_0$ and the conditional variance is $\operatorname{Var}(y|x) = \sigma^2$. Substitute these in, along with the the predicted response $\hat{y} = \hat{\beta}^T x + \hat{\beta}_0$:
$$MSE(y, \hat{y} \mid x) \ = \
\sigma^2 + \Big(
  (\hat{\beta}^T x + \hat{\beta}_0)
  - (\beta^T x + \beta_0)
\Big)^2$$
Marginal MSE
The overall (marginal) MSE is obtained by averaging the MSE for a fixed test point over all possible test points, weighted by the probability density of each:
$$MSE(y, \hat{y}) = E_x \Big[ MSE(y, \hat{y} \mid x) \Big]$$
Plug in our expression for $MSE(y, \hat{y} \mid x)$ and crank through some algebra:
$$= E_x \Big[
  \sigma^2
  + (\hat{\beta} - \beta)^T x x^T (\hat{\beta} - \beta)
  + 2 (\hat{\beta}_0 - \beta_0) (\hat{\beta} - \beta)^T x
  + (\hat{\beta}_0 - \beta_0)^2
\Big]$$
By linearity of expectationn:
$$= \sigma^2
+ (\beta - \hat{\beta})^T E_x[x x^T] (\beta - \hat{\beta})
+ 2 (\beta_0 - \hat{\beta}_0) (\beta - \hat{\beta})^T E_x[x]
+ (\beta_0 - \hat{\beta}_0)^2$$
According to the DGP, the distribution over test points has mean zero: $E_x[x] = \vec{0}$. Thus, $E_x[x x^T] = \Sigma$ is the covariance matrix of the predictors and the marginal MSE is:
$$MSE(y, \hat{y}) \ = \ \sigma^2
+ (\beta - \hat{\beta})^T \Sigma (\beta - \hat{\beta})
+ (\beta_0 - \hat{\beta}_0)^2$$
