Distribution of $\hat{y}=Hy$ I wish to find the distribution of $\hat{y}=Hy$ where $H$ is the hat matrix $X(X'X)^{-1}X'$ in which a dash represents the transpose. Also, $\epsilon$ is  $N(0, \sigma^2)$ distributed.
Thanks
 A: You have the assumption that the errors are multivariate normal $\epsilon \sim \mathcal MVN(0,\Sigma)$, where the errors is the complete vector $\epsilon = (\epsilon_1,...,\epsilon_N)$ assuming that there are $N$ observations. In some cases observations will be independent hence the covariance matrix does not need to allow for serial correlation hence all off-diagonal terms will be 0. In this case $\Sigma$ is then a diagonal matrix with the terms in the diagonal potentially varying allowing for heteroscedasticity. Under what is sometimes referred to as classical errors $\Sigma = \sigma^2 I_N$ all diagonal terms are the same and the errors are therefore all independent and identically distributed normal $\mathcal N(0,\sigma^2)$.
Nevertheless in any case if $\epsilon$ is $\mathcal MVN(\mu,\Sigma)$ then $A\epsilon$ is 
$\mathcal MVN(A\mu,A\Sigma A^\top)$ where $A$ is a constant matrix. Using this property consider linear regression model $y = X\beta + \epsilon$ in matrix version and premultiply with the hat matrix $H$ to get
$$\hat y = Hy = HX\beta + H\epsilon = X\beta + H\epsilon$$
Considering $X$ to be a constant matrix $X\beta$ is a constant vector simply changing the mean of the distribution of $H\epsilon$. The distribution of $H\epsilon$ is found simply by applying the above stated rule saying that because $\epsilon$ is mutivariate normal the linear transform $H\epsilon$ will also be multivariate normal and the mean is $\mathbb E[X\epsilon] = X \mathbb E[\epsilon] =\mathbf 0$ adding $X\beta$ therefore implies that $\mathbb E[X\beta + H\epsilon]= X\beta$. 
The variance is $H Var(\epsilon) H^\top = H \Sigma H^\top$. In the case where $\Sigma = \sigma^2I_N$ the expression for the variance reduces significantly to
$$H \Sigma H^\top = H \sigma^2I_N H^\top = \sigma^2 HH^\top = \sigma^2H$$
using idempotency and symmetry of $H$. The result is therefore
$$\hat y = X\beta + H\epsilon \sim \mathcal {MVN}(X\beta,H \Sigma H^\top )$$
and under assumption of classical errors (independence and homoscedasticity) this simplifies to
$$\hat y = X\beta + H\epsilon \sim \mathcal {MVN}(X\beta,\sigma^2H)$$
where $I_N$ is $N \times N$ identity matrix.
A: Assumptions:


*

*We are talking about a linear regression model

*$X$ is not random


Then
$$y = X\beta  + \epsilon$$
$$\hat{y} = Hy = X(X'X)^{-1}X'X\beta  + H\epsilon$$
$$\implies \hat{y} \sim N(X\beta ,H\sigma^2)$$
