Final update on 11/29/2019: I have worked on this a bit more, and wrote an article summarizing all the main findings. You can read it here.
I am looking for a solution where $f(z)\geq 0$ and $\int_{-\infty}^{\infty} f(z) dz = 1$. I believe there is only one solution. That solution also satisfies $f(z) = f(-z)$ and $f(0) = f(-1) = f(1)$.
In fact $f$ is the density of $Z$, with
$$ Z= X_1 + X_1 X_2 + X_1 X_2 X_3 +\cdots $$ with the $X_i$'s being i.i.d with the following distribution:
$$P(X_i= -1) = P(X_i =-0.5) = P(X_i = 0.5) = P(X_i = 1) = 0.25.$$
We also have
$$E(Z^k) = \int_{-\infty}^{\infty} z^k f(z) dz = \frac{1}{2} \Big(1+ \frac{1}{2^k}\Big) \int_{-\infty}^{\infty} (1+z)^k f(z)dz$$ if $k$ is even, or zero otherwise. Note that the characteristic function of $Z$ satisfies $$\frac{d^k \psi_Z}{dt^k}(0) = \frac{d^k E(\exp it Z)}{dt^k}(0)= i^k E(Z^k).$$ The reason why I am interested in this problem can be found here.
The empirical percentile distribution of $Z$ is pictured below:
Update
This is a particular case of a more general problem: solving $F_Z = F_{g(X,Z)}$ where $g$ is any function, $F$ represents the distribution associated with the density $f$, and $X, Z$ are independent random variables. In this case, $g(X,Z) = X(1+Z)$. Another case involving $g(X, Z) = \sqrt{X+Z}$ was discussed here. Despite the smooth appearance of $F_Z$, it is possible that $F_Z$ is nowhere differentiable. Similar distributions have been analyzed by David Bailey in his book Experimental Mathematics in Action, published in 2007. In particular, sections 5.2 and 5.3 (pages 114-137) are very relevant to this context. One of the densities he has studied, namely $2qf(x) = f(\frac{x-1}{q}) + f(\frac{x+1}{q})$, is very similar to the one I posted on Math.Stackexchange, here, leading to the same discussion as to when it is smooth or not, an unsolved problem in many (but not all) cases. All these problems end up in some functional equations like the one discussed in this question.
A possible way to find a numerical solution is as follows. Build a sequence of densities $f_n$ that are piecewise uniform on the support domain, starting in this very particular case with $f_1(x) = \frac{1}{2}$ if $x\in [-1, 1]$ and $f_1(x) = 0$ otherwise. In fact, it's quite possible that $f$ is constant on $[-1, 1]$ based on empirical evidence. Each $f_n$ must satisfy $f_n(x) \geq 0, f_n(x) = f_n(-x)$, and $\int_{-\infty}^{\infty} f(x) dx = 1$. Of course this assumes that the density exists, and that the solution satisfying these constraints is unique. It also assumes that the algorithm in question converges to the solution.
For instance, $f_n$ might be defined using $n$ intervals $I_{n1}, \cdots, I_{nn}$, with $f_n(x) = c_{nk}$ constant $(k=1,\cdots,n)$ if $x\in I_{nk}$. The intervals $I_{nk}$'s and the constants $c_{nk}$'s are chosen so as to minimize some error criterion $E_n$, for instance
$$E_n = \sup |F_{Z}^{(n)}(x) - F_{g(X,Z)}^{(n)}(x)| \mbox{ or } E_n = \int_{-\infty}^{\infty} |F_{Z}^{(n)}(x) - F_{g(X,Z)}^{(n)}(x)| dx.$$
Here $F^{(n)}$ is the distribution attached to $f_n$. Techniques about how to solve this problem are described in my article New Perspectives on Statistical Distributions and Deep Learning.
I will add an illustration in the next few days, if I find the time.
