Is it compulsory to make transformation to the econometric model in order to have only diagonal elements on variance-covariance matrix of errors?

I need some sharped and advanced advices for the following issue ...

Model and its assumptions

I'm working on the methodology of a two-way error component model. Here is the model: $y_{jis} = x_{jis} \beta + \upsilon_{jis}$. $j$ refers to school, $i$ refers to individual and $s$ to tested area/topic (mathematics or english, ...). $J$ is the number of schools and $S$ is the number of tested fields. $N$ is the number of students among all schools.

$\upsilon_{jis}$ is the error term and can be decomposed as follows: $\upsilon_{jis} = \theta_{j} + \phi_{js} + \epsilon_{jis}$

$\theta_{j}$ is the random effect for students attenting a school $j$, $\phi_{js}$ is the random effect for students attenting a school $j$ for a topic $s$ and $\epsilon_{jis}$ is the traditionnal idiosyncratic error.

In multivariate form it gives: $$Y = X \beta + \upsilon$$ where $\upsilon = R \theta + F \phi + \epsilon$

$R$ and $F$ are matrices that enable to correclty distribute their respective random effect. R: size $= NS \times J$. F: size $= NS \times JS$

$\theta$ follows a multivariate normal with mean $0$ and a variance covariance matrix $= \tau I_{J}$. $\tau$ is a scalar and $I_{J}$ is the identity matrix (size: $J \times J$).

$\phi$ follows a multivariate normal with mean $0$ and a variance covariance matrix $= \gamma I_{JS}$. $\gamma$ is a scalar and $I_{JS}$ is the identity matrix (size: $J S \times J S$).

$\epsilon$ follows a multivariate normal with mean $0$ and a variance covariance matrix $= \sigma I_{NS}$. $\sigma$ is a scalar and $I_{JS}$ is the identity matrix (size: $N S \times N S$).

Estimation of unknown parameters

I have to make estimation of the different variances: $\tau$, $\gamma$ and $\sigma$ before running FGLS method. I work step by step to delete the different random effect. First I made a within regression to get rid of $\phi$ effect. Since $\phi$ is an effect for a given school, that within operator ($W$) will automatically make the $\theta$ effect vanishing. $\epsilon$ is the only survivor and I can use the methodology in Hayashi to calculate the expectation of $( (We)' We | X) = . ..$ in order to correct with the right degree of freedom. It is easy since the variance of that element is only diagonal, it will then lead to calculate the trace ....

Now, I only want to make the term theta disappear with another within operator ($W_{2}$). By construction, it deletes $\theta$, so there are $\phi$ and $\epsilon$ left. Once again I try to perform the expectation of $( (W_{2}e)' W_{2}e | X)$ where $e$ is composed now of $\phi$ and $\epsilon$. The variance of that vector is not diagonal, thus we don't have a beautiful formula for the estimated variance of $\phi$ ... I carefully follow the good methodology but I don't have a good-looking formula.

Is it compulsory to make transformation to the model in order to have only diagonal elements on variance-covariance matrix of errors?

I can also provide more details about the development of Hayashi that I reuse ...