The explosive AR(1) process with $\varphi>1$, where was this first represented as a stationary, but non-causal, time-series? According to this question and answer Explosive AR(MA) processes are stationary? the AR(1) process (with $e_t$ white noise):
$$X_{t}=\varphi X_{t-1}+e_{t} \qquad , e_t \sim WN(0,\sigma)$$
is a stationary process if $\varphi>1$ because it can be rewritten as
$$X_t=\sum_{k=0}^\infty {\varphi}^{-k}u_{t+k}$$
But now the variable $X_t$ depends on the future.

I wonder where this representation (which I remember having seen in several places) and the derivation originally comes from.

I am confused about the derivation, and I wonder how it works. When I try to do the derivation myself I am failing.
I can rewrite the process $$X_{t+1}=\varphi X_{t}+e_{t+1}$$ as
$$X_{t}= \varphi^{-1} X_{t+1} -\varphi^{-1}  e_{t+1}$$
and replacing $\varphi^{-1} e_{t+1}$ by $u_{t}$ it becomes
$$X_{t}= \varphi^{-1} X_{t+1} + u_{t}$$
such that the expression is 'like' another AR(1) process but in reverse time and now the coefficient is below 1 so it seemingly is stationary (*).
From the above it would follow indeed that $$X_t=\sum_{k=0}^\infty {\varphi}^{-k}u_{t+k}$$
(*) But the $u_t$ is not independent from $X_{t+1}$, because it is actually $e_{t+1}$ times a negative constant.
 A: The question suggests some basic confusion between the equation and the solution
The Equation
Let ${\varphi} > 1$.
Consider the following (infinite) system of equations---one equation for each $t\in \mathbb{Z}$:
$$
X_{t}=\varphi X_{t-1}+e_{t}, \mbox{ where } e_t \sim WN(0,\sigma), \;\; t \in \mathbb{Z}. \quad (*)
$$
Definition
Given $e_t \sim WN(0,\sigma)$, a sequence of random variables $\{ X_t \}_{t\in \mathbb{Z}}$ is said to be a solution of $(*)$ if, for each $t$,
$$
X_{t}=\varphi X_{t-1}+e_{t},
$$
with probability 1.
The Solution
Define
$$
X_t= - \sum_{k=1}^\infty {\varphi}^{-k}e_{t+k},
$$
for each $t$.

*

*$X_t$ is well-defined: The sequence of partial sums
$$
X_{t,m} = - \sum_{k=1}^m {\varphi}^{-k}e_{t+k}, \;\; m \geq 1
$$
is a Cauchy sequence in the Hilbert space $L^2$, and therefore converges in $L^2$. $L^2$ convergence implies convergence in probability (although not necessarily almost surely). By definition, for each $t$, $X_t$ is the $L^2$/probability-limit of $(X_{t,m})$ as $m \rightarrow \infty$.


*$\{ X_t \}$ is, trivially, weakly stationary. (Any MA$(\infty)$ series with absolutely summable coefficients is weakly stationary.)


*$\{ X_t \}_{t\in \mathbb{Z}}$ is a solution of $(*)$, as can be verified directly by substitution into $(*)$.
This is a special case of how one would obtain a solution to an ARMA model:
first guess/derive an MA$(\infty)$ expression, show that it is well-defined, then verify it's an actual solution.
$\;$

...But the $\epsilon_t$ is not independent from $X_{t}$...

This impression perhaps results from confusing the equation and the solution.
Consider the actual solution:
$$
\varphi X_{t-1} + e_t = \varphi \cdot \left( - \sum_{k=1}^\infty {\varphi}^{-k}e_{t+k-1} 
\right) + e_t,
$$
the right-hand side is exactly $- \sum_{k=1}^\infty {\varphi}^{-k}e_{t+k}$, which is $X_t$ (we just verified Point #3 above). Notice how $e_t$ cancels and actually doesn't show up in $X_t$.
$\;$

...where this...derivation originally comes from...

I believe Mann and Wald (1943) already considered non-causal AR(1) case, among other examples. Perhaps one can find references even earlier. Certainly by the time of Box and Jenkins this is well-known.
Further Comment
The non-causal solution is typically excluded from the stationary AR(1) model because:

*

*It is un-physical.


*Assume that $(e_t)$ is, say, Gaussian white noise. Then, for every non-causal solution, there exists a causal solution that is observationally equivalent, i.e. the two solutions would be equal as probability measures on $\mathbb{R}^{\mathbb{Z}}$. In other words, a stationary AR(1) model that includes both causal and non-causal cases is un-indentified. Even if the non-causal solution is physical, one cannot distinguish it from a causal counterpart from data. For example, if innovation variance $\sigma^2 =1$, then the causal counterpart is causal solution to AR(1) equation with coefficient $\frac{1}{\varphi}$ and  $\sigma^2 =\frac{1}{\varphi^2}$.
A: Re-arranging your first equation and increasing the index by one gives the "reverse" AR(1) form:
$$X_{t} = \frac{1}{\varphi} X_{t+1} - \frac{e_{t+1}}{\varphi}.$$
Suppose you now define the observable values using the filter:
$$X_t = - \sum_{k=1}^\infty \frac{e_{t+k}}{\varphi^k}.$$
You can confirm by substitution that both the original AR(1) form and the reversed form hold in this case.  As pointed out in the excellent answer by Michael, this means that the model is not identified unless we exclude this solution by definition.
A: 
... the AR(1) process (with $e_t$ white noise):
$$X_{t}=\varphi X_{t-1}+e_{t} \qquad , e_t \sim WN(0,\sigma)$$
is a stationary process if $\varphi>1$ because ...

It seems me not possible as showed there: https://en.wikipedia.org/wiki/Autoregressive_model#Example:_An_AR(1)_process
for wide sense stationarity $-1 < \varphi < 1$ must hold.
Moreover, maybe I lose something here but it seems me that not only the process above cannot be stationary but it is entirely impossible and/or bad defined.
This because if we have an autoregressive process, we do not stay in a situation like $Y=\theta Z+u$ where $Z$ and $u$ can be two unrestricted random variables and $\theta$ an unrestricted parameter.
In a regression residuals and parameters are not free terms, given dependent and independent/s variables, they are given too.
So, in AR(1) case it is possible to show that $-1 \leq \varphi \leq 1$ must hold; like autocorrelation.
Moreover if we assume that $e_t$ (residuals) are white noise process ... we make a restriction on $X_t$ process too. If in the data we estimate an AR(1) and $e_t$ result as autocorrelated ... the assumption/restriction do not hold ... AR(1) is not a good specification.
