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)

mydata <- read.table("hprice1.csv", sep=",",  header = TRUE, na.strings = ".")

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))