Time-dependent autoregressive moving-average (TARMA) models, with parameters belonging to a subspace spanned by pre-defined time-domain functions, offer appropriate representation for a fairly wide class of non-stationary signals. The methods available for their estimation have, however, been limited to the subclass of pure autoregressive (TAR) models. In this paper, a polynomial-algebraic (P-A) method for the estimation of mixed TARMA models is introduced. The method is based upon a skew polynomial operator algebra, through which the AR and MA polynomial operators are related to the model’s inverse function, and through which appropriate filtering operations, essential for the construction of a restricted quadratic approximation of the prediction error (PE) criterion, are introduced. The P-A method is characterized by low computational complexity, no need for initial guess parameter values, and avoidance of local extrema problems associated with PE optimization. The performance characteristics of the P-A method are assessed via numerically synthesized non-stationary signals, while PE-based model refinement is also examined. A first application of the non-stationary TARMA approach and the P-A method to the problem of modeling and prediction of an actual non-stationary engineering signal is also presented, and critical comparisons with non-stationary ARIMA (integrated ARMA) and conventional ARMA methods made. The paper is divided into two parts. The polynomial-algebraic method is introduced in the first, whereas its application to the non-stationary modeling and prediction of an automotive active suspension power consumption signal paradigm is presented in the second.