Can a trend stationary series be modeled with ARIMA? I have a question / confusion about stationary series required for modeling with ARIMA(X). I am thinking of this more in terms of inference (effect of an intervention), but would like to know if forecasting versus inference makes any difference in the response.
Question:
All the introductory resources I have read state that the series needs to be stationary, which makes sense to me and that is where the "I" in arima comes in (differencing). 
What confuses me is the use of trends and drifts in ARIMA(X) and implications (if any) for stationary requirements. 
Does the use of either a constant/drift term and/or trend variable as an exogenous variable (i.e. adding 't' as a regressor) negate the requirement of the series being stationary? Is the answer different depending on if the series has a unit root (e.g. adf test) or has a deterministic trend but no unit root? 
OR
Does a series always have to be stationary, made so via differencing and/or detrending before using ARIMA(X)?
 A: As empirically shown below, using a trend variable in ARIMAX  negated the need for differencing and makes the series trend stationary. Here is the logic I used to verify.

*

*Simulated an AR process

*Added a deterministic trend

*Using ARIMAX modeled with trend as exogenous variable the above series without differencing.

*Checked the residuals for white noise and it is purely random

Below is the R code and plots:
    set.seed(3215)
    
    ##Simulate an AR process
    x <- arima.sim(n = 63, list(ar = c(0.7)));
    plot(x)

    ## Add Deterministic Trend to AR
    t <- seq(1, 63)
    beta <- 0.8
    t_beta <- ts(t*beta, frequency=1)
    ar_det <- x+t_beta
    plot(ar_det)
    
    ## Check with arima
    
    ar_model <- arima(ar_det, order=c(1, 0, 0), xreg=t, include.mean=FALSE)
    
    ## Check whether residuals of fitted model is random

    pacf(ar_model$residuals)

AR(1) Simulated Plot

AR(1) with deterministic trend

ARIMAX Residual PACF with trend as exogenous. Residuals are random, with no pattern left

As can be seen above, modeling deterministic trend as an exogenous variable in the ARIMAX model negates the need for differencing. At least in the deterministic case it worked. I wonder how this would behave with stochastic trend which is very hard to predict or model.
To answer your second question, YES all ARIMA including ARIMAX have to be made stationary. At least that's what text books say.
In addition, as commented, see this article. Very clear explanation on Deterministic Trend vs. Stochastic trend and how to remove them to make it trend stationary and also very nice literature survey on this topic. They use it in the neural network context, but it is useful for general time series problem. Their final recommendation is when it  is clearly identified as deterministic trend, the do linear detrending, else apply differencing to make the time series stationary. The jury is still out there, but most researchers cited in this article recommend differencing  as opposed to linear detrending.
Edit:
Below is random walk with drift stochastic process, using exogenous variable and difference arima. Both appear to give same answer and in essence they are same.
    library(Hmisc)
    
    set.seed(3215)
    
    ## ADD Stochastic Trend to simulated Arima this is AR(1) with unit root with non zero mean
    
    y = rep(NA, 63)
    y[[1]] <- 2
        
    for (i in 2:63)  {
    y[i] <- 3+1*y[i-1] + rnorm(1, mean = 0, sd = 1)
    } 
    
    plot(y, type="l")
    
    y_ts <- ts(y, frequency=1)
    
    ## Lag to create Xreg
    
    y_1 <- Lag(y, shift=1)
    
    ## Start from 2 value to avoid NA and make it equal length with xreg
    
    y <- window(y_ts, start =2, end=63)
    xreg1 <- y_1[-1]
    
    ## Check the values with ARIMA and xreg
    
    g <- arima(y, order=c(0, 0, 0), xreg=xreg1)
    
    pacf(g$residuals)
    
    ## Check the values with ARIM
    
    g1 <- arima(y, order=c(0, 1, 0))
    
    pacf(g1$residuals)
    
    ## 

    ARIMA(0, 0, 0) with non-zero mean 
    
    Coefficients:
          intercept   xreg1
             3.1304  0.9976
    s.e.     0.2664  0.0025

