Variance Estimate in linear regression

Question

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?

Answer 1

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$.