# How to model heteroscedasticity and get the coefficients?

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? • You can pass a variance model to the gls function in package nlme. Obviously, only common variance structures are implemented (a quadratic function is not common). – Roland Dec 19 '16 at 7:18 ## 1 Answer 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.