Variance Estimate in linear regression In a linear regression, $y=X\beta+\epsilon$, where $\epsilon\sim N(0, \sigma^2)$, $X\sim R^{N \times (p+1)}$. Assume the observations $y_i$ are uncorrelated and have constant variance $\sigma^2$, and that the $x_i$ are fixed. Then $\hat{\beta} = (X^T X)^{-1} X^T y$.
One estimate the variance $\sigma^2$ by
$\hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N (y_i-\hat{y}_i)^2$. 
How to prove $E(\hat{\sigma}^2) = \sigma^2$?
and why
$\hat{\beta}\sim N(\beta, (X^T X)^{-1}\sigma^2)$ ? I know how to get the mean and variance of $\hat{\beta}$, but why it follows a normal distribution?
 A: For a random vector $W \in \mathbb R^k$ with $\operatorname{E}(W) = \mu$ we have $\operatorname{var}(W) = \operatorname{E}((W-\mu)^T(W-\mu))$, and this variance is a $k\times k$ matrix.  If $A$ is a constant (i.e. not random) $\ell\times k$ matrix, then $AW$ is an $\ell\times 1$ random vector and its variance is an $\ell\times\ell$ matrix, and we have
$$
\operatorname{var}(AW) = A\Big( \operatorname{var}(W)\Big) A^T. \tag 1
$$
You have $\displaystyle y\sim \mathcal{N}_N( X\beta, \sigma^2 I_N)$ where $I_N$ is the $N\times N$ identity matrix. (Recall that $X$ is an $N\times (p+1)$ matrix and $\beta$ is a $(p+1)\times 1$ matrix, so $X\beta$ is an $N\times1$ matrix.)  So we have
\begin{align}
\operatorname{E}( \hat\beta ) & = \operatorname{E}((X^TX)^{-1}X^T y) \\[4pt]
& = (X^TX)^{-1} X^T \operatorname{E}(y) \\[4pt]
& = (X^T X)^{-1} X^T X\beta \\[4pt]
& = \beta \text{ because two matrices that are each other's inverses cancel.} \\[12pt]
\operatorname{var}(\hat\beta) & = \Big( (X^T X)^{-1} X^T\Big) \Big( \operatorname{var}(y)\Big)\Big( X(X^T X)^{-1}\Big) \\[6pt]
& = \Big( (X^T X)^{-1} X^T\Big) \Big( \sigma^2 I_N \Big)\Big( X(X^T X)^{-1}\Big) \\[8pt]
& = \sigma^2 \Big( (X^T X)^{-1} X^T \Big)\Big( X (X^T X)^{-1} \Big) \\[6pt]
& = \sigma^2 (X^T X)^{-1} (X^T X) (X^T X)^{-1} \\[6pt]
& = \sigma^2 (X^T X)^{-1}.
\end{align}
Why is it normally distributed?  Essentially it's because each of its components is a linear combination of $y_1,\ldots,y_n$, which are independent and normally distributed.  Recall that to say that a vector (such as $\hat\beta$) is normally distributed means that every constant (i.e. non-random) linear combination of its componenets has a one-dimensional normal distribution.  And every linear combination of the componenets of $\hat\beta$ is a linear combination of a bunch of linear combinations of $y_1,\ldots,y_n$, and thus is a linear combination of $y_1,\ldots,y_n$.
The matrix $X$ is $N\times(p+1)$ and has rank $p+1$, so the "hat matrix" $H = X(X^T X)^{-1}X^T$ is $N\times N$ and has rank $p+1$.
An $N\times 1$ vector $w$ is in the column space of $X$ if and only if for some $(p+1)\times1$ vector $\gamma$ we have $w = X\gamma$.  Then
$$ Hw = X(X^T X)^{-1} X^T w = X(X^T X)^{-1} X^T X\gamma = X\gamma = w, $$
but if an $N\times 1$ vector $w$ is orthogonal to the column space of $X$, then $Hw=0$ because $X^Tw=0$.  Therefore:

$Hw$ is always the orthogonal projection of $w$ onto the column space of $X$.

Now recall that $\hat y = Hy$ (indeed, that is why $H$ is called the "hat matrix"), and $\operatorname{E}(Hy) = X\beta = \operatorname{E}(y),$ so $\operatorname{E}(y - \hat y)=0$.  And by using $(1)$ above and the fact that $y-\hat y= (I_N-H)y$, you can show that $\operatorname{var}(y-\hat y) = \sigma^2 (I_N-H)$.
Now consider an orthonormal basis of $\mathbb R^N$ that consists of $p+1$ vectors in the column space of $X$ and $N-p-1$ vectors orthogonal to the column space of $X$.  Since $w\mapsto Hw$ is the orthogonal projection onto the column space of $X$, the matrix of $I_N - H$ with respect to this new basis is
$$
\begin{bmatrix} 0 \\ & 0 \\ & & 0 \\ & & & \ddots \\ & & & & 0 \\ & & & & & 1 \\ & & & & & & 1 \\ & & & & & & & \ddots \\ & & & & & & & & 1 \end{bmatrix}
$$
where the diagonal has $p+1$ entries equal to $0$ and the rest, $N-p-1$ of them equal to $1$.  Let the components of $y_1,\ldots,y_N$ in the new coordinate system by $u_1,\ldots,u_N$; the the components of $y-\hat y$ in the new coordinate system must be $\underbrace{0,\ldots,0}_{p+1},\,\underbrace{ u_{p+2}, u_{p+3}, \ldots, u_N}_{N-p-1}$.  And those last $N-p-1$ components have expected value $0$.  Let $u$ be the vector whose components are $u_1,\ldots,u_N$; then $u$ is normally distributed with an expected value whose last $N-p-1$ components are $0$ and whose variance is $\sigma^2 I_N$.  Thus
$$
u_{p+2}, u_{p+3}, u_{p+4}, \ldots, u_N \sim \text{i.i.d. }\mathcal N(0,\sigma^2),
$$
and hence
$$
\underbrace{ (y_1-\hat y_1)^2 + \cdots + (y_N - \hat y_N)^2}_{N \text{ terms}} = \underbrace{u_{p+2}^2 + \cdots + u_N^2}_{N-p-1 \text{ terms}} \sim \chi^2_{N-p-1}
$$
and that has expected value $N-p-1$. That is why $\operatorname{E}(\hat\sigma^2) = \sigma^2$.
