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