How do I derive prediction intervals for a general linear model?
My general linear model written in matrix form is,
$$ \mathbf{Y} = \mathbf{X} \mathbf{B} + \mathbf{R}$$
with each of the rows of the residual matrix $\mathbf{R}$, distributed as $\mathbf{r}_i \sim \mathcal{N}(\mathbf{0},\Sigma_{R})$. Alternatively, $row(\mathbf{R}) \sim \mathcal{N}(\mathbf{0}, \mathbf{I} \otimes \Sigma_{R})$.
I have a least squares estimator as
$$\hat{\mathbf{B}} = (\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T} \,, $$
from affine transformations on $\mathbf{R}$, this estimator is distributed as,
$$ row(\hat{\mathbf{B}}) \sim \mathcal{N}(\mathbf{B}, \Sigma_{\hat{\mathbf{B}}}) $$
What I wish to do is form a prediction interval on a new instance of $\mathbf{Y}$ that is within the range (an interpolation, not extrapolation). This new instance is $\tilde{\mathbf{y}}_{i}$. What I wish to do is find a prediction interval (or confidence interval?) on this new instance.
On my own, I have used affine transformations again on $\hat{\mathbf{B}}$ to find that
$$\tilde{\mathbf{y}}_{i} \sim \mathcal{N}(\tilde{\mathbf{x}} \hat{\mathbf{B}}, \mathbf{X}^{T} (\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X} \otimes \hat{\Sigma}_{R} ) $$
However I am uncertain in my solution. My textbooks on regular multiple linear regression have formulas for prediction intervals. But if the LS estimator has a distribution, why not just use this result? From it I can make $1\sigma$ or $3 \sigma$ intervals etc. In short,
$\textbf{How does one derive prediction intervals for the general linear model? }$