I am studying the factors influencing the annual salary for employees at a undisclosed bank. The regression model that I have decided to employ is as follows:
\begin{equation} Y_{k}=\beta_{1}+\beta_{2}E_{k}+\beta_{3}D_{gk}+\beta_{4}D_{mk}+\beta_{5}D_{2k}+\beta_{6}D_{3k}+\varepsilon_{k} \end{equation} where $Y_{k}$ is the logarithm of annual salary, $E$ is the number of years of education, $D_{g}$ is a gender dummy, $D_{m}$ is a minority dummy, and where
\begin{equation}
D_{2}=\begin{cases} 1 &\text{Custodial job} \\ 0 & \text{Otherwise} \end{cases}
\end{equation}
and
\begin{equation}
D_{3}=\begin{cases} 1 &\text{Management job} \\ 0 & \text{Otherwise} \end{cases}
\end{equation}
As you know, whenever one deals with GLS, $\Omega$ will almost surely be unknown and thus have to be estimated. In general there are $\frac{n(n+1)}{2}$ parameters to be estimated, which makes it pretty impossible to come up with a viable estimation out of $n$ observations. This is usually counteracted by imposing some structure on $\Omega$.
In my case, I would like to make the assumption that the disturbance terms $\varepsilon_{k}$ in the above regression model have variance $\sigma_{i}^{2}$ for $i=1,2,3$, according to whether the $i$-th employee has a job in category 1,2, or 3 respectively. Now, we may introduce the transformations $\gamma_{1}=\log (\sigma_{1}^{2}),\gamma_{2}=\log(\sigma_{2}^{2}/\sigma_{1}^{2})$, and $\gamma_{3}=\log(\sigma_{3}^{2}/\sigma_{1}^{2})$ so as to enable us to formulate the following model for
\begin{equation} \sigma_{k}^{2}= \exp \{ \gamma_{1}+\gamma_{2}D_{2k}+\gamma_{3}D_{3k} \} \end{equation} Since $\hat{\beta}_\rm{OLS}$ is a consistent estimate of $\beta$, even under the assumption of heteroscedasticity, we have that $\hat{\beta}_{\rm OLS} \xrightarrow[]{p}\beta$ as the number of observations increase. We may therefore argue that $e_{k}^{2} \approx \sigma_{k}^{2}$, and so we can regress upon information that we already possess.
Summary of procedure
(1) Calculate the OLS estimate.
(2) Calculate the OLS residual $\textbf{e}=\textbf{Y}-\textbf{X}\hat{\beta}$
(3) Calculate the OLS estimate of $\gamma$ from $e_{k}^{2}=f_{\gamma}(Z_{k})+\overline{\varepsilon}_{k}$.
(4) Calculate the FGLS estimate as the GLS estimate with $\hat{\Omega}=\Omega(\hat{\gamma})$ in place of $\Omega$.
What I would like to know is whether or not one can perform this estimation using a known function in R, say
gls? If the answer is yes, then how exactly should I write to ensure that that my heteroscedasticity assumption is taken into account? Thanks for taking the time! Have a great day!
