2

How Breusch-Pagan can't reject the null for a series like that?

> x = c(rnorm(10), rnorm(100,sd=10), rnorm(100,sd=25))
> mod = lm(x[-1]^2~x[-210]^2)
> plot(mod$res,type='l')
> bptest(mod)

        studentized Breusch-Pagan test

data:  mod 
BP = 1.1085, df = 1, p-value = 0.2924

The variance change drastically, could someone explain the reason?

enter image description here

Nick Sabbe
  • 12,119
  • 2
  • 35
  • 43
Dail
  • 2,147
  • 12
  • 40
  • 54
  • Maybe you should look at [some examples](http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/lmtest/html/bptest.html) and emulate them closely. – whuber Dec 09 '11 at 14:07
  • @whuber I do not need to see at the examples. I'm only asking why BP test can't reject the null with this kind of residuals. BP test works on the residuals of the model... I plot those points, and look drastically different from the beginning of the series. How is it possible that those residuals are homoscedastics? – Dail Dec 09 '11 at 14:23
  • 2
    You definitely do need to study the examples, because your model does not do what you think it does. *First* replicate the examples in the documentation, *then* put your data into the form used in the examples and re-apply the examples. – whuber Dec 09 '11 at 14:31

1 Answers1

5

Suppose that you have the model $y_i = \beta_0 + \beta_1 x_i + \epsilon_i$. You want to test whether $\text{Var}(\epsilon_i \mid x_i) = \sigma^2$, that is, constant across $i$. To test this, you need to write $\text{Var}(\epsilon_i \mid x_i) = h(x_i)$, there $h(\cdot)$ is some function, by default for the BP test given by $\alpha_0 + \alpha_1 x_i$.

Some of your $x$ variables have bigger variances than the others, but unless there is some other variable that lets you predict which $x$'s those are, you won't pick up on that fact.

@whuber is noting that your example is not set up like this. First, you presumably expect your x to be your outcome---that's the variable that has heteroskedasticity. But it's not heteroskedasticity that you can explain. So create a new variable z that, when represented as a factor in the regression, is a set of dummy variables for membership in the three possible groups in your sample:

set.seed(1109)
x <- c(rnorm(10), rnorm(100,sd=10), rnorm(100,sd=25))
z <- c(rep(1, 10), rep(2, 100), rep(3, 100))
z <- as.factor(z)

The regression of interest here would be lm(x ~ z). The bptest() would be

bptest(x ~ z)

studentized Breusch-Pagan test

data:  x ~ z 
BP = 37.1871, df = 1, p-value = 1.073e-09

Since x has a heteroskedasticity pattern that is predictable using covariates, we can detect it. But if we didn't have that covariate, we wouldn't detect the heteroskedasticity.

To sum up in a different way, the BP test only has power against (that is, it can only detect) heteroskedasticity that is predictable using the covariates. If you can't do that, then you can't detect the heteroskedasticity.

Charlie
  • 13,124
  • 5
  • 38
  • 68
  • Thank you for your reply, the problem is that you subdivide the serie is three different sample, as the example I made, BUT it is just and example. I don't know WHERE there is a change in variance (where you put the rep()) – Dail Dec 09 '11 at 16:44
  • It does not make sense to me because you KNOW the series and so you created the factor for that series...but if you only have the serie x how do you detect heteroscedasticity? – Dail Dec 09 '11 at 18:20
  • The only heteroskedasticity that you can detect is heteroskedasticity that is predictable by the $x$ variables. That's the theory. If you are worried about heteroskedasticity that you can't detect (that is, you run the BP test, can't reject the null, but think that there is heteroskedasticity), then forget about the test and just use robust standard errors and move on. – Charlie Dec 09 '11 at 18:39
  • thank you about the explanation, what do you mean with "robust standard errors" ? could you give me a link (R related) to test it? thank you so much! – Dail Dec 09 '11 at 19:16
  • Remember, heteroskedasticity doesn't bias your estimates; you just get the wrong standard errors. Please go through those notes that I sent you---I discuss robust standard errors there. In R, use the function `vcovHC()` in the package `sandwich` (http://cran.r-project.org/web/packages/sandwich/index.html). – Charlie Dec 09 '11 at 20:37
  • @dail, try Charlie's example with `z – naught101 Apr 28 '12 at 14:30