T_M's answer addresses the first part of the question, namely, how (3.6) implies (3.8). I will address the second part, why use

$ \hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2. $

(This was probably answered many times elsewhere, but it's easier to repeat it in the convenient notation than to translate other answers.)

If we make additional modelling assumptions (see below), then $ \hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2$ is an unbiased estimator of the true variance. This justifies using this expression even when no such assumptions are made, because it is better than nothing. The sharp statement is as follows.

Claim. Let the data be generated by the true model $Y = \sum_{j=0}^m \beta_j X_j + \mathcal{E}$ with $E[\mathcal{E}]=0$, $\mathrm{var}[\mathcal{E}]=\sigma^2$, and uncorrelated errors (meaning that $\mathrm{var}[\mathbf{y}] = \sigma^2\mathbf{I}$). Then the least-squares estimate $\hat{\mathbf{y}}:=\mathbf{X}\hat{\boldsymbol{\beta}}$ satisfies $E[\sum_{i=1}^N(Y_i-\hat{Y}_i)^2]=(N-p-1)\sigma^2$, where $p+1$ is the number of linearly independent columns in $\mathbf{X}$.

For simplicity, let all $m+1$ columns of $\mathbf{X}$ be linearly independent, so $m=p$. Each row of $N\times(p+1)$ matrix $\mathbf{X}$ has the form $[1, X_1, \dots, X_p]$. I will use $\mathbf{X}^+$ to denote the (Moore-Penrose) pseudo-inverse of $\mathbf{X}$. If you are not comfortable with it, just remember that for a matrix with independent columns we write $\mathbf{X}^+ = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$.

Proof. The least-squares solution for parameters is given by $\hat{\boldsymbol{\beta}}=\mathbf{X}^+\mathbf{y}$, so $\hat{\mathbf{y}}=\mathbf{X}\mathbf{X}^+\mathbf{y}$. We can write $$ \begin{align} \sum_{i=1}^N(Y_i-\hat{Y}_i)^2 & = (\mathbf{y} - \hat{\mathbf{y}})^T (\mathbf{y} - \hat{\mathbf{y}}) = \mathbf{y}^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} \\ & = \mathrm{tr}\, \mathbf{y}^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} = \mathrm{tr}\, (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} \mathbf{y}^T. \end{align} $$ In the last transition we used the property of trace $\mathrm{tr}\,AB = \mathrm{tr}\,BA$ and that $P:=\mathbf{I} - \mathbf{X}\mathbf{X}^+$ is an orthogonal projection: orthogonal projections satisfy $P = P^T = P^2$. We will also require another property: that $P$ is an orthogonal projection onto the orthogonal complement of the column-space (range) of $\mathbf{X}$. That is, for any linear combination of columns, $\mathbf{X}\mathbf{w} = \sum_{j=0}^{p} \mathbf{x}_p w_p$, we have $P \mathbf{X}\mathbf{w} = \mathbf{0}$. Also, $\mathrm{tr}\,P = N-p-1$. All these properties can be demonstrated without a reference to the pseudo-inverse and orthogonal projections, but they are best understood at this level of abstraction.

The observed data is $\mathbf{y} = \mathbf{f}(\mathbf{X}) + \mathbf{e}$, where $\mathbf{f}(\mathbf{X})\equiv \mathbf{f} := E[\mathbf{y}]$. From the assumptions about the errors and $X$'s being "fixed", we have $E[\mathbf{y}\mathbf{y}^T] = \mathbf{f}\mathbf{f}^T + \sigma^2\mathbf{I}$. We use this in the expectation of the sum of squared residuals:

