# Breusch-Pagan test for heteroskedasticity

The program below performs the Breusch-Pagan test step-by-step, then uses bptest from the lmtest package.

In order to use the built-in heteroskedasticity tests in R, you’ll have to install the lmtest package.

install.packages("lmtest")


## References

Problems, examples, and data from Wooldridge, Jeffrey M. Introductory Econometrics: A Modern Approach, 4th ed. South-Western Pub. 2009.

## Listing of the R program.

#
#  30 Oct 2012 D.S.Dixon
#
# Reproduces Example 8.4 in
#
#  Wooldridge, Jeffrey M. Introductory Econometrics: A Modern Approach. 4th ed. South-Western Pub. 2009.
#

# HPRICE1.DES
#
# price     assess    bdrms     lotsize   sqrft     colonial  lprice    lassess
# llotsize  lsqrft
#
#   Obs:    88
#
#   1. price                    house price, $1000s # 2. assess assessed value,$1000s
#   3. bdrms                    number of bedrooms
#   4. lotsize                  size of lot in square feet
#   5. sqrft                    size of house in square feet
#   6. colonial                 =1 if home is colonial style
#   7. lprice                   log(price)
#   8. lassess                  log(assess
#   9. llotsize                 log(lotsize)
#  10. lsqrft                   log(sqrft)

library(car)
library(lmtest)

formula0 <- mydata$lotsize + mydata$sqrft + mydata$bdrms myfit0 <- lm(price ~ lotsize + sqrft + bdrms, data=mydata) output0 <- summary(myfit0) print(output0) # regress the residuals squared against the same model myfit1 <- lm(I(myfit0$residuals^2) ~ lotsize + sqrft + bdrms, data=mydata)
output1 <- summary(myfit1)
print(output1)

# H0: u^2 does not depend on any parameter (e.g. all coefficients are statistically zero)
# This would be the F-test if it's homoskedastistic (which we don't know 'till after we
#   test it) so we use the LM test where LM ~ chi-squared(df = number of parameters in the model)
Rsq <- output1$r.squared N <- length(output0$residuals)
df <- length(myfit0\$coefficients) - 1
LM <- N * Rsq
# P [X > x]
pvalue <- pchisq(LM, df, lower.tail = FALSE)

cat("LM = ", LM, "\n")
cat("df = ", df, "\n")
cat("p = ", pvalue, "\n")

# now use the canned version of Breusch-Pagan from the lmtest package
print(bptest(myfit0))