I think it would be worth exploring exponential smoothing models as well. Exponential smoothing models are a fundamentally different class of models from ARIMA models, and may yield different results on your data.
This sounds like a valid approach, and is very similar to the time series cross-validation method proposed by Rob Hyndman.
I would aggregate the cross-validation error from each forecast (exponential smoothing, ARIMA, ARMAX) and then use the overall error to compare the 3 methods.
You may also want to consider a "grid search" for ARIMA parameters, rather than using auto.arima. In a grid search, you would explore each possible parameter for an arima model, and then select the "best" ones using forecast accuracy.
To become better at time-series forecasting, it is no doubt beneficial to expand the number of forecasting methods or models available at your fingertips. Being able to model time-series data using ARIMA and exponential smoothing models is a good notch to have under your belt. Delving into non-linear models, regime switching models, and varying parameter models can only be a good thing for you.
The important thing to keep in mind is that we'd normally like to build simple linear models and not necessarily complicate matters by building non-linear models. This is one thing that you should definitely consider. That is, before considering building non-linear models, first test to see if non-linearity is actually present in the data under consideration.
You mention that you need to take into account weekends and public holidays. This suggests that maybe you need to make sure that the model you choose contains the appropriate seasonal (and non-seasonal, of course) components. Have you tried a seasonal non-seasonal multiplicate arima model? Or, a seasonal non-seasonal additive arima model?
I can only speak in general terms since you haven't given much details about the data you're working with, so a general word of advice would be to let the data do the talking and rather than choose a forecasting model a priori, let the data lead you to the appropriate model. If the data suggests that a non-linear model is worth considering, for sure, consider it.
In addition to F. Tusell's suggestion, another R package worth checking out is tsDyn. It is available at http://cran.r-project.org/web/packages/tsDyn/index.html.
Best Answer
You need the future values of the covariate to make ARIMAX (or perhaps regression with ARIMA errors – see The ARIMAX model muddle by Rob J Hyndman) feasible. If you do not have these values, you may need to forecast them. This could be done separately or jointly with the dependent variable. In the latter case, a multivariate time series model such as VAR (vector autoregression) could be used.