Find the Taylor series of $f^{\circ n}=\underset{n\text{ times}}{\underbrace{f\circ f \circ f\ldots \circ f}}$ Define $f(x)=ln(1+x)$. Then $f^{\circ 2}(x)=ln(1+ln(1+x))$, and $f^{\circ 3}(x)=ln(1+ln(1+ln(1+x)))$, etc.
Find the Taylor series of $f^{\circ n}=\underset{n\text{ times}}{\underbrace{f\circ f \circ f\ldots \circ f}}$ about $x=0$ up to order $x^2$.

How am I supposed to take a derivative of this function? If I use the chain rule the derivatives gets unwieldy very fast, unless I'm mistaken.
Furthermore, how do I go about finding the Taylor series expansion for this composition? 
 A: After computing the first few,$^{(\dagger)}$ we "guess" the expansion to order $x^2$ (at $0$) will be of the form:
$$
f^{\circ n}(x) = x - \frac{n}{2} x^2 + o(x^2)\tag{1}
$$
so let us prove it by induction on $n\geq 1$.


*

*For $n=1$, it is clear from the known expansion of $\ln(1+x)$.

*Assume it holds for some $n\geq 1$. then
$$\begin{align}
f^{\circ (n+1)}(x) &= \ln(1+ f^{\circ n}(x)) \stackrel{(1)}{=} \ln( 1+x - \frac{n}{2} x^2 + o(x^2)) \\&= 
x - \frac{n}{2} x^2 + o(x^2) - \frac{1}{2}\left(x - \frac{n}{2} x^2 + o(x^2)\right)
= 
x - \frac{n}{2} x^2 - \frac{x^2}{2} + o(x^2)
\\&= 
x - \frac{n+1}{2} x^2 + o(x^2)
\end{align}$$
and by the induction principle we can conclude.

$(\dagger)$ I personally computed the first $6$ with Mathematica, and from there the pattern seems to be 
$$
f^{\circ n}(x) = x-\frac{n}{2}x^2+\frac{n(3n+1)}{12}x^3 + o(x^3)
$$
A: According to Faà di Bruno's formula and given
$f(x) = \ln(1 + x) = \sum_{j=1}^{\infty} (-1)^{j-1} \frac{x^j}{j!}$,
$$(f \circ f)^{(n)} = \sum_{k=1}^{n} f^{(k)}(f) B_{n,k}(f',f'',...,f^{n-k+1})$$
where
$$B_{n,k}(y_1,...,y_{n-k+1}) = \sum_{\begin{matrix}j_1,...,j_{n-k+1} \geq 0\\j_1 + \cdots + j_{n-k+1} = k\\ j_1 + 2 j_2 + \cdots + (n-k+1) j_{n-k+1} = n\end{matrix}} \frac{n!}{j_1! \cdots j_{n-k+1}!} \prod_{i = 1}^{n-k+1} \left(\frac{y_i}{i!}\right)^{j_i}$$
is the Bell polynomial. Hence,
\begin{align*}
(f \circ f)^{(n)}(0) &= \sum_{k=1}^{n} f^{(k)}(f(0)) B_{n,k}(f'(0),f''(0),...,f^{n-k+1}(0))\\
&= \sum_{k=1}^{n} f^{(k)}(0) B_{n,k}(1,-1,...,(-1)^{i-1}(i-1)!...,(-1)^{n-k} (n-k)!)\\
&= \sum_{k=1}^{n} (-1)^{k-1}(k-1)! \left(\sum_{\begin{matrix}j_1,...,j_{n-k+1} \geq 0\\j_1 + \cdots + j_{n-k+1} = k\\ j_1 + 2 j_2 + \cdots + (n-k+1) j_{n-k+1} = n\end{matrix}} \frac{n!}{j_1! \cdots j_{n-k+1}!} \prod_{i = 1}^{n-k+1} \left(\frac{(-1)^{i-1}(i-1)!}{i!}\right)^{j_i}\right)\\
&= \sum_{k=1}^{n} (-1)^{k-1}(k-1)! \left(\sum_{\begin{matrix}j_1,...,j_{n-k+1} \geq 0\\j_1 + \cdots + j_{n-k+1} = k\\ j_1 + 2 j_2 + \cdots + (n-k+1) j_{n-k+1} = n\end{matrix}} \frac{n!}{j_1! \cdots j_{n-k+1}!} (-1)^{\overbrace{\sum (i-1) j_i}^{n-k}} \prod_{i = 1}^{n-k+1} \left(\frac{1}{i}\right)^{j_i}\right)\\
&= (-1)^{n-1} n! \sum_{k=1}^{n}  (k-1)! \left(\sum_{\begin{matrix}j_1,...,j_{n-k+1} \geq 0\\j_1 + \cdots + j_{n-k+1} = k\\ j_1 + 2 j_2 + \cdots + (n-k+1) j_{n-k+1} = n\end{matrix}} \frac{1}{j_1! \cdots j_{n-k+1}! \prod_{i = 1}^{n-k+1} i^{j_i}} \right)
\end{align*}
so the Taylor series for $f \circ f$ is already complicated. I do not think there is closed form (i.e. in term of $n$, without summation) for the above formula.
A: We can write this recursively:
$$(f^{\circ (n-1)} \circ f) (x) = \sum_{k=0}^{2} c_k{f^{\circ (n-1)}(x)}^k$$
where $c_k$ is the coefficients for taylor expansion of $\log(1+x)$ : $c = [0,1,-1/2,1/3]$, (index starting at 0).
Now each $f^{\circ (n-1)}(x)^k$ is an iterated convolution $k$ times of the vector of coefficients. This may sound nasty, but keep in mind that we are free to at each convolution truncate at our $x^3$ terms. This reduces complexity a lot. We can also do this in the Fourier domain by multiplication instead where it explicitly will become a sum of a sum of monomials of the FFT of the coefficients for $\log(x+1)$. And since an FFT is just a linear switch of basis and sum of monomials is a polynomial this makes the other answers' findings theoretically reasonable that we should be getting polynomials in $n$.
First few acquired in Matlab / Octave:
$$\left[\begin{array}{rrrr}
0&1&-0.5&0.\overline {3}\cdots\\
0&1&0.5&-1.1\overline{6}\cdots\\
0&1&1.5&-0.\overline{6}\cdots\\
0&1&2.5&1.8\overline{3}\cdots\\
0&1&3.5&6.\overline{3}\cdots
\end{array}\right]$$
with $\overline a$ meaning digit $a$ is repeated infinitely. Now, by looping the line 
    f = f + c(2)*conv(f,f)(1:4) + c(3)*conv(conv(f,f),f)(1:4) + c(4)*conv(conv(conv(f,f),f),f)(1:4);

We see the $x$ constant, $x^2$ linear and if we do polynomial regression on $x^3$ coefficient we will probably find it being second degree polynomial.
