# Breusch-Pagan Test for ARIMA Model in R [closed]

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")])


## closed as off-topic by Michael Chernick, Stephan Kolassa, Peter Flom♦Aug 3 '17 at 11:21

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Michael Chernick, Stephan Kolassa, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

• 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. – rvl May 18 '15 at 22:43
• You can apply any script or function to the residuals ar1$residuals. – Stephan Kolassa Aug 3 '17 at 10:45 ## 2 Answers 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 is 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 <- 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 require(dynlm) fit <- dynlm(y ~ xreg + L(y,1)) # lmtest::bptestreturns the Breusch and Pagan test statistic require(lmtest) bptest(fit)  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)

• 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. – Ray May 21 '15 at 18:07
• 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. – javlacalle May 22 '15 at 7:16
• 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). – javlacalle May 22 '15 at 7:17

i think this can help:

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)
#cat("S.sq=",S.sq,"\n")
.U <- (residuals^2)/S.sq
#cat(".U=",.U,"\n")
mod <- lm(.U ~ fitted.values(model))
varnames <- "fitted.values"
var.formula <- ~ fitted.values
df <- 1
SS <- anova(mod)\$"Sum Sq"
#cat("SS=",SS,"\n")
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"
#result
}

arima_model=auto.arima(time_series_data)
summary(arima_model)
ncvTest.arima(arima_model)