Why does the estimation of the variance of the outputs in linear regression include the fitted values? I was revisiting my notes about the classical linear regression model, $Y = X \beta$.
If we want to estimate the variance of the least squares estimator we usually suppose that the outputs $y_{i}$ are uncorrelated and have fixed variance $\sigma ^{2}$.
We then have $Var(\hat{\beta} ^{2}) = (X^{T}X)^{-1} \sigma ^{2}$ and we estimate $\sigma^{2}$ using $\hat{\sigma}^{2} = \frac{1}{N-p-1}\sum_{i=1}^{N}{y_{i}-\hat{y_{i}}}$. Where $N$ is the number of training examples and $p$ the number of predictors.
Why does the estimation of the variance of the outputs contains elements from the model ($p$ and $\hat{y}$)?
 A: $\DeclareMathOperator{\X}{\mathbf X^\mathsf T\mathbf X}\DeclareMathOperator{\ep}{\boldsymbol\varepsilon}$
Let $\mathbf M := \mathbf I-\mathbf X(\X)^{-1}\mathbf X^\mathsf T$ be the residual maker. Then
$$\mathbf e =\mathbf M\mathbf y. \tag 1$$
Now $\mathbf e^\mathsf T\mathbf e=\ep^\mathsf T\mathbf M\ep.$ Therefore its expectation
\begin{align}\mathbb E\left[\mathbf e^\mathsf T\mathbf e|\mathbf X\right]&= \mathbb E\left[\ep^\mathsf T\mathbf M\ep|\mathbf X\right]\\&= \mathbb E\left[\operatorname{tr}(\ep^\mathsf T\mathbf M\ep)|\mathbf X\right]\\ &= \mathbb E\left[\operatorname{tr}(\mathbf M\ep\ep^\mathsf T)|\mathbf X\right]\\ &= \operatorname{tr} \left(\mathbf M~\mathbb E\left[(\ep\ep^\mathsf T)|\mathbf X\right]\right)\\&= \operatorname{tr}(\mathbf M~\sigma^2\mathbf I)\\&=(n-P)\sigma^2\tag 2\label 2\end{align}
From $\eqref{2},$ one can deduce a natural (unbiased) estimator of $\sigma^2$ that is $(n-P)^{-1}(\mathbf e^\mathsf T\mathbf e). $

Reference:
$[\rm I]$ Econometric Analysis, William H. Greene, Pearson Education, $2018, $ chapter $4, $ p. $62.$
