The variance-covariance matrix of the least squares parameter estimation I'm learning Linear Regression for Regression from "The Elements of Statistical Learning".
Why

The variance-covariance matrix of the least squares parameter
estimates is easily derived from (3.6) and is given by
$$ Var(\hat{\beta}) = (X^TX)^{-1}\sigma^2. \tag{3.8} $$
Typically one estimates the variance $\sigma^2$ by
$$ \hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2. $$

Could someone explain the above formulas in detail, like the following:
$$
E(\hat{\beta})=E((X^TX)^{-1}X^Ty)   \\
=(X^TX)^{-1}X^TE(y)  \\
=(X^TX)^{-1}X^TX\beta  \\
=\beta
$$

As mentioned in the comments, Equation (3.6) is

$\hat\beta = (\mathbf{X}^T \mathbf{X})^{−1}\mathbf{X}^T \mathbf{y} \tag{3.6}$

 A: The other answer here and the answers on a later version of this question here [Covariance matrix of least squares estimator $\hat{\beta}$  are not correct
In the book you are referencing, the data $x_1,\dots,x_N$ ($x_i^{\top}$ is  the ith row of $\mathbf{X}$) are not random. The authors say that the $y_i$ are uncorrelated with constant variance. And we have the formula
$$
\hat{\beta} = (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{y}.
$$
That's really all they say. There is no assumption that the real distribution of $Y$ is a linear function of $X$ plus a noise. And there is no explicit assumption that $\mathbb{E}(Y|X) = 0$. So, if you try to work with the information you are actually given in the book, you'll do something like this:
First we compute the expectation:
$$
\mathbb{E}(\hat{\beta}) = (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y}) = (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y}) 
$$
So
\begin{align}
\mathbb{E}(\hat{\beta})\mathbb{E}(\hat{\beta})^T &= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y}) \Bigl((\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y})\Bigr)^{\top} \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y}) \mathbb{E}(\mathbf{y})^{\top} \mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1} \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y}) \mathbb{E}(\mathbf{y})^{\top} \mathbf{X}\bigl((\mathbf{X}^{\top}\mathbf{X})^{-1}\bigr)^{\top}  \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y}) \mathbb{E}(\mathbf{y})^{\top} \mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}
\end{align} 
And 
\begin{align}
\mathbb{E}(\hat{\beta}\hat{\beta}^T) &= \mathbb{E}\biggl((\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{y}\Bigl( (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{y}\Bigr)^{\top} \biggr)\\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}(\mathbf{y} \mathbf{y}^{\top}) \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} 
\end{align}
The variance-covariance matrix is the difference as usual, which comes out as
\begin{align}
&(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\bigl(\mathbb{E}(\mathbf{y} \mathbf{y}^{\top}) - \mathbb{E}(\mathbf{y}) \mathbb{E}(\mathbf{y})^{\top}  \bigr) \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1}  \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\bigl(\sigma^2 I_{N\times N}  \bigr) \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\sigma^2
\end{align}
So the only assumption that we had, I use explicitly at the end: We know the variance-covariance matrix of $\mathbf{y}$ is just $\sigma^2$ multiplied by the identity matrix.
A: Because $x_i$ are fixed, so
$$\mathrm{Var}[\hat{\beta}] = (\mathbf{X}^T\mathbf{X})^{-1}\mathrm{Var}[\mathbf{y}] $$
And
$$\mathrm{Var}[\mathbf{y}] = \mathrm{Cov}[\mathbf{y}]=
\left[\begin{array}{ccc}
\sigma_{11} & \cdots & \sigma_{1 n} \\
\vdots & \ddots & \vdots \\
\sigma_{p1} & \cdots & \sigma_{n n}
\end{array}\right]
=\sigma^2 \mathbf{I}
$$
Because $y_i$ are uncorrelated and have constant variance $\sigma^2$,
$$\sigma_{ij}=\mathrm{Cov}[y_i, y_j]=E[y_iy_j]-E[y_i]E[y_j]=\sigma^2\delta_{ij}$$
Therefore,
$$\mathrm{Var}[\hat{\beta}] = (\mathbf{X}^T\mathbf{X})^{-1} \sigma^2$$
A: 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.

