Variance of Y in OLS linear regression If in a linear regression of Y=XB+U in matrix form, where U is the vector of error terms, each of which is normally distributed with (0, sigma^2), and it satisfies all of those classical assumptions (E(U|X)=0, etc), then what is the variance of Y?
 A: First you have to differentiate between 'a priori' (assumption before regression) variance $\sigma_0^2$ and 'a posteriori' (computed value with regression) variance $\hat{\sigma}_0^2$. Your a priori variance is by your own definition $\sigma_0^2$ and is used to build the cofactor (weight) matrix for the $n$ uncorrected observations $l$ (in your notation it's $Y$).
$\Sigma_{ll} =
\frac{1}{\sigma _0^2}\begin{bmatrix}
\sigma _{l_1}^2 & \sigma _{l_1} \sigma _{l_2} \rho _{l_{1,2}} & ... & \sigma _{l_1} \sigma _{l_{n-1}} \rho _{l_{1,n-1}} & \sigma _{l_1} \sigma _{l_n} \rho _{l_{1,n}} \\
\sigma _{l_2} \sigma _{l_1} \rho _{l_{2,1}} & \sigma _{l_2}^2 & \text{} & \text{} & \sigma _{l_2} \sigma _{l_n} \rho _{l_{2,n}} \\
\vdots & \text{} & \text{} & \ddots & \vdots \\
\sigma _{l_{n-1}} \sigma _{l_1} \rho _{l_{n-1,1}} & \text{} & \text{} & \text{} & \sigma _{l_{n-1}} \sigma _{l_n} \rho _{l_{n-1,n}} \\
\sigma _{l_n} \sigma _{l_1} \rho _{l_{n,1}} & \sigma _{l_n} \sigma _{l_2} \rho _{l_{n,2}} & ... & \sigma _{l_n} \sigma _{l_{n-1}} \rho _{l_{n,n-1}} & \sigma _{l_n}^2
\end{bmatrix}$
Here $\sigma_{l_i}$ are variances and $\rho_{l_i,l_i}$ correlations for your observations. If your observations are not correlated, $\Sigma_{ll}$ is an diagonal matrix. If all your observations are equally precise then $\Sigma_{ll}$ is the unit matrix. As you see $\sigma_0^2$ can be used 'like a' normalization factor to achieve better numeric stability.
Once your regression is computed you can compute the residual vector $r$ and use them to compute the corrected $\hat{l}$ values: $\hat{l}=l+r$. If you use 'common' least squares, which looks something like this:
$\hat{x} = \underbrace{\left(A^T\Sigma_{ll}^{-1}A\right)^{-1}}_{\Sigma{\hat{x}\hat{x}}}A^{T}\Sigma_{ll}^{-1}l$
$r = A\hat{x}-l$
the a posteriori variance factor can be computed as 
$\hat{\sigma}_0=\sqrt{\frac{r^T\Sigma_{ll}^{-1}r}{n-u}}$($u=2$ since we have two unknowns).
Using this approach your $\Sigma_{\hat{x}\hat{x}}$ matrix contains the variances and covariances ($\Sigma_{ll}$ equivalent for parameters) for the computed parameter vector $\hat{x}$. You can use the $\Sigma_{\hat{x}\hat{x}}$ matrix to compute the variance not only for $\hat{l}$ but also for any function where $\hat{x}$ is used. The approach is called 'Propagation of uncertainty' and is computed as follows:
$\Sigma_{\hat{l}\hat{l}}=J\Sigma_{\hat{x}\hat{x}}J^{T}$
where the standard deviation for each observation $l_i$ is
$\sigma_{l_i} = \hat{\sigma}_0\sqrt{\Sigma_{\hat{l_i}\hat{l_i}}}$
where $J$ is the 'Jacobian matrix' of the desired function (derivative with respect to each observation and parameter). In your case $J$ is $A$ since it's by definition the Jacobian matrix for $l$.
