Koenker学习化BP测试

Question

我正在尝试为Breusch Pagan测试编写代码，包括学生化和非学生化。下面的代码可以工作，并与R的lmtest“包中的bptest()函数相匹配。

x = rnorm(50,5,2)
y = 5*x + rnorm(50)
dat <- data.frame(x=x,y=y)

## Get the residuals
mod<-lm(y ~ x, data=dat)
res<-residuals(mod)
dat$res2<-res^2

#SSE = Sum of squared error
sse <- sum(res^2)
sse

mod2<-lm(res2~x,data=dat)
#SSR = Sum of squared regresion = explained sum of squares (ESS)
yhat<-fitted(mod2)
ybar<-mean(dat$res2)
ssr = sum((yhat-ybar)^2)

##BP Test
(.5*ssr)/(sse/nrow(dat))^2
library("lmtest")
bptest(mod,studentize = F)

lambda = (nrow(dat) - 1) / nrow(dat) * var(res^2) / (2 * ((nrow(dat) - 1) / nrow(dat) * var(res))^2)
##Studentized BP Test
(.5*ssr)/(lambda*(sse/nrow(dat))^2)
bptest(mod,studentize = T)

然而，当我尝试简化代码时，它不会返回相同的值。

> n=nrow(dat)
> v=var(res^2)
> (n - 1) / n * v / (2 * ((n - 1) / n * v)^2)
[1] 0.4129854
> lambda
[1] 1.081115

Answer 1

正如评论中指出的，我犯了一个很难看出来的打字错误。

> n=nrow(dat)
> v1=var(res^2)
> v2=var(res)
> (n - 1) / n * v2 / (2 * ((n - 1) / n * v1)^2)
> lambda