model <- Arima(y , xreg=cbind(x1, x2), order=(p,d,q))
If I am forecasting $x_1$ and $x_2$, then for forecasting $y$:
1) If I use expected forecasts for $x_1$ and $x_2$ (single numbers), I simply do:
forecast(model, xreg=cbind(E(future x1) , E(future x2))
2) How about if I want to use forecast intervals for $x_1$ and $x_2$?
This post suggests that: you can draw (a lot of) random numbers from each predictive density, plug them into the model and get a predictive distribution for $y$. Then I guess taking the average prediction interval to come up with the one forecast interval for $y$. Does this make sense?