The Arima function in the forecast package can fit a regression model to the data with an ARIMA model for the errors. The order argument specifies the orders of the ARIMA model, while the argument xreg defines which data object contains the observations of the predictors. E.g., if xreg is a matrix of predictors:
model = Arima(series, order = c(1,1,0), xreg = covariates)
To find the order of the ARIMA process, you can simply use the auto.arima function also found in the forecast package. It automatically locates the best-fitting ARIMA model to the data, “fit” defined by one of three possible information criteria in the ic argument: the AIC (given by aic), the AICc (aicc), or the BIC (bic). E.g.,
model = auto.arima(series, ic = “aic”)
I think you may find this page really helpful, especially the section about R.