High level overview of `auto.arima` with `xreg` predictors I'm trying to understand how auto.arima with covariates in the xreg parameter works. I'm familiar with regression and I'm starting to work on forecasting.
My understanding of forecasting is that you look for patterns in the past time series and then project those paterns onto the future.  
My uderstanding of regression is that you use predictors to try to generate an output value and minimize the difference between your created value and the real value.  
So how does forecasting auto.arima with xreg work? Do you create a forecast for a timeseries based on past data and regression model based on the input time series and input xreg, and then forecast each data point in the time series and for each forecasted data point use the regression model you built and future xreg values to adjust the forecasted values?
I'm a former physics grad student, so I'm not allergic to math but I'm just looking for a high level overview of the process here to understand how forecasting auto.arima works.  
For example like, 


*

*step 1: build forecast model on input time series, and regression model on input time series and input xreg values

*step 2: forecast model into future one step, and predict value with regression model and future xreg values

*step 3: algorithm combines forecasted value and regression model prediction to get combined value
This is just a guess at how it works, but it's an example of the kind of high level explanation I'm looking for.
I've included some code below that I've been working on trying to forecast time in to out TiTo for customers at a restaurant with predictor count of customers in the restaurant CustCount.
OV<-zoo(SampleData$TiTo, 
    order.by=SampleData$DateTime)


eDate <- ts(OV, frequency = 24)

Train <-eDate[1:15000]
Test <- eDate[15001:22773]

xregTrain <- SampleData[1:15000,]$CustCount
    xregTest <- SampleData[15001:22773,]$CustCount

Arima.fit <- auto.arima(Train, xreg = xregTrain)

Acast<-forecast(Arima.fit, h=7772, xreg = xregTest)

accuracy(Acast$mean,Test)

 A: I think the real question is not about auto.arima but rather about regression with ARMA errors as implemented in arima with exogenous regressors. 

So how does forecasting auto.arima with xreg work? Do you create a forecast for a time series based on past data and regression model based on the input time series and input xreg, and then forecast each data point in the time series and for each forecasted data point use the regression model you built and future xreg values to adjust the forecasted values?

No, this is not how it works. It is simpler than that. Let us call the time series of interest $y$ and the exogenous regressors $X$. In a regression with ARMA errors as implemented in arima (and auto.arima) with xreg, the model equations are
$$
\begin{align}
y_t &= \beta' X_t + u_t \\
u_t &= \varphi_1 u_{t-1} + \dotsc + \varphi_p u_{t-p} + \varepsilon_t + \theta_1\varepsilon_{t-1} + \dotsc + \theta_q\varepsilon_{t-q}
\end{align}
$$
where $u_t$ is not $i.i.d.$ but $\varepsilon_t$ is.
The "main" equation is the upper one which specifies a linear relationship between $y_t$ and $X_t$ and looks very much like a regular regression model. The difference from a regular regression is made by the second equation that specifies the error term $u_t$ to follow an ARMA($p,q$) process (instead of $u_t$ being uncorrelated as would be the case in a regular regression). The second equation serves to "fix the problem with the first equation's errors". I like to think of this model as basically the good old regression but with unpleasant error pattern that needs to be accounted for.
How does auto.arima differ from arima in this case? auto.arima selects the optimal autoregressive and moving-average orders $p$ and $q$ based on a chosen information criterion (AICc by default, alternatively AIC or BIC) from a local search over a few regions of values. Meanwhile, arima requires the user to specify $p$ and $q$ manually.
How are these models fit? Generally, by maximum likelihood estimation. Refering to your three-step example algorithm, arima can be done in one step where all the parameters ($\beta$s, $\varphi$s and $\theta$s) are estimated at once together producing the fitted values $\hat y_t$ and the residuals $\hat\varepsilon_t$. auto.arima would do that multiple times as it tries out different values of $p$ and $q$.
See more in the Hyndman & Khandakar "Automatic time series forecasting: the forecast package for R" (2008). Also, try Rob J. Hyndman's blog post "The ARIMAX model muddle".
