There are many *local* solutions of this equation.  For example, suppose that one starts with a $C^2$ function $f$ on an interval $I\subset\mathbb{R}$ such that $f'$ is positive on $I$ and $f(I)$ is disjoint from $I$.  Then an inverse $g:f(I)\to I$ of $f:I\to f(I)$ exists and is $C^2$.  Now define $f$ on the interval $f(I)$ so that $f(y) = f''(g(y))$ for $y\in f(I)$.  Then for $x\in I$, we will have $x = g(y)$ for some $y\in f(I)$ and, of course, $y = f(x)$.  Then $f''(x) = f''(g(y)) = f(y) = f(f(x))$ for all $x\in I$.

These sorts of 'rough' solutions are constructed without any fixed points.  Solutions with fixed points are much more rigid.  A $C^2$ solution on an open domain $D$ such that $f(D)\subset D$ must be smooth on $D$, since $f''=f{\circ}f$, implying that if $f$ is $C^k$, then $f$ must be $C^{k+2}$.  In fact, with a little effort, one can show that a $C^2$ solution with an contracting fixed point must be real-analytic in a neighborhood of the fixed point, since the equation $f''=f{\circ}f$ allows one to prove an estimate of the form $|f^{(k)}|\le C^k\,k!$ for some constant $C$ on a neighborhood of the fixed point.

**Note 1:**  For every constant $b\in\mathbb{C}$, there is a unique formal power series with lowest order term $bz$ that satisfies $f''(z) = f(f(z))$.  The first few terms are 
$$
f(z) = bz+\frac{{b}^{2}}{3!}\,{z}^{3}
+{\frac {{b}^{3} \left( {b}^{2}{+}1 \right)}{5!}}\,{z}^{5}
+{\frac {{b}^{4} \left( {b}^{6}{+}{b}^{4}{+}11\,{b}^{2}{+}1\right)}{7!}}\,{z}^{7}+\cdots.\tag1
$$
When $|b|<1$, this series converges absolutely and uniformly on the disk $|z|^2\le 6\bigl(1{-}|b|\bigr)$, and satisfies $|f(z)|\le |z|$ there.  See the Addenedum below for a sharper (but still not sharp) estimate of the radius of convergence.

Though I don't (yet) have a proof, numerical calculations indicate that, when $b$ is a small negative real number, the above function $f$ extends real-analytically and *periodically* to the entire real line and gives a solution $f:\mathbb{R}\to\mathbb{R}$.

**Addendum to Note 1:**  One can show that, when $0<b<1$, the real-analytic odd function $f$ that equals the power series (1) on its interval of convergence extends real-analytically to a bounded interval $\bigl(-r(b),r(b)\bigr)\subset\mathbb{R}$ on which $|f(x)|<|x|$ and that $\lim_{x\to r(b)^-}f(x)=r(b)$.  In particular, $f:\bigl(-r(b),r(b)\bigr)\to\bigl(-r(b),r(b)\bigr)$ is a real-analytic diffeomorphism with a single contracting fixed point at $x=0$.  Moreover, the formal series converges to $f$ uniformly on compact subsets of $\bigl(-r(b),r(b)\bigr)$, and $f$ cannot be extended real-analytically to any larger interval.  (There are some indications that $f$ may extend *smoothly* beyond $x = r(b)$, in which case, $x=r(b)$ would become an expanding fixed point of $f$.)  Also, $r:(0,1)\to(0,\infty)$ is a continuous, decreasing bijection, and
$$
\frac{\sqrt{6\bigl(1{-}b\bigr)}}{b}
> r(b)> \begin{cases}
\sqrt{\displaystyle\frac3{2b}} & \text{for}\ 0<b\le\tfrac12,\\
\\
\sqrt{6(1{-}b)} & \text{for}\ \tfrac12\le b<1,
\end{cases}
$$
from which it follows that, for $b=1$, the radius of convergence of the series is $0$.

**Note 2:**  More generally, for any two constants $a,b\in\mathbb{C}$, there is a formal power series
$$
f(z) = a+b\,(z{-}a) 
+\frac{a}{2!}\,(z{-}a)^2
+\frac{b^2}{3!}\,(z{-}a)^3
+\frac{ab(b{+}1)}{4!}\,(z{-}a)^4
+\cdots\tag2
$$
that has $a$ as a formal fixed point, i.e., $f(a) = a$, so that the composition $f(f(z))$ makes sense as a power series centered at $z = a$
and, formally, $f'(a) = b$, that satisfies $f''(z) = f(f(z))$ as formal power series centered at $z = a$.  Moreover, this is the unique power series centered at $z=a$ that has $f(a) = a$ and $f'(a) = b$ and satisfies $f''(z) = f(f(z))$ as formal power series.  

As in the case $a=0$, when $|b|<1$, so that $f$ is a 'formal contraction' on a neighborhood of $a$, it turns out that the series converges absolutely and uniformly on a disc of the form $|z-a| \le r(a,b)$ for some $r(a,b)>0$, so this gives a two-parameter family of local solutions with a contracting fixed point.  It remains to be seen whether there are values of $(a,b)$ (other than $(0,0)$) for which the corresponding $f$ extends to an entire holomorphic function on $\mathbb{C}$, or whether there exist nontrivial values of $(a,b)\in\mathbb{R}^2$ for which $f$ extends analytically to a neighborhood of $\mathbb{R}\subset\mathbb{C}$.

**Note 3:** The (two) multivalued solutions described by Michael Engelhardt have fixed points and hence are (analytic continuations of) solutions of the type (2).  One can see this as follows:  These (multivalued) solutions can be written in the form 
$$
f(x) = i\sqrt{2}\,\left(\frac{x}{i\sqrt{2}}\right)^b,\qquad
\text{where}\  b = \tfrac12(1\pm i\sqrt{7}).
$$
Clearly, $a\in\mathbb{C}$ will be a fixed point, i.e., $f(a) = a$ if and only if 
$$
1 = \left(\frac{a}{i\sqrt{2}}\right)^{b-1},
$$ 
and this happens (for $b = \tfrac12(1+i\sqrt7)$) when, for some integer $k$, 
$$
a = a_k = i\sqrt{2}\, e^{i\pi k(1+i\sqrt7)/2} = i^{k+1}\sqrt{2}\,\left(e^{-\pi\sqrt7}\right)^{k/2}.
$$
Moreover, we have
$$
f'(a_k) = b\left(\frac{a_k}{i\sqrt2}\right)^{b-1} = b,
$$
so $|f'(a_k)| = |b| = \sqrt 2>1$, which implies that the fixed point is a repelling fixed point.  

This is interesting because it implies that the formal power series (2) for $(a_k,b)$ must have a positive radius of convergence, even though $|b|>1$.  This led me to speculate that maybe the formal power series (2) might have a positive radius of convergence for any $(a,b)\in\mathbb{C}$, but Will Sawin (in a comment below) pointed out that this cannot be true.