# Rectify Heteroskedasticity in R?

library(lmtest)
library(MASS)
library(car)
library(caret)
x      = rep(1:100)
a      = 10
b      = 2
sigma2_2 = x*10
eps    = rnorm(x, mean=0, sd=sqrt(sigma2_2))
y2      = a+b*x + eps
model2    = lm(y2 ~ x)

sigma2_3 = x^2*10
eps    = rnorm(x,mean=0,sd=sqrt(sigma2_3))
y3      = a+b*x + eps
model3    = lm(y3 ~ x)

par(mfcol=c(1,2))
plot(x,y2,main="mild heteroscedasticity")
abline(coef(model2), col="red")
plot(x,y3,main="severe heteroscedasticity")
abline(coef(model3), col="red")
par(mfcol=c(1,1))

bptest(model2)
bptest(model3)

model2_new <- lm(y2~x, weights = 1/sqrt(x*10))
summary(model2_new)
plot(model2_new)
bptest(model2_new)

model3_new <- lm(y3~x, weights = 1/x)
summary(model3_new)
plot(model3_new)
bptest(model3_new)

y2BCMod <- BoxCoxTrans(y2)
print(y2BCMod)
y2_new=predict(y2BCMod, y2)
Mod2_bc <- lm(y2_new ~ x)
bptest(Mod2_bc)

y3BCMod <- BoxCoxTrans(y3+450)
print(y3BCMod)
y3_new=predict(y3BCMod, y3)
Mod3_bc <- lm(y3_new ~ x)
bptest(Mod3_bc)


I want to compare the box-cox transformation and WLS method in solving the heteroskedasticity problem, but no matter how I simulate the data, it seems the bptest always rejected the null hypothsis. I am confused. Is there a problem in choosing the weights or other problem in my code?

• You are simulating models with strong heteroskedasticity, and the bptest rejects the null of homoskedasticity, as it should. What is wrong? Please extend your post by including plots and some output from your code, and explain where you are surprised. – kjetil b halvorsen Jan 8 at 12:33
• Box-Cox doesn't attempt to rectify heteroscedasticity, does it? – Sal Mangiafico Jan 9 at 15:05

The problem you are going through is the same one I had when reading about this topic in most introductory books. Most books indicate that WLS will address heteroskedasticity by "eliminating" it from the model. This is technically done by using weights equal to 1 over the source of the heteroskedasticity. The alternative is to transform all data accordingly. What the books do not mention, but we expect, is for heteroskedasticity to now be eliminated. It happens to be that only with the transformed model can one obtain this result. For instance, try the following code

   y2t=y2/(10*x)^.5
xt=x/(10*x)^.5
ct=1/(10*x)^.5
model2_new=lm(y2~x, weights = 1/(x*10))
model2_new2=lm(y2t~0+ct+xt)
summary(model2_new)
summary(model2_new2)
bptest(model2_new)
bptest(model2_new2)


You will see that the one with transformed variables does indeed suggest there is no heteroskedasticity in the model HTH

• Thank you very much.What if I use robust standard error?Do I need to use bptest after using that? – zhi li Jan 14 at 3:15
• That is a completely different approach to address heteroskedasticity. Instead of "correcting" for heteroskedasticity, Robust standard errors "account" for heteroskedasticity by estimating the coefficient standard errors assuming some unknown form of heteroskedasticity. This means. a) if you do a BPTEST, you will still find heteroskedasticity. b) but your standard errors will be correctly estimated, so that t-stats and F-Stats are valid. – Fcold Jan 14 at 13:40
• Also. Generally speaking, the "right" procedure should be: 1. test for heteroskedasticity and 2. if you find evidence of it correct it using WLS or Robust standard errors. – Fcold Jan 14 at 13:44