This is in scope of the simple linear model. Im trying to prove that $\mathbb{E}\left(\widehat{\sigma}^2\right) = \sigma^2$ for
$$\widehat{\sigma}^2 = \frac{1}{n-2}\sum^n_{i=1} \left(y_i-\widehat{y}_i\right)^2$$
where $$Y_i\sim N(\beta_0+\beta_1x_i,\sigma^2)$$ and $\widehat{y_i},i=1,2,3,...,n$ are pedricted values and $y_1,y_2,...,y_n$ is a sample from $Y_i$
$\newcommand{\b}{\begin{bmatrix}}\newcommand{\eb}{\end{bmatrix}}$ The vector of fitted values $$ \hat Y = \b \hat y_1 \\ \vdots \\ \hat y_n\eb $$ is the orthogonal projection of $$ Y = \b y_1 \\ \vdots \\ y_n \eb $$ onto the column space of the design matrix $$ X = \b 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n \eb. $$ The vector $\hat\varepsilon$ is $Y-\hat Y$. That means $\hat\varepsilon$ is the projection of $Y$ onto the $(n-2)$-dimensional orthogonal complement of that column space. Notice that this latter orthogonal projection maps the expected value $$ \mathbb E Y = \b 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n \eb \b \beta_0 \\ \beta_1 \eb $$ to the zero vector. Let $H$ (conventionally called the "hat matrix") be the matrix of the first orthogonal projection above, so that $\hat Y= HY$. Then $\varepsilon = (I-H)Y$. Recall that the $n\times n$ matrix $H$ of rank $2$ is $X(X^TX)^{-1}X^T$ and this is symmetric. It is also idempotent, i.e. $H^2=H$, or in other words, if you project onto a space, and then project that projection onto that same space, then you just get the point you had when you first projected onto the space. Similarly $I-H$ is symmetric and idempotent. So $$ \varepsilon \sim N_n(0, (I-H)\Big(\sigma^2 I\Big)(I-H)^T) = N_n(0, \sigma^2(I-H)). $$
You have a normally distributed random vector in an $(n-2)$-dimensional Euclidean space. Its expected value is $0$ and the normal distribution is spherically symmetric. Hence the square of its norm is distributed as $\sigma^2\chi^2_{n-2}$, so the expected value of the square of the norm is $\sigma^2(n-2)$.
Indeed, if you project $Y$ onto that $(n-2)$-dimensional space, getting $\hat\varepsilon$, and then let $U_1,\ldots, U_{n-2}$ be the coordinates of $\hat\varepsilon$ with respect to an orthonormal basis of that space then you have $$ U_1,\ldots,U_{n-2}\sim\mathrm{i.i.d.}\ N(0,\sigma^2) $$ and $$\varepsilon_1^2+\cdots+\varepsilon_n^2=U_1^2+\cdots+U_{n-2}^2.$$