What is the variance of residual of regression (matrix form)?

Question

I am confused at the variance of residual in regression (matrix form)

Given $Y = X*\beta \hat{}+ \epsilon$

Y is a nx1 predicted vector

X is a nxp matrix

$\beta\hat{}$ is a px1 coefficient vector = $(X^TX)^{-1}X^TY$

$\epsilon$ is a nx1 residual vector

Then $var(\epsilon)=var(Y-X\beta\hat{}$) = var($Y$ - $X(X^TX)^{-1}X^TY$) (1)

Let projection matrix $P$ = $X(X^TX)^{-1}X^T$ . It is idempotent $(P^n=P)$

(1) becomes $var(Y - PY)$ (2)

Method 1: (2) = $var ((I-P)Y) = (I-P)Var(Y)(I-P)^T = \sigma^2(I-P)$

Method 2: (2) = $var(Y)+var(PY) = \sigma^2I + P\sigma^2P^T = \sigma^2(I+P)$

I expect that both methods should result in the same formula but they did not. I am not sure what is missed here

Answer 1

The second formula you have assumed independence (or at least zero covariance) of $Y$ and $PY$, that does not hold.

See here, you forgot two important terms.

$$var (Y - PY) = var (Y) + var(PY) + cov(Y, -PY) + cov(-PY, Y) $$