Breusch-Pagan hettest following weighted least squares vs R code

[Robert Gulotty]

I am trying to replicate the behavior of STATA's Breusch-Pagan test following running a weighted least squares regression.

sysuse auto
reg price mpg foreign
hettest

This code gives the same result (3.81) as the following code in R (using the autos data exported as a CSV):

autos <- read.csv("~/Downloads/autos.csv")
mod <- lm(price~mpg+foreign, data=autos)
resi <- mod$residuals
sigma2 <- sum(resi^2)/nobs(mod)
bpregSTATA <- lm(resi^2/sigma2~fitted(mod))
0.5 * sum((bpregSTATA$fitted.values-mean(resi^2/sigma2))^2)

However, while I get the same regression output as R when I apply weights, the output of hettest changes relative to the R code.

sysuse auto
reg price mpg foreign [aweight=weight] 
hettest

Produces: 1.06

modwt <- lm(price~ mpg + foreign, weights = weight, data=autos)
resiwt <- modwt$residuals
sigma2wt <- sum(resiwt^2)/nobs(modwt)
bpregSTATAwt <- lm(resiwt^2/sigma2wt ~ fitted(modwt))
0.5 * sum((bpregSTATAwt$fitted.values - mean(resiwt^2/sigma2wt))^2)

Produces: 1.342911

What I would like to know is what hettest is doing following a weighted regression to help understand how to correct my formulas. Thank you for your time.

[FernandoRios]

I think for your R replications you need to include the weights on the estimation of sigma2wt, the regression, and the statistic.

[Robert Gulotty]

Dear Fernando, I appreciate your help, I think you are correct. However it is difficult for me to figure out how to include the weights. Unfortunately, this is not the correct forum for getting help with R. Here is what I have figured out:

wts<-autos$weight/mean(autos$weight)
modwt<-lm(price~mpg+foreign, weights = wts, data=autos)

resiwt<-sqrt(wts)*residuals(modwt)

sigma2wt<-sum(resiwt^2)/nobs(modwt)

ui<-resiwt^2/sigma2wt

bpregSTATAwt<-lm(ui~modwt$fitted.values, weights=wts, data=autos)

0.5 * sum(wts*(bpregSTATAwt$fitted.values-(1/sum(wts))*sum(wts*ui))^2)

The STATA code is a bit obscure to me, as some of it is baked into the e(mss) following the weighted regression.

[Robert Gulotty]

It turns out there are three things to know. First, it is necessary to rescale the weights

wts<-autos$weight/mean(autos$weight)

modwt<-lm(price~mpg+foreign, weights = wts, data=autos)
resiwt<-sqrt(wts)*residuals(modwt)
sigma2wt<-sum(resiwt^2)/nobs(modwt)

Second, the squared residuals need to be unweighted, but normalized by the weighted rss.

ui<-residuals(modwt)^2/sigma2wt

Third it is necessary to extract a weighted TSS from the auxiliary regression to calculate the MSS.

bpregSTATAwt<-lm(ui ~ modwt$fitted.values, weights=wts, data=autos)
.5*(sum(wts*(ui-weighted.mean(ui, wts))^2)-deviance(bpregSTATAwt))

Hopefully this helps someone else in the future.