Yes. Here is an example. The true data is simulated as

data t1;
  err0=0;
  do i=-100 to 1000;
    x=rannor(123);
        err=0.7*err0+rannor(123);
    y=2+0.8*x+err;
        if i>=1 then output;
        err0=err;
   end;
   keep x y;
run;

The intercept is 2 and coef. of x is 0.8 and ar of error is 0.7. And z is a noise variable added to illustrate it supported multiple x-variables.  

simudata<- read.csv("c:\\temp\\test.csv");
y<- simudata$y
x<- simudata$x
z<- rnorm(length(simudata$x))
xreg=cbind(x,z)

arima(y, order = c(1,0,0), xreg = x)
arima(y, order = c(1,0,0), xreg = xreg)

*********************
Call:
arima(x = y, order = c(1, 0, 0), xreg = x)

Coefficients:
         ar1  intercept       x
      0.6592     1.8110  0.7897
s.e.  0.0237     0.0918  0.0278

sigma^2 estimated as 0.9822:  log likelihood = -1410.27,  aic = 2828.54
> arima(y, order = c(1,0,0), xreg = xreg)

Call:
arima(x = y, order = c(1, 0, 0), xreg = xreg)

Coefficients:
         ar1  intercept       x        z
      0.6595     1.8098  0.7895  -0.0174
s.e.  0.0237     0.0919  0.0278   0.0268

sigma^2 estimated as 0.9818:  log likelihood = -1410.06,  aic = 2830.12

多谢~~