Breusch-Pagan Test for ARIMA Model in R 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")])

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

A: 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)
    #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)