$$ \begin{align} E[\sum_{i=1}^N(Y_i-\hat{Y}_i)^2] & = E[ \mathrm{tr}\, (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} \mathbf{y}^T] = \mathrm{tr}\, (\mathbf{I} - \mathbf{X}\mathbf{X}^+) E[ \mathbf{y} \mathbf{y}^T] \\ & = \mathrm{tr}\, \{ (\mathbf{I} - \mathbf{X}\mathbf{X}^+) \mathbf{f}\mathbf{f}^T \} + \sigma^2 \mathrm{tr}\, \{ (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\} \\ & = \mathrm{tr}\, \{ \underbrace{(\mathbf{I} - \mathbf{X}\mathbf{X}^+) \mathbf{f}}_{=\mathbf{0}}\mathbf{f}^T \} + (N-p-1) \sigma^2. \end{align},$$ where we again used the properties of the projection. $\blacksquare$.

T_M's answer addresses the first part of the question, namely, how (3.6) implies (3.8). I will address the second part, why use

$ \hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2. $

(This was probably answered many times elsewhere, but it's easier to repeat it in the convenient notation than to translate other answers.)

If we make additional modelling assumptions (see below), then $ \hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2$ is an unbiased estimator of the true variance. This justifies using this expression even when no such assumptions are made, because it is better than nothing. The sharp statement is as follows.

Claim. Let the data be generated by the true model $Y = \sum_{j=0}^m \beta_j X_j + \mathcal{E}$ with $E[\mathcal{E}]=0$, $\mathrm{var}[\mathcal{E}]=\sigma^2$, and uncorrelated errors (meaning that $\mathrm{var}[\mathbf{y}] = \sigma^2\mathbf{I}$). Then the least-squares estimate $\hat{\mathbf{y}}:=\mathbf{X}\hat{\boldsymbol{\beta}}$ satisfies $E[\sum_{i=1}^N(Y_i-\hat{Y}_i)^2]=(N-p-1)\sigma^2$, where $p+1$ is the number of linearly independent columns in $\mathbf{X}$.

For simplicity, let all $m+1$ columns of $\mathbf{X}$ be linearly independent, so $m=p$. Each row of $N\times(p+1)$ matrix $\mathbf{X}$ has the form $[1, X_1, \dots, X_p]$. I will use $\mathbf{X}^+$ to denote the (Moore-Penrose) pseudo-inverse of $\mathbf{X}$. If you are not comfortable with it, just remember that for a matrix with independent columns we write $\mathbf{X}^+ = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$.

Proof. The least-squares solution for parameters is given by $\hat{\boldsymbol{\beta}}=\mathbf{X}^+\mathbf{y}$, so $\hat{\mathbf{y}}=\mathbf{X}\mathbf{X}^+\mathbf{y}$. We can write $$ \begin{align} \sum_{i=1}^N(Y_i-\hat{Y}_i)^2 & = (\mathbf{y} - \hat{\mathbf{y}})^T (\mathbf{y} - \hat{\mathbf{y}}) = \mathbf{y}^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} \\ & = \mathrm{tr}\, \mathbf{y}^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)^T (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} = \mathrm{tr}\, (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} \mathbf{y}^T. \end{align} $$ In the last transition we used the property of trace $\mathrm{tr}\,AB = \mathrm{tr}\,BA$ and that $P:=\mathbf{I} - \mathbf{X}\mathbf{X}^+$ is an orthogonal projection: orthogonal projections satisfy $P = P^T = P^2$. We will also require another property: that $P$ is an orthogonal projection onto the orthogonal complement of the column-space (range) of $\mathbf{X}$. That is, for any linear combination of columns, $\mathbf{X}\mathbf{w} = \sum_{j=0}^{p} \mathbf{x}_p w_p$, we have $P \mathbf{X}\mathbf{w} = \mathbf{0}$. Also, $\mathrm{tr}\,P = N-p-1$. All these properties can be demonstrated without a reference to the pseudo-inverse and orthogonal projections, but they are best understood at this level of abstraction.

The observed data is $\mathbf{y} = \mathbf{f}(\mathbf{X}) + \mathbf{e}$, where $\mathbf{f}(\mathbf{X})\equiv \mathbf{f} := E[\mathbf{y}]$. From the assumptions about the errors and $X$'s being "fixed", we have $E[\mathbf{y}\mathbf{y}^T] = \mathbf{f}\mathbf{f}^T + \sigma^2\mathbf{I}$. We use this in the expectation of the sum of squared residuals:

$$ \begin{align} E[\sum_{i=1}^N(Y_i-\hat{Y}_i)^2] & = E[ \mathrm{tr}\, (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\mathbf{y} \mathbf{y}^T] = \mathrm{tr}\, (\mathbf{I} - \mathbf{X}\mathbf{X}^+) E[ \mathbf{y} \mathbf{y}^T] \\ & = \mathrm{tr}\, \{ (\mathbf{I} - \mathbf{X}\mathbf{X}^+) \mathbf{f}\mathbf{f}^T \} + \sigma^2 \mathrm{tr}\, \{ (\mathbf{I} - \mathbf{X}\mathbf{X}^+)\} \\ & = \mathrm{tr}\, \{ \underbrace{(\mathbf{I} - \mathbf{X}\mathbf{X}^+) \mathbf{f}}_{=\mathbf{0}}\mathbf{f}^T \} + (N-p-1) \sigma^2. \end{align},$$ where we again used the properties of the projection. $\blacksquare$.

Note that:

• the assumption of uncorrelated homoskedastic errors ($\mathrm{var}\,\mathbf{y} = \sigma^2 \mathbf{I}$) is weaker than then assumption of i.i.d. errors;
• there is no assumption of anything being Gaussian or normally distributed.