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.