How to model heteroscedasticity and get the coefficients?

Question

I would like to model heteroscedasticity or better, the standard deviation of the mean as a function of x and retrieve the coefficients that says something about how standard deviation behaves with respect to x.

here is some data:

n <- 1000
x <- runif(n)
y <- 1 - 2*x + x*rnorm(n)
plot(y~x)
r <- lm(y~x)
qplot(x,r$residuals)

so the standard deviation of the mean (linear regression line of -2*x ) seems to be a function of x + x^2; how can I retrieve the coefficients of the standard deviations in R?

Answer 1

First, you need to formally define the model you want to estimate. Based on your example, you'll want to set a model like this: $$Y_i = \beta_0 +\beta_1 X_i + \epsilon_i$$ where $\epsilon_i \sim N(0, X_i^2 \sigma^2)$ and each error is iid.

Similar to a Simple Linear Regression model, there are three parameters to estimate, $\beta_0, \beta_1$ and $\sigma^2$. However, the model we've set up is a special case of a Generalized Linear Model. This means that estimating parameters can't be done as easily - at least, there are no simple formulas for them. The most common approach used to estimate the parameters is to emply Maximum Likelihood Estimation. In other words, find the three values which maximize the likelihood function.

This can easily be implemented in R:

# Define the data
n <- 1000
x <- runif(n)
y <- 1 - 2*x + x*rnorm(n)

# Define the negative log-likelihood function
LL <- function(theta) {
  fvec <- dnorm(y - (theta[1] + theta[2] * x), sqrt(x^2 * theta[3]^2))
  loglike <- sum(-log(fvec))

  return(loglike)
}

# Define initial parameter guesses and then optimize the parameters
theta0 <- c(0, 0, 0)
opt0 <- optim(theta0, LL)
theta1 <- opt0$par
theta1

Unfortunately, this numerical estimator won't have nice statistical properties like unbiasedness or even consistency. But it'll do a passable job of fitting the model and allow you to avoid overthinking things. However, if you want the best possible estimator for your model, check out https://en.wikipedia.org/wiki/Heteroscedasticity-consistent_standard_errors.