Inference in linear model with conditional heteroskedasticity

Question

Suppose I observe independent variable vectors $\vec{x}$ and $\vec{z}$ and dependent variable $y$. I would like to fit a model of the form: $$y = \vec{x}^{\top}\vec{\beta_1} + \sigma g\left(\vec{z}^{\top} \vec{\beta_2}\right) \epsilon,$$ where $g$ is some positive-valued twice-differentiable function, $\sigma$ is an unknown scaling parameter, and $\epsilon$ is a zero-mean, unit-variance Gaussian random variable (assumed to be independent of $\vec{x}$ and $\vec{z}$). This is essentially the setup of Koenker's test of heteroskedasticity (at least as far as I understand it).

I have $n$ of observations of $\vec{x}, \vec{z}$ and $y$, and I would like to estimate $\vec{\beta_1}$ and $\vec{\beta_2}$. I have a few problems, though:

1. I am not sure how to pose the estimation problem as something like least squares (I assume there is a well-known trick). My first guess would be something like $$min_{\vec{\beta_1}, \vec{\beta_2}} \left(\sum_{i=1}^n \frac{\left(y_i - \vec{x_i}^{\top}\vec{\beta_1}\right)^2}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)\left(\sum_{i=1}^n \frac{1}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)^{-1},$$ but I am not sure how to solve that numerically (perhaps an iterative quasi-Newton method might do).
2. Assuming I can pose the problem in a sane way and find some estimates $\hat{\beta}_1, \hat{\beta}_2$, I would like to know the distribution of the estimates so that e.g. I can perform hypothesis tests. I would be fine with testing the two coefficient vectors separately, but would prefer some way to test, e.g. $H_0: \vec{w_1}^{\top} \vec{\beta_1} + \vec{w_2}^{\top} \vec{\beta_2} \le c$ for given $\vec{w_1}, \vec{w_2}, c$.

Answer 1

In a slightly more general context with $Y$ an $n$-dimensional vector of $y$-observations (the responses, or dependent variables), $X$ an $n \times p$ matrix of $x$-observations (covariates, or dependent variables) and $\theta = (\beta_1, \beta_2, \sigma)$ the parameters such that $Y \sim N(X\beta_1, \Sigma(\beta_2, \sigma))$ then the minus-log-likelihood is $$l(\beta_1, \beta_2, \sigma) = \frac{1}{2}(Y-X\beta_1)^T \Sigma(\beta_2, \sigma)^{-1} (Y-X\beta_1) + \frac{1}{2}\log |\Sigma(\beta_2, \sigma)|$$ In the OP's question, $\Sigma(\beta_2, \sigma)$ is diagonal with $$\Sigma(\beta_2, \sigma)_{ii} = \sigma^2 g(z_i^T \beta_2)^2$$ so the determinant becomes $\sigma^{2n} \prod_{i=1}^n g(z_i^T \beta_2)^2$ and the resulting minus-log-likelihood becomes $$\frac{1}{2\sigma^2} \sum_{i=1}^n \frac{(y_i-x_i^T\beta_1)^2}{ g(z_i^T \beta_2)^2} + n \log \sigma + \sum_{i=1}^n \log g(z_i^T \beta_2)$$ There are several ways to approach the minimization of this function (assuming the three parameters are variation independent).

• You can try to minimize the function by a standard optimization algorithm remembering the constraint that $\sigma > 0$.
• You can compute the profile minus-log-likelihood of $(\beta_1, \beta_2)$ by minimizing over $\sigma$ for fixed $(\beta_1, \beta_2)$, and then plug the resulting function into a standard unconstrained optimization algorithm.
• You can alternate between optimizing over each of the three parameters separately. Optimizing over $\sigma$ can be done analytically, optimizing over $\beta_1$ is a weighted least squares regression problem, and optimizing over $\beta_2$ is equivalent to fitting a gamma generalized linear model with $g^2$ the inverse link.

The last suggestion appeals to me because it builds on solutions that I already know well. In addition, the first iteration is something I would consider doing anyway. That is, first compute an initial estimate of $\beta_1$ by ordinary least squares ignoring the potential heteroskedasticity, and then fit a gamma glm to the squared residuals to get an initial estimate of $\beta_2$ $-$ just to check if the more complicated model seems worthwhile. Iterations incorporating the heteroskedasticity into the least squares solution as weights might then improve upon the estimate.

Regarding the second part of the question, I would probably consider computing a confidence interval for the linear combination $w_1^T\beta_1 + w_2^T\beta_2$ either by using standard MLE asymptotics (checking with simulations that the asymptotics works) or by bootstrapping.

Edit: By standard MLE asymptotics I mean using the multivariate normal approximation to the distribution of the MLE with covariance matrix the inverse Fisher information. The Fisher information is by definition the covariance matrix of the gradient of $l$. It depends in general on the parameters. If you can find an analytic expression for this quantity you can try plugging in the MLE. In the alternative, you can estimate the Fisher information by the observed Fisher information, which is the Hessian of $l$ in the MLE. Your parameter of interest is a linear combination of the parameters in the two $\beta$-vectors, hence from the approximating multivariate normal of the MLE you can find a normal approximation of the estimators distribution as described here. This gives you an approximate standard error and you can compute confidence intervals. It's well described in many (mathematical) statistics books, but a reasonably accessible presentation I can recommend is In All Likelihood by Yudi Pawitan. Anyway, the formal derivation of the asymptotic theory is fairly complicated and rely on a number of regularity conditions, and it only gives valid asymptotic distributions. Hence, if in doubt I would always do some simulations with a new model to check if I can trust the results for realistic parameters and sample sizes. Simple, non-parametric bootstrapping where you sample the triples $(y_i,x_i,z_i)$ from the observed data set with replacement can be a useful alternative if the fitting procedure is not too time consuming.