# Heteroscedastic error estimation

We have random, heteroscedastic, uncorrelated effects $u_i \sim N(0,\sigma_i^2)$. We use that $u=[u_1,u_2,\ldots,u_N] \sim N(0, \Sigma_u)$. We observe weighted sums of $u$, we call these residuals $\epsilon$. We denote the weights by $A$, a known matrix, with rows $A_k$ and elements $a_{k,j}$. \begin{align} \epsilon_k = A_k u \sim N(0, \sigma_k^2 A_k A_k^T) = N(0, \sum\limits^N_{j=1} \sigma_j^2 a^2_{k,j} ) \end{align} An important thing to note is that all residuals are independent, thus $\operatorname{cov}(\epsilon_i,\epsilon_j)=0$ for $i\ne j$.

My question is; how do we make (unbiased) estimators for the variance of $u_i$'s? I have tried looking at a matrix representation of the problem: \begin{align} \epsilon \sim N(0, A \Sigma_u A \ \circ I) \end{align} Here we use $B \circ C$ as the element-wise product of matrices $B$ and $C$. Furthermore we can estimate $\operatorname{var}(\epsilon)$ by: $$\hat\Sigma_\epsilon = \operatorname{var}(\epsilon) = E[\epsilon \epsilon^T] = \begin{bmatrix} \epsilon_1^2& & & \\ & \epsilon_2^2& & \\ & & \ddots & \\ & & & \epsilon_M^2 \end{bmatrix} = \epsilon \epsilon^T \circ I$$ Can we now use the following (as a possible solution): \begin{align} \Sigma_\epsilon = \epsilon \epsilon^T \circ I & = A \Sigma_u A \ \circ I \\ \epsilon \epsilon^T & = A \Sigma_u A\\ \Sigma_u &= (A^TA)^{-1}A^T \epsilon \epsilon^T A (A^TA)^{-1} \end{align} For the last part to hold we require that $A^TA$ has full rank $N$ and therefore is $A^TA$ invertible.

## 1 Answer

I think I found a solution method by maximization of the likelihood. We get the following steps: \begin{align} \mathcal{L}(\epsilon |\Sigma_u) &= \frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{|\Sigma_u|}}\exp(-\frac12\epsilon (A \Sigma_u A^T \circ I)^{-1}\epsilon^T)\\ &= \frac{1}{\sqrt{2\pi}} (\prod\limits_{i=1}^N \sigma_i^{-1} ) \exp(-\frac12\sum\limits_{j=1}^{M}\epsilon_j^2 (\sum\limits_{i=1}^{N}a^2_{i,j}\sigma_i^2)^{-1})\\ \log\mathcal{L}(\epsilon |\Sigma_u) &= -\frac{1}{2}\log(2\pi) - \sum\limits_{i=1}^N \log(\sigma_i) - \frac12\sum\limits_{j=1}^{M}\epsilon_j^2 (\sum\limits_{i=1}^{N}a^2_{i,j}\sigma_i^2)^{-1} \\ \frac{\partial \log\mathcal{L}(\epsilon |\Sigma_u)}{\partial \sigma_k}&= -\frac{1}{\sigma_k} +\frac12 \sum\limits_{j=1}^{M}\epsilon_j^2 2a^2_{k,j}\sigma_k (\sum\limits_{i=1}^{N}a^2_{i,j}\sigma_i^2)^{-2})= 0\\ \frac{1}{\sigma_k} &= \sigma_k\sum\limits_{j=1}^{M}\epsilon_j^2 a^2_{k,j}(\sum\limits_{i=1}^{N}a^2_{i,j}\sigma_i^2)^{-2}\\ \sigma_k^2 &= \frac{1}{\sum\limits_{j=1}^{M}\epsilon_j^2 a^2_{k,j}(\sum\limits_{i=1}^{N}a^2_{i,j}\sigma_i^2)^{-2}}\\ \end{align} From here I would suggest to apply fixed-point iteration, iteratively improving out estimates by filling in our previous estimates in the last equation. To continue analytically from here we can maybe apply lagrange multipliers.