Unit root tests and stationarity Two common methods of testing whether a time series is stationary are the KPSS and ADF tests. If my understanding is correct, these tests essentially work by measuring the residuals of fitting the time-series to an autoregressive model which is linear.
So my question is this, if the time series is possibly of a non-linear nature are the results of the above tests still valid?
 A: The short answer is no at least for the ADF. I suspect a similar reasoning as the one outlined below applies for the KPSS, but I have not investigated this. The reason that the ADF will not work is that it is based on the notion of integration of order d, $d\geq 0$ to capture nonstationarity. (KPSS also is based on this notion, hence my suspicion that a similar problem as the one below will arise for KPSS.) (1) recaps the definition of integration order, while (2) explains why/how this definition makes no sense as a basis for a test in nonlinear settings using the DF test. 
The derivations here are only schematic and aimed at making conceivable the intuition of the (A)DF test family. I hope they will answer the question satisfactory.
(1) Nonstationarity: The ADF and the KPSS are defined for linear nonstationary processes. For such processes, the 'degree' of nonstationarity is traditionally captured by the order of integration. One says that a process $X_t \sim I(1)$ (integrated of order 1) $\Longleftrightarrow$ $\Delta X_t \sim I(0)$ (integrated of order 0, i.e. stationary). Similarly, $X_t \sim I(d)$ $\Longleftrightarrow$ $\Delta X_t \sim I(d-1)$ $\Longleftrightarrow$ $\Delta^d X_t \sim I(0)$. (Where $\Delta^d $ denotes that we have differenced $d$ times.)
For example, with $\varepsilon_i   \overset{iid}{\sim} N(0, \sigma^2), \; 1\leq i \leq N$, define the random variable $S_i \equiv \sum_{j=1}^i\varepsilon_j$ and define $R_i \equiv \sum_{j=1}^i S_j$ for $1\leq i \leq N$. Then clearly, $\varepsilon_i \sim I(0)$. Thus, by construction $\Delta S_i = S_i - S_{i-1} = \varepsilon_i \sim I(0)$ which implies that $S_i \sim I(1)$. Using similar reasoning, one can see that $R_i \sim I(2)$.
(2) Example (A)DF: Conceptually, there is no major difference between the ADF and the DF. Both are based on the convergence of the autoregressive coefficient to a functional of Brownian Motion. For mathematical simplicity, I will use a simple nonlinear process in the DF (rather than the ADF) framework to show why it is generally not feasible to apply the ADF test for nonlinear time series. I progress by first demonstrating what the DF test does in a linear setting and then give an example where it fails in a nonlinear setting.
(2a) Linear setting: let $\varepsilon_t   \overset{iid}{\sim} N(0, \sigma^2), \; - \infty \leq t \leq N$ and consider the process
\begin{align}
X_t &= X_{t-1} + \varepsilon_t =  \sum_{i=0}^{\infty} \varepsilon_{t-i}\\
\end{align}
Where the last stop follows because we can always express $X_t$ in deviations from the initial condition and so wlog impose $x_0 = 0$. Clearly, this process is $I(1)$. The DF test now works by estimating $\beta$ in the statistical model 
\begin{align}
X_t &= \beta X_{t-1} + \varepsilon_t
\end{align}
by using OLS. Suppose $N>1$ and we have observations $x_i:1\leq i \leq N$. Then the OLS estimator can simply be written as
\begin{align}
\hat{\beta} &= \frac{\sum_{i=1}^N x_t x_{t-1}} {\sum_{i=1}^Nx_{t-1}^2} \\
&= \frac{N^{-1}\sum_{i=1}^N (x_{t-1} + \varepsilon_t) x_{t-1}} {N^{-1}\sum_{i=1}^Nx_{t-1}^2} \\
&= 1 + \frac{N^{-1}\sum_{i=1}^N  \varepsilon_t x_{t-1}} {N^{-1}\sum_{i=1}^Nx_{t-1}^2} \\
\end{align}
Now note that the functional central limit theorem implies two convergence results (in distribution) for the fraction's numerator and the denominator:
\begin{align}
N^{-1}\sum_{t=1}^N  \varepsilon_t x_{t-1}  =  
\sum_{t=1}^N  \left( \frac{\varepsilon_t}{\sqrt{N}} \cdot \frac{\sum_{j=1}^{t-1}\varepsilon_{j-1}}{\sqrt{N}}\right) 
&\Longrightarrow \sigma^2\int_0^1W(r) dW(r) \\
N^{-2}\sum_{t=1}^Nx_{t-1}^2 = 
\sum_{t=1}^N  \left( \frac{\sum_{j=1}^{t}\varepsilon_{j-1}}{\sqrt{N}}
\right)^2
&\Longrightarrow \sigma^2\int_0^1W(r)^2 dr 
\end{align}
Where $W(r):r\in[0,1]$ denotes a Standard Brownian Motion/Wiener Process (i.e., $W(r) \sim N(0,r)$. When putting this all together, the continuous mapping theorem implies that
\begin{align}
N(\hat{\beta} -1) &\Longrightarrow \frac{\int_0^1W(r) dW(r)}{\int_0^1W(r)^2 dr} \equiv DF
\end{align}
which is known as the Dickey-Fuller (DF) distribution. The DF test now computes $N(\hat{\beta} -1)$ in the model above. Under the alternative of stationarity, $\hat{\beta}$ will consistently estimate $\beta$ and so the statistic will diverge to $\infty$. The critical values of the DF distribution are then used to reject/accept the null hypothesis of $X_t \sim I(1)$.
Clearly, this setting only works in the simplistic first order autoregression model. For more involved statistical models of the linear process model class, the ADF is used. The only further adjustment the ADF test makes is to restructure the regression equation such that the parameter $\phi = (1-\beta)$ is estimated directly, see e.g. Wikipedia (https://en.wikipedia.org/wiki/Augmented_Dickey%E2%80%93Fuller_test). Conceptually, this is the same as a normal DF test, but it converges to different functionals of Brownian motions in the limit (i.e., the appropriate critical values are different).
The take away for either the DF or the ADF is that what the test statistic really relies on is the meaningfulness of the first difference of $X_t$. The parameter $\beta$ (or, equivalently, the parameter $\phi = (1-\beta)$) converge in distribution to a functional of Brownian Motion ONLY if $X_t \sim I(1)$. This is the null hypothesis of the test, and if $X_t$ is not integrated of order one (but order 2, 0, or not integrated of any order while still nonstationary) the test statistics are not useful anymore.
Summary: The (A)DF tests against integration of order 1, NOT against general nonstationarity!
(2b) Nonlinear Setting: Suppose the data are generated by the following statistical model with nonlinear function $f$:
\begin{align}
X_t = f(X_{t-1})+\varepsilon_t
\end{align}
Using recursive substitution, one can observe that $X_t = f(f(...f(X_{t-r}) + \varepsilon_{t-r+1})+ \varepsilon_{t-r+2}) ... \varepsilon_{t-1}) + \varepsilon_t$. Applying the DF or ADF test to data that has a data generating process captured by the above equation amounts to fitting a regression to a non-linear process. Consequently, the estimated parameter $\hat{\beta}$ would not have the same interpretation/meaning as in the linear world*. It would also depend on the exact functional form of $f(\cdot)$ whether and where to the OLS estimate $\hat{\beta}$ would converge. Generally however, the (A)DF test would be invalid unless $f$ would be linear in $X_{t-1}, ..., X_{t-p}$ for some $p\geq 0$
*In particular, it would have the interpretation of the 'best linear projection' of $X_{t-1}, ..., X_{t-p}$ on $X_t$.
