Arima time series forecast (auto.arima) with multiple exogeneous variables in R

I would like to conduct a forecast based on a multiple time series ARIMA-model with multiple exogeneous variables. Since I am not that skillfull with regards to neither statistics nor R I want to keep is as simple as possible (Trend forecast for 3 months is sufficient).

I have 1 dependent time series and 3-5 predictor time series, all monthly data, no gaps, same time “horizon”.

I encountered the auto.arima function and asked myself if this would be a suitable solution for my problem. I have different commodity prices and prices of products made from them. All raw-data are non-stationary but via first-order differencing they all become stationary data. ADF, KPSS indicate this. (This means that I have tested for integration, right?).

My question now is: How do I apply this with the auto.arima function AND is ARIMA the right approach anyways? Some ppl already adviced me to use VAR, but is it possible with ARIMA too?

The following table is my data. Actually the data-set goes up til 105 observations, but the first 50 will do. Trend as well as seasonality are obviously of interest here.

Thanks for any advices and help!
Georg

If your external regressors are causal for $y$, but not the other way around and do not cause each other, then ARIMA is definitely appropriate. VAR makes sense if your different time series all depend on each other.
For auto.arima() to work with external regressors, collect your regressors into a matrix X, which you feed into the xreg parameter of auto.arima(). (Of course, X must have the same number of rows as the time series y you are modeling.)
For forecasting, you will need the future values of your regressors, which you then again feed into the xreg parameter of forecast.
The help pages are ?auto.arima and ?forecast.Arima (note the capital A – this is not a typo. Don’t ask me…).