A functional equation involving the inverse function $\newcommand\ep\epsilon\newcommand\R{\mathbb R}$Let $P$ denote the set of all continuous probability density functions (pdf's) $p$ on $\R$ vanishing at $\pm\infty$. Let us say that a pdf $p\in P$ is good if for each small enough $\ep>0$ the functional equation
$$g(x)-g^{-1}(x)=\ep\, p(x)\quad\forall x\in\R\tag1$$
has a solution $g\colon\R\to\R$, which is an increasing continuous function such that $g(x)\ge x$ for all real $x$. Here $g^{-1}$ is of course the compositional inverse of $g$.
It is clear that, if a pdf $p\in P$ is good, then for any real $a$ and any real $b>0$ the pdf $p_{a,b}$ given by the formula $p_{a,b}(x):=b\,p(a+bx)$ for real $x$ is good as well.

The problem here is to characterize the set of all good pdf's $p\in P$.

Of course, there is always a tautological characterization: a pdf $p\in P$ is good if and only if it is good. Any non-tautological characterization would be of interest, including incomplete ones, such as conditions that are only sufficient or only necessary for the goodness. In particular, it would be of interest to know if the "triangular" pdf $p_\triangle$ given by the formula $p_\triangle(x):=\max(0,1-|x|)$ for real $x$ is good.
This question is related to this answer.

Comment 1: Geometrically, the left-hand side $g(x)-g^{-1}(x)$ of the functional equation (1) is the length $l(x)$ of the cross-section by the vertical line with abscissa $x$ of the region
$$A_g:=\{(u,v)\in\R^2\colon v\le g(u)\ \&\ u\le g(v)\},$$
symmetric about the diagonal $D:=\{(u,v)\in\R^2\colon u=v\}$.
Here is a picture of $A_g$ for $g(x)\equiv x+\frac{1}{2 \left(x^2+1\right)}$:

We want to choose $g$ so that this vertical length $l(x)$ be proportional to $p(x)$.
Comment 2: One may want to try to solve equation (1) using successive approximations given by the dynamical system
$$g_n(x)=g_{n-1}^{-1}(x)+\ep\, p(x)$$
for natural $n$ and all real $x$. The initial condition $g_0(x)\equiv x+\ep\,p(x)/2$ seems reasonable for small $\ep>0$.
However, at least for the "triangular" pdf $p=p_\triangle$ and $\ep=1/10$, this dynamics seems to be going haywire. Here are the corresponding graphs $\{\big(x,g_n(x)-g_0(x)\big)\colon|x|\le1\}$ for $n=1$ (red), $n=2$ (yellow), $n=3$ (green), $n=4$ (blue), $n=5$ (magenta):

Edit: Please feel free to replace the phrase "for each small enough $\ep>0$" in the definition of a good pdf in the first paragraph of this question by "for some real $\ep>0$". A characterization of either version of goodness, even a partial one, would be welcome.
 A: Here is the uniqueness part. Suppose you have two functions $g_1$ and $g_2$ satisfying $g(x)-g^{-1}(x)=p(x)$ with the same $p$ (red and blue lines above the diagonal; the lines below the diagonal correspond to $g_1^{-1}$ and $g_2^{-1}$).

Then for every $x$, we have $|g_1(x)-g_2(x)|=|g_1^{-1}(x)-g_2^{-1}(x)|$, so the (non-oriented!) areas between the red and the blue curves squeezed between two green lines above and below the diagonal are the same by the Cavalieri principle (vertical cross-sections have the same length). On the other hand, the areas squeezed between the orange lines are also the same by symmetry. Thus, the part $A_+$ of the area squeezed between the green and the orange lines above the diagonal at the right point on the picture, the part $A_-$ of the area squeezed between the green and the orange lines below the diagonal at the right point on the picture, and the corresponding areas $B_{\pm}$ for the left point satisfy the equation
$$
A_+-B_+=B_--A_-
$$
(just look at what is added and what is subtracted when you go from the areas between the orange lines to those between the green ones).
Hence $A_++A_-=B_++B_-$. But the right hand side tends to $0$ when the left point is moved to $-\infty$, so $A_+=A_-=0$, which, in particular, means that $g_1=g_2$ at the right point. Since the right point is arbitrary, we are done.
Existence  As Mike correctly noticed, if you have some point $Z_0=(x_0,y_0)$ on the graph of $g$, it generates a sequence of points $Z_k=(x_k,y_k)$ lying alternately on the graphs of $g$ and $g^{-1}$ and given by the recursion $x_{k+1}=x_k, y_{k+1}=y_k-p(x_k)$ when $k$ is even and $x_{k+1}=x_k-p(y_k), y_{k+1}=y_k$ when $k$ is odd.

Such a sequence can be considered for any point $Z_0$ and the property we want is that all intervals $[Z_k,Z_{k+1}]$ cross the diagonal $x=y$. We are going to show now that if $\gamma$ is any curve that goes top left to bottom right (so when you go along $\gamma$, the $x$ and the $y$ coordinates are changing in opposite directions) and $p$ is non-decreasing, then for each $n\ge 0$ the set of points such that all intervals $[Z_k,Z_{k+1}],\quad k=0,\dots,n$ intersect the diagonal is an interval on $\gamma$ and its endpoints are determined by the conditions that $Z_{n}$ or $Z_{n+1}$ hit the diagonal. The proof is just induction in $n$. Assume for simplicity that everything is smooth and the point $Z_0(t)$ moves along $\gamma$ as above.

Then, for $n=0$ notice that $\dot x_1=\dot x_0$ and $\dot y_1=\dot y_0-p'(x_0)\dot x_0$, so $\dot y_1-\dot x_1=\dot y_0-(1+p'(x_0))\dot x_0$, i.e., the coordinates $x_1,y_1$ also move in opposite directions and the "cross-diagonal component" of the speed of $Z_1$ goes in the samne direction as that of $Z_0$ and is larger in absolute value than that of $Z_0$. Starting with $Z_0$ on the diagonal and moving it up/left, we see that there will be a unique moment when $Z_1$ intersects the diagonal and the interval between these two moments is what we are looking for.
Now assume that the statement is true for some $n$. Then we apply the same argument to the points $Z_{n+1}$ and $Z_{n+2}$ and just notice that the endpoint of the addmissible interval for $Z_0$ that comes from $Z_{n+1}$ hitting the diagonal is still there, but when $Z_{n}$ hits the diagonal, the point $Z_{n+2}$ is already on the wrong side of it because each next interval is shorter than the previous one (again, due to the fact that $p$ is increasing), so the new admissible interval is determined by $Z_{n+1}$ and $Z_{n+2}$.
Finally, since the cross-diagonal speed increases all the time and the admissible length tends to $0$ as se go left and down along the diagonal, the admissible intervals shrink to one point as $n\to\infty$. Consider now the set $G$ of points $Z_0$ such that the generated sequence $Z_k$ is good in the sense that all intervals $[Z_k,Z_{k+1}],\quad k=0,1,\dots$ intersect the diagonal. Each horizontal and vertical line up to $0$ will have exactly one point in $G$. Also $G$ is closed and locally bounded. Thus, $G$ is a graph of a continuous increasing function. By symmetry, the reflection of $G$ about the diagonal is also a set of good points (but now with the first horizontal move), so we have our equation due to the same uniqueness.
That finishes what I wanted to post but the question of Mateusz if there are (reasonable) functions $p$ for which the equation has a solution for $\varepsilon p$ with "many" $\varepsilon$ still remains unanswered. His own example is degenerate (the resulting function is discontinuous), so it certainly shows that the problem is not trivial but doesn't yet settle it to my satisfaction. As to Mike's doubt about why the left and the right solutions for the triangular functions do not fit together, the story is quite simple: the values of $p$ for $x<0$ determine uniquely the red part of the graphs of $g$ and $g^{-1}$ on the picture below while the values for $x>0$ determine the blue parts and near $0$ there is a clear disagreement area.

