Consider a vector of observations $\mathbf{Y}$ that can be modeled as \begin{equation} \mathbf{Y} \sim \mathcal{N}( \mathbf{H}\boldsymbol{\beta} , \boldsymbol{\Sigma} ) \end{equation} with $\mathbf{H}$ a linear parametric model, and $\boldsymbol{\Sigma}$ a diagonal matrix whose elements are given by $a_1 \mathbf{H}\boldsymbol{\beta} + a_2$, and where the constants $a_1$ and $a_2$ are known.
To properly estimate $\boldsymbol{\beta}$, and since the variance is heteroscedastic but depends on the regression parameters, a first estimation can be obtained applying OLS: $$ \widetilde{\boldsymbol{\beta}}^{\mathrm{ols}} = \left( \mathbf{H}^\mathrm{T} \mathbf{H} \right)^{-1} \mathbf{H}^\mathrm{T} \mathbf{Y} . $$ Then, a rough estimation of the noise variance is obtained from $\widetilde{\boldsymbol{\beta}}^{\mathrm{ols}}$ by: $$ \widetilde{\boldsymbol{\Sigma}}^{\mathrm{ols}} = \mathbf{I}\times (a_1 \mathbf{H}\, \widetilde{\boldsymbol{\beta}}^{\mathrm{ols}} + a_2 ), $$ where $\mathbf{I}$ denotes the identity matrix. This rough estimation of the covariance is thus reused to update the estimation of the expectation using WLS: \begin{equation} \label{eq:estimated_expectation_detect} \begin{cases} \widetilde{\boldsymbol{\beta}} = \left( \mathbf{H}^\mathrm{T} {\widetilde{\boldsymbol{\Sigma}}^{{\mathrm{ols}}^{-1}}} \mathbf{H} \right)^{-1} \mathbf{H}^\mathrm{T} {\widetilde{\boldsymbol{\Sigma}}^{{\mathrm{ols}}^{-1}}} \mathbf{Y} , \\ \widetilde{\boldsymbol{\Sigma}} = \mathbf{I}\times (a_1 \mathbf{H}\, \widetilde{\boldsymbol{\beta}} + a_2 ). \end{cases} \end{equation}
We can prove that the estimation $\widetilde{\boldsymbol{\beta}}^{\mathrm{ols}}$ is unbiased. However, is the estimation $\widetilde{\boldsymbol{\beta}}$ unbiased?
Another question is, is there a better way to estimate $\boldsymbol{\beta}$?
