I am testing my model using the Breusch-Pagan Test, but have not been able to find anything online regarding how to calculate it for an ARIMA Model.

My AR1 Model is:

ar1 <- arima(x = datasource[, "Y"], order = c(1, 0, 0), 
    xreg = datasource[, c("X1", "X2", "X3", "X4", "X5")])
  • 3
    $\begingroup$ Welcome to SE. I edited your post so that the code is readable (just space it over 4 positions, or use the "{}" icon in the compose window). Also, this might be better suited for StackExchange, as it's primarily about coding/programming and not about a statistical issue. $\endgroup$
    – Russ Lenth
    May 18, 2015 at 22:43
  • $\begingroup$ You can apply any script or function to the residuals ar1$residuals. $\endgroup$ Aug 3, 2017 at 10:45

2 Answers 2


Your question makes me wonder why the Breusch-Pagan (BP) test is not available for the output from a fitted ARIMA model. Is the BP test applicable to ARIMA models? (This question fits better within the scope of this site.)

I don't see theoretical reasons that would invalidate the test in the context of an ARIMA model. The purpose and interpretation of the test are nonetheless more appropriate in the context of a regression model rather than in a model for the autocorrelation structure of a time series model.

You did not mention it in the question but I noticed in the sample code that you are including exogenous regressors in the ARIMA model, X1,...,X5. In this case, the BP test can be interesting in order to test whether there is a relationship between these regressors and the variance of the error.

Some ideas on how to code or obtain the BP test statistic with an ARIMA model in R (the software that you are using according to your sample code):

If you are using an AR model, then you can fit it as a linear regression model where lags of the dependent variable are included as regressors (along with the other explanatory variables). Then you can use the function ncvTest in package car or bptest in package lmtest. Example for an AR(1) model:

y <- datasource[, "Y"]
n <- length(y)
xreg <- as.matrix(datasource[, c("X1", "X2", "X3", "X4", "X5")])
fit1 <- lm(y[2:n] ~ 0 + xreg[2:n,] + y[1:(n-1)])
# alternatively the "dynlm" interface can be used
fit <- dynlm(y ~ xreg + L(y,1))
# lmtest::bptestreturns the Breusch and Pagan test statistic

Be aware of the differences between the linear regression model and the specification of AR model with exogenous regressors. See this post for details.

If you want to use the ARIMA model specification or if the model includes an MA part, you could include as regressors in the linear model lags of the residuals obtained in a previous step, but it is easier to implement the test upon the residuals of the fitted ARIMA model:

ar1 <- arima(y, order=c(1,0,0), xreg=xreg)
e <- residuals(ar1)
fitaux <- lm(e ~ 1 + xreg)
df <- ncol(xreg) - 1
bp.statistic <- (length(e) - df) * summary(fitaux)$r.squared
bp.pvalue <- pchisq(bp.statistic, df, lower.tail=FALSE)
  • $\begingroup$ Thanks a lot for your feedback! It helped me get the BP statistic. Now in the case of an AR(2) model, would the formulas you provided be different? I'm guessing n-1 would turn into n-2 and 2:n would become 3:n? Apologies for my ignorance in statistics. $\endgroup$
    – Ray
    May 21, 2015 at 18:07
  • $\begingroup$ For an AR(2) you should change the indices, lm(y[3:n] ~ 0 + xreg[3:n] + y[2:(n-1)] + y[1:(n-2)]). It may be easier and safer to use the dynlm interface. $\endgroup$
    – javlacalle
    May 22, 2015 at 7:16
  • $\begingroup$ Notice also that I have fixed in the answer the indices for xreg, they should be the same that for the dependent variable (assuming they are not lagged regressors). $\endgroup$
    – javlacalle
    May 22, 2015 at 7:17

i think this can help:

adapted from: https://github.com/cran/car/blob/master/R/ncvTest.R

ncvTest.arima <- function(model, ...) {
    data <- getCall(model)$data
    model <- model
    sumry <- summary(model)
    residuals <- residuals(model)

    S.sq <- (length(residuals)-2)*(model$sigma2)/length(residuals)
    .U <- (residuals^2)/S.sq
    mod <- lm(.U ~ fitted.values(model))
    varnames <- "fitted.values"
    var.formula <- ~ fitted.values
    df <- 1
    SS <- anova(mod)$"Sum Sq"
    RegSS <- sum(SS) - SS[length(SS)]
    Chisq <- RegSS/2
    result <- list(formula=var.formula, formula.name="Variance",ChiSquare=Chisq, Df=df, p=pchisq(Chisq, df, lower.tail=FALSE), test="Non-constant Variance Score Test")
    class(result) <- "chisqTest"


Not the answer you're looking for? Browse other questions tagged or ask your own question.