A: The triangular function:
as fedja remarks in his comment,  when $p(x)$ is supported in $[0,1]$ and symmetric around 1/2, there is a solution from $\infty$ and a solution from $-\infty$.    The condition that they coincide is  that  $1-g^{-1}(1-x)$ also satisfies  (1) so  $g^{-1}(x) = 1-g(1-x)$ and therefore $$\epsilon p(x) = g(x) + g(1-x) -1   (2)$$.
Using this equation find  $g(1/2) = (1 + \epsilon p(1/2))/2$.  This incidentally gives us $$ g^{-1}((1 + \epsilon p(1/2))/2)$$ and using equation (1) find $g(g(1/2)$.  Iterate, finding the value of g at all the points $g(1/2), g(g(1/2),...$  and iterating backwards, which you can using $g^{-1}$ in place of g and making the appropriate changes, find g at all values $g^n(1/2)$.  Then, I claim, $g$ and $g^{-1}$ are linear in between.
$$$$  Interpolate g linearly on $(1/2, g(1/2)$.  p(x) is linear on this interval, so this wants to force $g^{-1}$ to be linear on the same interval. However, there is already a definition of $g^{-1}$, and I need to show that they are the same.  To solve the equation $g^{-1}(x) = y$ for $x \in (1/2, g(1/2)$ we need to know for what values of $y$ is $g(y) \in (1/2, g(1/2)$, and that is the interval $(g^{-1}(1/2), 1/2)$  but we have already defined $g^{-1} $ to be linear on this interval.  As the correct relations hold at the  endpoints, and  all functions are linear in between,  this definition of $g$ works. $$$$
There is a lacuna I don't know how to deal with, this all seems good provided that $g^n(1/2) \rightarrow 1$ as $n \rightarrow \infty$.  This can't always be true, but I don't know when what I have sketched above fails.
