We propose a global smoothing method based on polynomial splines for the estimation of functional coefficient regression models for non-linear time series. Consistency and rate of convergence results are given to support the proposed estimation method. Methods for automatic selection of the threshold variable and significant variables (or lags) are discussed. The estimated model is used to produce multi-step-ahead forecasts, including interval forecasts and density forecasts. The methodology is illustrated by simulations and two real data examples.