Function $g(x)$, such that $\int_{\mathbb{R}} f(g(x))dx = \int_{\mathbb{R}} f(x)dx$ for all $f\in L^1(\mathbb{R})$ It's not too hard to show that for all $f\in L^1(\mathbb{R})$ and for $g(x) = x-\frac{1}{x}$ we have $\int_{\mathbb{R}} f(g(x))dx = \int_{\mathbb{R}} f(x)dx$:use the substitution $y = x - \frac{1}{x}$ and split the integral in two parts: $x<0$ and $x>0$. On one of them make the choice $x = \frac{y + \sqrt{y^2 + 4}}{2}$, on the other $x = \frac{y - \sqrt{y^2 + 4}}{2}$ and add the resulting integrals. 
The existence of even one function $g(x)$, such that the above works is astonishing. Are there other choices of $g$, such that $\int_{\mathbb{R}} f(g(x))dx = \int_{\mathbb{R}} f(x)dx$ for all $f\in L^1(\mathbb{R})$?
Edit: $g(x) = x + a$ for any a would be a solution, but I was wondering if there are more complicated $g$ which work, or if $x - \frac{1}{x}$ is just a single example.
 A: The set of all such functions is closed under composition. As a result, all the functions $g^k$ (in the sense of composition) work.
E.g. $$g^2(x)=x-\frac1x -\frac{1}{x-\frac1x}$$
Adding to that the fact that you can also compose them with affine functions of slope $1$ (lisyarus' answer) yields already a pretty large class.

Incorporating @J.G.'s comment, one gets for instance all rational functions of the form $$\frac{x^2+ax+b}{x+c}$$ with the only requirement that $b<c^2$ ($a$ and $c$ are arbitrary).
A: We have $\int f\circ g=\int f$ for every integrable $f$ if and only if the same holds for $f=\chi_E$, which says precisely $$m(g^{-1}(E))=m(E)$$for every measurable set $E$, which is to say that $g$ is measure-preserving.
If we assume in addition that $g$ is a smooth bijection this is equivalent to $|g'|=1,$so the only smooth bijections with this property are $$g(x)=\pm x+c.$$
This seemed so obvious that it seemed clear to me at first that $x-1/x$ cannot have the stated property. But of course that function is not injective; if $g(x)=x-1/x$ you  can calculate $g^{-1}((a,b))$ explicitly, and sure enough you get two intervals with total length $b-a$.
(This is actually consistent  with the analysis above, if you look at it right. The relevant condition for a smooth bijection is actually $|({g^{-1})}'|=1$,  which is equivalent to $|g'|=1$. But now if $g(x)=x-1/x$ and you let $y_1$ and $y_2$ denote the two "branches" of $g^{-1}$ you easily calculate that $|y_1'|+|y_2'|=1$.)
One can easily concoct discontinuous examples. For example if $A\subset(0,\infty)$ let $$g(x)=\begin{cases}x,&(|x|\in A),
\\-x,&(|x|\notin A).\end{cases};$$then $g$ is a measure-preserving bijection.
A: There is a not that well known result about these kinds on integrals that quite often comes up on this site. It was already mentioned in the comments by achille hui, but I'll add this as a CW as it's a shame not having it as an answer:
Glasser's master theorem$^{[1]}$ says that if $f(x)$ is integrable then
$$\int_{\mathbb{R}}f(x)\,{\rm d}x = \int_{\mathbb{R}}f(\phi(x))\,{\rm d}x$$
for all $\phi(x) = x - \sum_{n=1}^N\frac{|a_n|}{x-b_n}$ where $a_n,b_n$ are arbitrary constants and where the integrals are to be interpreted in a principal value sense.
[1]: Glasser, M. L. "A Remarkable Property of Definite Integrals." Math. Comput. 40, 561-563, 1983
A: Clearly if $g: \mathbb{R}\to\mathbb{R}$ is a monotonic, we must have have, under $u=g(x)$,
$$ \int_{\mathbb{R}}f(g(x))dx=\int_{\mathbb{R}}f(u)(g^{-1})'(u)du $$
for all $f\in L^1(\mathbb{R})$. Thus
$$(g^{-1})'(u)=1$$
which implies $g(x)=x+a$. If $g: (-\infty,0)\to\mathbb{R}$ and $g: (0,\infty)\to\mathbb{R}$ are monotonic, namely $g^{-1}$ has two branches $x=g_1^{-1}(u):\mathbb{R}\to(-\infty,0)$ and $x=g_2^{-1}(u):\mathbb{R}\to(0,\infty)$, then
\begin{eqnarray}
\int_{\mathbb{R}}f(g(x))dx&=&\int_{(-\infty,0)}f(g(x))dx+\int_{(0,\infty)}f(g(x))dx\\
&=&\int_{\mathbb{R}}f(u)(g_1^{-1})'(u)du+\int_{\mathbb{R}}f(u)(g_2^{-1})'(u)du\\
&=&\int_{\mathbb{R}}f(u)\bigg[(g_1^{-1})'(u)+(g_2^{-1})'(u)\bigg]du
\end{eqnarray}
holds for all $f\in L^{-1}(\mathbb{R})$. Thus we must have
$$ (g_1^{-1})'(u)+(g_2^{-1})'(u)=1 $$
or
$$ g_1^{-1}(u)+g_2^{-1}(u)=u. \tag{1} $$
or all $u\in\mathbb{R}$. Assume $g(x)$ has the form $g(x)=ax+\frac{b}{x+c}$. Thus if $g(x)$ is monotonic, then $a$ and $b$ must satisfy $ab<0$. WOLG, assume $a>0,b<0$. Using (1), it is easy to obtain $a=1$. So
$$ g(x)=x+\frac{b}{x+c}, b<0 $$
A: As mentioned by David C. Ullrich, such $g$ is called a measure-preserving transformation (MPT) on $\mathbb{R}$ (w.r.t. the Lebesgue measure, of course). The comment by achille hui provides an important family of MTPs given by
$$ g(x) = \pm \left( x - c - \sum_{k=1}^{n} \frac{\mu_k}{x - \lambda_k} \right) $$
for $c, \lambda_k \in \mathbb{R}$ and $\mu_k \in (0, \infty)$, which is the statement of the Glasser's master theorem (1983). More generally, a theorem by Letac (1977) tells that any function of the form
$$ g(x) = \pm \left( x - c - \int_{\mathbb{R}} \left( \frac{1}{x - \lambda} + \frac{\lambda}{1+\lambda^2} \right) \, \mu(d\lambda) \right)$$
for $c \in \mathbb{R}$ and a singular Borel measure $\mu$ with $\int_{\mathbb{R}} \frac{\mu(d\lambda)}{1+\lambda^2} < \infty$ is MPT. This theorem includes interesting examples such as
$$g(x) = x - \cot x$$
corresponding to $\mu = \sum_{k \in \mathbb{Z}} \delta_{\pi k}$. Notice also that the Glasser's master theorem is a special case of this theorem with $\mu = \sum_{k=1}^{n} \mu_k \delta_{\lambda_k}$, although the proof techniques involved are different.
A: Consider $g(x) = x+a$ for a fixed $a\in \mathbb R$. Then, using the substitution $y = x+a$, we get $dy = dx$, and thus
$$\int\limits_{\mathbb R}f(g(x))dx = \int\limits_{\mathbb R}f(y)dy$$
A: The functions of the form 
$$\phi(x) = x - \sum_{i=1}^{n-1} \frac{\rho_i}{x- \alpha_i} -\beta $$
where $\rho_i>0$ for all $1\le i \le {n-1}$ invariate the Lebesgue measure.
The set of such functions is closed under composition. 
Sketch of proof:
Assume $\alpha_1 < \ldots <\alpha_{n-1}$. On each of the intervals $(-\infty, \alpha_1)$, $(\alpha_1, \alpha_2)$, $\ldots$, $(\alpha_{n-2}, \alpha_{n-1})$,$(\alpha_{n-1},\infty)$ the function $\phi$ is strictly increasing from $-\infty$ to $\infty$. Therefore, for every $u \in \mathbb{R}$ the equation
$$\phi(x)=u$$ has exactly $n$ distinct real roots. From Viete's relations we see that the sum of the roots equals 
$$u +\sum_{i=1}^{n-1}\alpha_i + \beta$$
Therefore, for every $I$ interval of $\mathbb{R}$ the set $\phi^{-1}(I)$ is a disjoint union of interval of length $|I|$. This implies
$$\int_{\mathbb{R}}(\chi_I\circ \phi)\  d\mu= \int_{\mathbb{R}}\chi_I \ d \mu$$
for the characteristic function of the interval $I$. From here one concludes that 
$$\int f \circ\phi \ d\mu = \int f \ d\mu$$ for all $f \in L^1(\mathbb{R})$. 
Consider two such functions $\phi= x-\sum_1^{n-1} \frac{\rho_i}{x-\alpha_i} - \beta$ and 
$\chi = x - \sum_1^{n'-1}\frac{\rho'_i}{x-\alpha'_i} - \beta'$. 
Each of the intervals $(\infty, \alpha'_1)$, $(\alpha'_1,\alpha'_2)$, $\ldots$,$(\alpha_{n'},\infty)$ get divided into intervals that map to one of $(-\infty, \alpha_1)$, $(\alpha_1, \alpha_2)$, $\ldots$, $(\alpha_{n-2}, \alpha_{n-1})$,$(\alpha_{n-1},\infty)$.
 We conclude that $\mathbb{R}$ gets divided into $n\cdot n'$ intervals that are mapped strictly increasing onto $(-\infty, \infty)$ by $\phi\circ \psi$ . 
The rational function $\phi\circ \psi$ has therefore at least $n\cdot n'-1$ real finite poles and a pole at infinity. 
Since it is of degree $n\cdot n'$ it must have a partial fraction  decomposition of form 
$$x - \sum_1^{n n'-1} \frac{\rho''_i}{x- \alpha''_i}- \beta''$$
Now one checks that the $\rho''_i$ must be positive. 
