Is the OLS estimator the UMVUE (assuming Normality)? Suppose
$$
\mathbf{y} = \mathbf{X} \mathbf{b} + \mathbf{e} \, ,
\\
\mathbf{e} \sim \mathcal{N}(0,\mathbf{I}_P) \, .
$$
We know that $\mathbf{\hat{b}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}$ is the BLUE.
Is it also the UMVUE? I can only find a single source (page 6) that claims this, so I'm unsure.
In case $\mathbf{X}=\mathbf{1}$ it is true ($\mathbf{\hat{b}}$ becomes the sample mean).
But other, related results like Stein's example make me cautious.
And if it's true, then why isn't it more famous?
 A: Under the assumptions
$$
\begin{align}
&\mathbf{y} = \mathbf{X} \mathbf{b} + \mathbf{e}, \;\mathbf X \;\text{full column rank},\\
&\mathbf e \mid \mathbf{X}  \sim \mathop{\mathcal{N}}\left(\mathbf 0,\sigma^2\mathbf{I}\right),\;\sigma^2 \in \mathbb{R}_{>0}
\end{align}
$$
the OLS estimator $\hat{\mathbf{b}}=\left(\mathbf X^{\mathsf T}\mathbf X\right)^{-1}\mathbf X^{\mathsf T}\mathbf y$ is the UMVUE of $\mathbf b$.
This is clear from the facts that $\hat{\mathbf b}$ is unbiased and that  $\mathop{\mathbb{V}}\left(\hat{\mathbf b}\right) = \sigma^2 \left(\mathbf X^{\mathsf T}\mathbf X\right)^{-1}$ is the inverse expected Fisher information of $\mathbf b$, i.e., $\hat{\mathbf b}$ attains the Cramér–Rao lower bound.
This result is in a sense more general than the Gauss–Markov theorem in that it's not restricted to linear estimators. On the other hand it's about linear regression with i.i.d. normal errors only.
Interestingly, under the more general assumptions
$$
\begin{align}
&\mathbf{y} = \mathbf{X} \mathbf{b} + \mathbf{e}, \;\mathbf X \;\text{full column rank},\\
&\mathbf e = \left(e_1,\ldots, e_n\right)^\mathsf{T},\\
&e_1,\ldots, e_n  \mid \mathbf X \overset{\text{(c.)i.i.d.}}{\sim} \left(0, \sigma^2\right),\;\sigma^2 \in \mathbb{R}_{>0},
\end{align}
$$
the OLS estimator $\hat{\mathbf{b}}=\left(\mathbf X^{\mathsf T}\mathbf X\right)^{-1}\mathbf X^{\mathsf T}\mathbf y$ is also the UMVUE of $\mathbf b$ if $\hat{\mathbf{b}}$ is unbiased for all regression models that satisfy
$$
\begin{align}
&\mathbf{y} = \mathbf{X} \mathbf{b} + \mathbf{e}, \;\mathbf X \;\text{full column rank},\\
&\mathop{\mathbb{E}}\left(\mathbf e\mid \mathbf{X}\right)=\mathbf 0,\\
&\mathbf e = \left(e_1,\ldots, e_n\right)^\mathsf{T},\\
&e_1,\ldots, e_n  \mid \mathbf X \;\text{(conditionally) independent},\\
&\mathop{\mathbb{V}}\left(\mathbf e \mid \mathbf{X}\right)=\mathop{\mathrm{diag}}\left(\sigma^2_1,\ldots,\sigma^2_n\right),\; \sigma^2_i \in \mathbb{R}_{>0},\\
\end{align}
$$
i.e., for all linear regression models with independent and homo- or heteroscedastic errors (in particular, not only for the data-generating class of linear regression models with independent and homoscedastic errors).

References

*

*Hansen, B. E. (2022). A modern Gauss–Markov theorem. Econometrica, 90(3), 1283–1294.

*Pötscher, B. M., & Preinerstorfer, D. (2022). A Modern Gauss-Markov Theorem? Really?. arXiv:2203.01425v3

