  We propose a nonlinear function-on-function regression model where both the
covariate and the response are random functions. The nonlinear regression is
carried out in two steps: we first construct Hilbert spaces to accommodate the
functional covariate and the functional response, and then build a second-layer
Hilbert space for the covariate to capture nonlinearity. The second-layer space
is assumed to be a reproducing kernel Hilbert space, which is generated by a
positive definite kernel determined by the inner product of the first-layer
Hilbert space for $X$--this structure is known as the nested Hilbert spaces. We
develop estimation procedures to implement the proposed method, which allows
the functional data to be observed at different time points for different
subjects. Furthermore, we establish the convergence rate of our estimator as
well as the weak convergence of the predicted response in the Hilbert space.
Numerical studies including both simulations and a data application are
conducted to investigate the performance of our estimator in finite sample.
