Functional data analysis tools, such as function-on-function regression models, have received considerable attention in various scientific fields because of their observed high-dimensional and complex data structures. Several statistical procedures, including least squares, maximum likelihood, and maximum penalized likelihood, have been proposed to estimate such function-onfunction regression models. However, these estimation techniques produce unstable estimates in the case of degenerate functional data or are computationally intensive. To overcome these issues, we proposed a partial least squares approach to estimate the model parameters in the functionon-function regression model. In the proposed method, the B-spline basis functions are utilized to convert discretely observed data into their functional forms. Generalized cross-validation is used to control the degrees of roughness. The finite-sample performance of the proposed method was evaluated using several Monte-Carlo simulations and an empirical data analysis. The results reveal that the proposed method competes favorably with existing estimation techniques and some other available function-on-function regression models, with significantly shorter computational time.Recent advances in computer storage and data collection have enabled researchers in diverse branches of science such as, for instance, chemometrics, meteorology, medicine, and finance, recording data of characteristics varying over a continuum (time, space, depth, wavelength, etc.). Given the complex nature of such data collection tools, the availability of functional data, in which observations are sampled over a fine grid, has progressively increased. Consequently, the interest in functional data analysis (FDA) tools is significantly increasing over the years. Silverman (2002, 2006), Ferraty and Vieu (2006), Horvath and Kokoszka (2012) and Cuevas (2014) provide excellent overviews of the research on theoretical developments and case studies of FDA tools.Functional regression models in which both the response and predictors consist of curves known as, function-on-function regression, have received considerable attention in the literature. The main goal of these regression models is to explore the associations between the functional response and the functional predictors observed on the same or potentially different domains as the response function. In this context, two key models have been considered: the varying-coefficient model and the function-on-function regression model (FFRM). The varying-coefficient model assumes that the functional response Y (t) and functional predictors X (t) are observed in the same domain. Its estimation and test procedures have been studied by numerous authors, including Fan and Zhang among many others. In contrast, the FFRM considers cases in which the functional response Y (t) for a given continuum t depends on the full trajectory of the predictors X (s). Compared with the varying-coefficient model, the FFRM is more natural; therefore, we restrict our attent...