Groupwise heteroskedasticity I know how to derive the GLS estimator of beta (theoretical GLS), but there a slight change to the question and i am not quite sure how to go about it.
A researcher has reason to believe that the disturbances in the standard model are heteroscedastic, with a threefold block-diagonal structure, such that:
$$Y_i = X_i\beta + U_i, \;\;\ i=1,2,3$$
where $Y_i$ is a $T_i\times 1$ column vector $X_i$ is a $T_i\times K$ regressor matrix and $U_i$ isa $T_i\times 1$ column vector of disturbances. $\beta$ is a $K \times 1$ of population regression coefficients, assumed to be the same for each of the blocks.  
Assume that the only departure from the standard assumptions is that:
$$E(U_iU_i') = \sigma^2_iI_{T_i}$$   
where $I_{T_i}$  is the identity matrix of order $T_i$ .
My question is simply how to obtain the GLS estimator of $\beta$ assuming that the $\sigma^{2}_{i}$'s are known.
 A: This is essentially a "groupwise heteroskedastic" model, and it doesn't matter whether the three $Y_i$'s and the three $X_i$'s may contain different random variables (and not just different realizations of the same random variables). All that matters is that
a) The coefficient vector is assumed identical,
b) The regressor matrices have the same column dimension, and
c) The $\mathbf U_1, \mathbf U_2, \mathbf U_3$ are assumed independent.  
Then it is business as usual.  
To solve it, we stack our data set per equation, i.e. all observations pertaining to $Y_1$ equation first, etc, obtaining a $T\times K$ regressor matrix, $T= T_1+T_2+T_3$. We have the model
$$\mathbf y = \mathbf X\beta + \mathbf u \tag{1}$$
Mind you, we don't really care if each equation has different regressors
The model is heteroskedastic, with $T \times T$ conditional covariance matrix
$$E(\mathbf u \mathbf u' \mid \mathbf X)=\mathbf V = \text{diag}(\sigma^2_1,...,\sigma^2_2,...,\sigma^2_3,...) \tag{2}$$ with appropriate sub-dimensions. This matrix is symmetric and positive definite, so, as a general matrix-algebra result, there exists a decomposition
$$\mathbf V^{-1} = \mathbf C'\mathbf C = \text{diag}(1/\sigma^2_1,...,1/\sigma^2_2,...,1/\sigma^2_3,...) \tag{3}$$
which also gives
$$\mathbf V = \mathbf C^{-1} (\mathbf C')^{-1} \tag{4}$$
So in our case 
$$\mathbf C = \text{diag}(1/\sigma_1,...,1/\sigma_2,...,1/\sigma_3,...) \tag{5}$$
where $\mathbf C$ is non-singular $T \times T$ square matrix. Consider the linear transformation of the model by premultiplication by $\mathbf C$:
$$\mathbf {\tilde y} = \mathbf C \mathbf y, \;\; \mathbf {\tilde X}= \mathbf C \mathbf X,\;\; \mathbf {\tilde u}= \mathbf C \mathbf u$$
This transformation makes the error term conditionally homoskedastic, since
$$E[\mathbf {\tilde u}\mathbf {\tilde u}'\mid \mathbf X] = \mathbf C E(\mathbf u\mathbf u'\mid \mathbf X) \mathbf C' = \mathbf C \mathbf V \mathbf C' = \mathbf C \mathbf C^{-1} (\mathbf C')^{-1} \mathbf C' = \mathbf I_T$$
and one can verify that all other properties of the linear regression model are satisfied. So running OLS on the transformed model is optimal (in the usual sense).
We will obtain
$$\hat \beta^{(tr)}_{OLS} = \left(\mathbf {\tilde X}'\mathbf {\tilde X}\right)^{-1}\mathbf {\tilde X}'\mathbf {\tilde y}  = \left(\mathbf X'\mathbf C' \mathbf C\mathbf X\right)^{-1}\mathbf X' \mathbf C' \mathbf C\mathbf y = \left(\mathbf X'\mathbf V^{-1}\mathbf X\right)^{-1}\mathbf X'\mathbf V^{-1}\mathbf y $$
But this is equivalent to run GLS on the un-transformed model,
$$\beta_{GLS} = \left(\mathbf X'\mathbf V^{-1}\mathbf X\right)^{-1}\mathbf X'\mathbf V^{-1}\mathbf y $$
and its conditional variance is easy to derive as
$$\text{Var}\left[\hat \beta_{GLS}\right] = \left(\mathbf X'\mathbf V^{-1}\mathbf X\right)^{-1}$$
