When using GLS, how do I compute the conditional variance of the error term? I want to use generalized least squares (GLS) to estimate $\beta$ in
$$r = X\beta + u$$
where $r$ is a vector of stock returns and $X$ is a matrix of factor exposures.  I choose GLS because the relationship
$$\hat{\beta} = (X^T\Omega^{-1}X)^{-1}X^T\Omega^{-1}r,$$
where $\Omega = Var(u | X)$ is the conditional variance of the error term, will allow me to express the factors in terms of the original vector $r$:
$$\beta_k = \sum_{n = 1}^{k} c_{k,n}  r_n.$$
My question is:  When using GLS, how can I compute $\Omega$?
EDIT: This was previously a question on how to do it in R.  I've changed it to a more theoretical question:  I'll figure out the R way after I understand how to get $\Omega$.  
 A: From my perspective, there are two possibilities for you.
Solution from classical Econometrics
Either you choose the path of standard econometrics (i.e. the Wooldridge solution) and you try to approximate the weights within $\Omega$ by first applying a linear regression $\hat{\beta}=(X'X)^{-1}X'y$ and using the residuals $\hat{u}=y-X\hat{\beta}$ in order to approximate the weights in via $\hat{\Omega}=\text{diag}\{\hat{u}_i^2\}_{i=1,...,N}$. This approach is called weighted least-squares and is helpful in presence of heteroscedasticity. 
Another classical idea from econometrics is using the residuals from the first stage (the $\hat{u}$'s) for a regression framework that allows  the residuals to depend on the covariates. I.e. using something like $\hat{u}^2=\exp(Z\gamma)$, where $Z$ is a subset of $X$. After that you can use the fitted values from this regression (call them $\hat{\omega}$) in order to define your weights $\hat{\Omega}=\text{diag}\{\hat{\omega}_i\}_{i=1,...,N}$. 
Solution from Statistics
If you want to stay within the framework of generalized linear models, I would recommend you the application of hierarchical generalized linear models (HGLM, see https://journal.r-project.org/archive/2010-2/RJournal_2010-2_Roennegaard~et~al.pdf for a nice implementation in R).
These are very flexible in modelling covariance structures and they do not rely on a two-step estimation procedure. The basic structure is given by
$$Y|\beta,\gamma,X \sim N (X\beta,\exp(Z\gamma)) .$$
With this approach it is possible to fit the dispersion term very flexible. And the nice thing of this approach is the possiblity to employ other distributions from the exponential family (i.e. you are not restricted to gaussian). On top of that, they also allow for a very convenient inclusion of random effects.