A: Remember that there are different kinds of non-stationarity and different ways on how to deal with them. Four common ones are:
1) Deterministic trends or trend stationarity. If your series is of this kind de-trend it or include a time trend in the regression/model. You might want to check out the Frisch–Waugh–Lovell theorem on this one.
2) Level shifts and structural breaks. If this is the case you should include a dummy variable for each break or if your sample is long enough model each regimé separately.
3) Changing variance. Either model the samples separately or model the changing variance using the ARCH or GARCH modelling class.
4) If your series contain a unit root. In general you should then check for cointegrating relationships between the variables but since you are concerned with univariate forecasting you shoud difference it once or twice depending on the order of integration.
In order to model a time series using the ARIMA modelling class the following steps should be appropriate:
1) Look at the ACF and PACF together with a time series plot to see wheter or not the series is stationary or non-stationary. 
2) Test the series for a unit root. This can be done with a wide range of tests, some of the most common being the ADF test, the Phillips-Perron (PP) test, the KPSS test which has the null of stationarity or the DF-GLS test which is the most efficient of the aforementioned tests. NOTE! That in case your series contain a structural break these tests are biased towards not rejecting the null of a unit root. In case you want to test the robustness of these tests and if you suspect one or more structural breaks you should use endogenous structural break tests. Two common ones are the Zivot-Andrews test which allows for one endogenous structural break and the Clemente-Montañés-Reyes which allows for two structural breaks. The latter allows for two different models. An additive outlier model which accounts for sudden changes in the slope of the series and an innovative outlier model which takes gradual changes into account and allows a break in the intercept and slope.
3) If there is a unit root in the series then you should difference the series. Afterwards you should run look at the ACF, PACF and the time series plot and probably check for a second unit root to be on the safe side. The ACF and PACF will help you decide on how many AR and MA terms you should be including. 
4) If the series does not contain a unit root but the time series plot and the ACF show that the series has a deterministic trend you should add a trend when fitting the model. Some people argue that it is completely valid to just difference the series when it contains a deterministic trend although information may be lost in the process. Never the less its a good idea to difference it in order to see have many AR and/or MA terms you will need to include. But a time trend is valid.
5) Fit the different models and do the usual diagnostic checking, you might want to use an information criterion or the MSE in order to select the best model given the sample you fit it on.
6) Do in sample forecasting on the best fitted models and calculate loss functions such as MSE, MAPE, MAD to see which of them actually perform best when using them to forecast because that is what we want to do!
7) Do your out of sample forecasting like a boss and be pleased with your results!
A: Determining whether the trend (or other component such as seasonality) is deterministic or stochastic is part of the puzzle in time series analysis. I will add a couple of points to what has been said.
1) The distinction between deterministic and stochastic trendsis important because if a unit root is present in the data (e.g. a random walk) then test statistics used for inference do not follow the traditional distribution. See this post for some details and references.
We can simulate a random walk (stochastic trend where first differences should be taken), test for the significance of deterministic trend and see the percentage of cases in which the null of deterministic trend is rejected. In R, we can do:
require(lmtest)
iter <- 10000
cval <- 0.05
n <- 120
rejections <- 0
set.seed(123)
for (i in seq.int(iter))
{
  x <- cumsum(rnorm(n)) # random walk
  fit <- lm(x ~ seq(n))
  if (coeftest(fit)[2,"Pr(>|t|)"] < cval)
    rejections <- rejections + 1
}
100 * rejections / iter
#[1] 88.67

At the 5% significance level, we would expect to reject the null in the 95% of cases, however, in this experiment it was rejected only in ~89% of cases out of 10,000 simulated random walks.
We can apply unit root tests to test whether a unit root is present. But we must be aware that a linear trend may in turn lead to failure to reject the null of a unit root. In order to deal with this, the KPSS test considers the null of stationarity around a linear trend.
2) Another issue is the interpretation of the deterministic components in a process in levels or first differences. The effect of an intercept is not the same in a model with a linear trend as in a random walk. See this post for illustration.
Analytically, let's take a random walk with drift:
$$
y_t = \mu + y_{t-1} + \epsilon_t \,,\quad \epsilon_t \sim NID(0, \sigma^2) \,.
$$
If we substitute repeatedly $y_{t-i}$ by lagged versions of $y_t$:
\begin{eqnarray*}
y_t &=& \mu + \underbrace{y_{t-1}}_{\mu + y_{t-2} + \epsilon_{t-1}} + \epsilon_t \\ 
&=& 2\mu + \underbrace{y_{t-2}}_{\mu + y_{t-3} + \epsilon_{t-2}} + \epsilon_{t-1} + \epsilon_t \\
&=& 3\mu + y_{t-3} + \epsilon_{t-2} + \epsilon_{t-1} + \epsilon_t \\
&...&
\end{eqnarray*}
We arrive to:
$$
y_t = y_0 + \mu t + \sum_{i=1}^t \epsilon_i
$$
where $y_0$ is some arbitrary initial value. Thus, we see that the accumulation of shocks and the long memory of the random walk makes the intercept $\mu$ to have the effect of a linear trend with slope $\mu$ (in this case the constant term $\mu$ is called a drift).
If the graphical representation of a series shows a relatively clear linear trend, we cannot be sure whether it is due to the presence of a deterministic linear trend or to a drift in a random walk process. Complementary graphics and tests statistics should be applied.
There are some caveats to bear in mind since an analysis based on unit root 
and other test statistics is not foolproof. Some of these tests may be affected by the presence of outlying observations or level shifts and require the selection of a lag order which is not always straightforward.
As a workaround to this puzzle, I think that the common practice is to take differences of the data until the series looks stationary (for example looking at the autocorrelation function, which should go to zero fast) and then choose an ARMA model. 
A: Looking at the comments it seems that we didn't address the question about how to choose between a deterministic or stochastic trend. That is, how to proceed in practice rather than the consequences or properties of each case.
One way to proceed is the following: Start by applying the ADF test.


*

*If the null of a unit root is rejected we are done. The trend (if any) can be represented by a deterministic linear trend.

*If the null of the ADF test is not rejected then we apply the KPSS test 
(where the null hypothesis is the opposite, stationarity or stationarity around a linear trend).
o If the null of the KPSS test is rejected then we conclude that there is a unit root and work with the first differences of the data. Upon the first differences of the series we can test the significance of other regressors 
or choose an ARMA model.
o If the null of the KPSS test is not rejected then we would have to say that the data are not much informative because we weren't able to reject none the of the null hypotheses. In this case it may be safer to work with the first differences of the series.
As mentioned in a previous answer, remember that these tests may be affected by the presence of outliers (e.g. an outlier at a single time point due to an error when recording the data or a level shift due for example to a policy change that affects the series from a given time point on). Thus, it is advisable to check these issues as well and repeat the previous analysis after including regressors for some potential outliers.
