# Why is the Breusch-Pagan test significant on simulated data designed not to be heteroscedastic?

I'm testing the residuals of a linear regression using Breusch-Pagan Test to detect Heteroscedasticity.

This is the plot of the residuals:

and this is the R code:

> library(lmtest)
>
> mod <- lm(rnorm(1000)~1)
>
> bptest(mod)

studentized Breusch-Pagan test

data:  mod
BP = 0, df = 0, p-value < 2.2e-16


Could someone tell me why it rejects the null hypothesis of homoscedastic errors?

The plot doesn't look heteroscedastic.

EDIT:

However the plot is an example, I have two list of prices (priceA and priceB), I need to check if the residuals generated by a linear regression of these two list: lm(priceA ~ priceB + 0) I need zero intercept are homescedastic or not. Could someone give me a small example? The length of each price list is 750.

EDIT:

I also get:

BP = 67.4362, df = 1, p-value < 2.2e-16

with this chart

Does it be homoscedastic? I have plotted the residuals.

@Wolfgang, I get this result following the example you posted.

> summary(mod)$r.squared * 750 [1] 681.0114  - Unless the index is some meaningful quantity (e.g., directly related to a time variable), this plot doesn't tell you very much. You might start by looking at a plot of the residuals versus the fitted values. – cardinal Aug 27 '11 at 12:12 @cardinal I need a method to check if my residuals are homo/hetero scedastic. My goal is not to see the result plotting something. I plot the chart above because I get a strange p.value, but I Only need a response if the model is homoscedastic. Could this method good for me? i have a matrix with prices so i do: lm(prices[,1] ~ prices[,2]+0) I need zero intercept....if I put +0 I always get p-value < 2.2e-16 – Dail Aug 27 '11 at 12:32 Why do you want the intercept to be 0? – Wolfgang Aug 27 '11 at 15:08 @Wolfgang, because if one stock ir rising +2 also the other should be the same. – Dail Aug 27 '11 at 21:30 Unfortunately, I don't really understanding your reasoning there. If you want to examine whether the variance in the residuals (when regressing prices[,1] on prices[,2]) changes as a function of prices[,2], then I would suggest you do bptest(lm(prices[,1] ~ prices[,2])). – Wolfgang Aug 28 '11 at 12:34 show 3 more comments ## 3 Answers No, the data are not heteroscedastic (by way of how you simulated them). Did you notice the 0 degrees of freedom of the test? That is a hint that something is going wrong here. The B-P test takes the squared residuals from the model and tests whether the predictors in the model (or any other predictors you specify) can account for substantial amounts of variability in these values. Since you only have the intercept in the model, it cannot account for any variability by definition. Take a look at: http://en.wikipedia.org/wiki/Breusch-Pagan_test Also, make sure you read help(bptest). That should help to clarify things. One thing that is going wrong here is that the bptest() function apparently does not test for this errant case and happens to throw out a tiny p-value. In fact, if you look carefully at the code underlying the bptest() function, essentially this is happening: format.pval(pchisq(0,0), digits=4)  which gives "< 2.2e-16". So, pchisq(0,0) returns 0 and that is turned into "< 2.2e-16" by format.pval(). In a way, that is all correct, but it would probably help to test for zero dfs in bptest() to avoid this sort of confusion. EDIT There is still lots of confusion concerning this question. Maybe it helps to really show what the B-P test actually does. Here is an example. First, let's simulate some data that are homoscedastic. Then we fit a regression model with two predictors. And then we carry out the B-P test with the bptest() function. library(lmtest) n <- 100 x1i <- rnorm(n) x2i <- rnorm(n) yi <- rnorm(n) mod <- lm(yi ~ x1i + x2i) bptest(mod)  So, what is really happening? First, take the squared residuals based on the regression model. Then take$n \times R^2$when regressing these squared residuals on the predictors that were included in the original model (note that the bptest() function uses the same predictors as in the original model, but one can also use other predictors here if one suspects that the heteroscedasticity is a function of other variables). That is the test statistic for the B-P test. Under the null hypothesis of homoscedasticity, this test statistic follows a chi-square distribution with degrees of freedom equal to the number of predictors used in the test (not counting the intercept). So, let's see if we can get the same results: e2 <- resid(mod)^2 bp <- summary(lm(e2 ~ x1i + x2i))$r.squared * n
bp
pchisq(bp, df=2, lower.tail=FALSE)


Yep, that works. By chance, the test above may turn out to be significant (which is a Type I error since the data simulated are homoscedastic), but in most cases it will be non-significant.

-
I'm consufed at the moment. I have two vectors with many prices (stock prices). Can I test if their residuals are constant with this method or not? could you give me a small example? thank you so much! –  Dail Aug 27 '11 at 11:56
When the variance in the residuals is heteroscedastic, it means that the variance changes as a function of one of the predictor variables, as a function of some other variable not included in the model, and/or it changes over time. For example, as Karl already pointed out, you could use the index of the observations as an explanatory variable to test whether the variance changes over time. –  Wolfgang Aug 27 '11 at 12:27
with index do you mean a vector like: c(1:N) where N is the length of my residuals vector? –  Dail Aug 27 '11 at 13:07
Yes. Like Karl showed. –  Wolfgang Aug 27 '11 at 15:05
:Wolfgang BP tests the hypthesis that all the "information" has been extracted from the predictors. That is a far cry from thoroughly testing for non-constant variance as non-constant variance can arise from a number of sources. Please see my comment at stats.stackexchange.com/questions/14842/… –  IrishStat Aug 28 '11 at 1:24

The results are not meaningful without some predictor (note df=0). Heteroscedastic means that the variance is not constant, but not constant with respect to what? Perhaps you have in mind the index (order of measurement)? Then you should do

y <- rnorm(1000)
x <- 1:1000
mod <- lm(y~x)
bptest(mod) # I get p=0.59


If you just have a vector of numbers, there's not a whole lot of meaning to the question "Is the variance constant?" For example, consider a mixture of two normal distributions with different variances:

v <- sample(c(1,10), 100, repl=TRUE)
y <- rnorm(100, 0, v)


$\text{var}(y|v)$ is not constant, but depends on $v$. But unconditionally, $\text{var}(y)$ is just a number.

-
what are the y elements? –  Dail Aug 27 '11 at 11:48
I was following from your example; edited to be more clear –  Karl Aug 27 '11 at 12:09
the last example seems constant, I get: BP = 0.0133, df = 1, p-value = 0.9083 –  Dail Aug 27 '11 at 12:36
however as I told, I need this test to check the residuals generated from a regression with two list of prices lm(pricesA~PricesB) How can I use this method to check their residuals? –  Dail Aug 27 '11 at 13:36
what about if i do: newMod <- lm(resid(mod)~c(1,750)) where 750 are the length of the vector of prices. Then, I tried lm(resid(mod)~rep(1,750)) I always get a p-value = 1 surelly is a stupid questione but, could you explain the reason? Thank you so much! –  Dail Aug 27 '11 at 21:40