# MLE of regression coefficients with non-constant variance

In the simplest case of the statistical view of linear regression, we have that

$$y = f(\mathbf{x};\mathbf{w}) + \nu, \text{ where }\nu \sim \mathcal{N}(0,\sigma^2)$$

where $\mathbf{x}$ is an input vecto, $\mathbf{w}$ the vector of the coefficients of the linear function $f$ and $\nu$ is so called "white noise". Finding the MLE for $\mathbf{w}$ is pretty much a standard procedure, no real tricks involved.

What I am wondering about is what happens in the following scenario

$$y = f(\mathbf{x};\mathbf{w}) + \nu, \text{ where }\nu \sim \mathcal{N}(0,\sigma_{\mathbf{x}}^2)$$

namely, when the variance $\sigma^2$ is a function of the input vector $\mathbf{x}$. In the case of constant variance, we create the likelihood function

$$\displaystyle p(\mathbf{y}|\mathbf{x}_i; \mathbf{w}, \sigma^2) = \prod_{i}^{N}p(y_i|\mathbf{x}_i; \mathbf{w}, \sigma^2)$$

and maximize it using standard techniques. Something like this does not seem to be possible under the new model with non-constant variance, since the variance is not really a fixed parameter of the model but rather a function of the input. It is also important to stress that we do not know the value of $\sigma_{\mathbf{x}}^2$ on each $\mathbf{x}$, we only know that it varies with $\mathbf{x}$.

Now, if we happen to know the variance on each input vector $\mathbf{x}_i$, we (I think) can use the exact same methodology as the case with constant variance, namely maximize:

$$\displaystyle p(\mathbf{y}|\mathbf{x}_i; \mathbf{w}, \sigma) = \prod_{i}^{N}p(y_i|\mathbf{x}_i; \mathbf{w}, \sigma^2_{\mathbf{x}_i})$$

I would appreciate comments on the following two questions:

1. Is the first problem solvable without using Bayesian methods or some very advanced estimation theory?
2. Am I missing something in the proposed solution for the case where we know the variance for each input vector?
• You could start out by considering simpler variance functions such as $\sigma^2_{X}=X_1 \sigma^2$. Not sure you can get too far considering all possible variance functions. A Bayesian method would also need some specificity. – JimB Oct 13 '17 at 2:03
• Is this different from generalized least squares? – Dave Feb 21 '20 at 3:20

You are using a linear model and assume normally distributed errors. With constant variance you assume that $y_i$ is normally distributed but with different means for each individual, and that mean depends on some explanatory variables $x$. So you can write the distribution of $y_i$ as $y_i \sim N(\mu_i, \sigma)$. This becomes solvable if we replace $\mu_i$ with $\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots+ \beta_k x_{ki}$ (as long as $k<N$.

You can follow the same logic with the residual variance. Now both $\mu_i$ and $\sigma_i$ differ between individuals, and both depend on explanatory variable $x$ and $z$. $x$ and $z$ can be the same, overlapping, completely different. So $y_i$ is now distributed $y_i \sim N(\mu_i,\sigma_i)$ and $\mu_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots+ \beta_k x_{ki}$, and $\sigma_i = \gamma_0 + \gamma_1 x_{1i} + \gamma_2 x_{2i} + \cdots+ \gamma_k x_{ki}$. This might work, but this specification does not rule out negative variances. So it is safer to specify $\ln(\sigma_i) = \gamma_0 + \gamma_1 x_{1i} + \gamma_2 x_{2i} + \cdots+ \gamma_k x_{ki}$

For example, this is how I would implement such a model on the fly in Stata. mlexp expects the log likelihood function, which is the log of the distribution of $y_i$, so that explains the ln(normalden()) part. The first argument of normalden() is the dependent variabe, so that explains mpg. The second and third arguments of normalden are the mean and the variance. Here {some_name:var1 var2 _cons} is short for $\beta_0 + \beta_1 var_1 + \beta_2 var_2$. The some_name is used for labeling, which is expecially important when the same variable is used in both the mean and the variance equation.

. sysuse auto
(1978 Automobile Data)

. mlexp (ln(normalden(mpg , {xb:weight foreign _cons}, exp({lnsigma:foreign _cons}))))

initial:       log likelihood =     -<inf>  (could not be evaluated)
feasible:      log likelihood = -6441.8168
rescale:       log likelihood = -368.12496
rescale eq:    log likelihood =  -257.3948
Iteration 0:   log likelihood =  -257.3948  (not concave)
Iteration 1:   log likelihood = -231.78128  (not concave)
Iteration 2:   log likelihood = -226.86965  (not concave)
Iteration 3:   log likelihood = -224.24518  (not concave)
Iteration 4:   log likelihood = -198.21111
Iteration 5:   log likelihood = -192.59299
Iteration 6:   log likelihood = -183.30713
Iteration 7:   log likelihood = -183.25379
Iteration 8:   log likelihood = -183.25373
Iteration 9:   log likelihood = -183.25373

Maximum likelihood estimation

Log likelihood = -183.25373                     Number of obs     =         74

------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
xb           |
weight |  -.0061135   .0004522   -13.52   0.000    -.0069998   -.0052272
foreign |  -1.175046   1.211334    -0.97   0.332    -3.549217    1.199126
_cons |   40.10603   1.532606    26.17   0.000     37.10217    43.10988
-------------+----------------------------------------------------------------
lnsigma      |
foreign |   .8022486   .1808845     4.44   0.000     .4477215    1.156776
_cons |    .818957   .0982276     8.34   0.000     .6264344     1.01148
------------------------------------------------------------------------------