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