Recent technological developments have enabled us to collect complex and high-dimensional data in many scientific fields, such as population health, meteorology, econometrics, geology, and psychology. It is common to encounter such datasets collected repeatedly over a continuum. Functional data, whose sample elements are functions in the graphical forms of curves, images, and shapes, characterize these data types. Functional data analysis techniques reduce the complex structure of these data and focus on the dependences within and (possibly) between the curves. A common research question is to investigate the relationships in regression models that involve at least one functional variable. However, the performance of functional regression models depends on several factors, such as the smoothing technique, the number of basis functions, and the estimation method. This paper provides a selective comparison for function-on-function regression models where both the response and predictor(s) are functions, to determine the optimal choice of basis function from a set of model evaluation criteria. We also propose a bootstrap method to construct a confidence interval for the response function. The numerical comparisons are implemented through Monte Carlo simulations and two real data examples.where Φ(t) = {φ 1 (t), · · · , φ K Y (t)} and Ψ(s) = {ψ m1 (s), · · · , ψ mK m,X (s)} denote the vectors of the basis functions, c i = {c i1 , · · · , c iK Y } and d im = {d im1 , · · · , d imK m,X } are the coefficient vectors, and K Y and K m,X are the number of basis functions used for approximating the functional response and functional predictors, respectively. Similarly, the bivariate coefficient function is defined as: β m (s, t) = ∑ j,k ψ mj (s)b mjk φ k (t) = Ψ m (s)B m Φ(t),where B m = (b mjk ) j,k denote a K m,X × K Y dimensional coefficient matrix. Accordingly, the regression