A new computationally simple and precise model approximation method is described for large-scale linear discrete-time systems. By least squares matching of a suitable number of time moment proportionals and Markov parameters about (z–1) of the original higher order system within the approximate model, stable denominator polynomial coefficients of the approximate model are determined. To improvise the accuracy of the approximate model, numerator polynomial coefficients are determined by minimizing the integral squared error (ISE) between the unit impulse responses of the original system and its approximate model. A matrix formula is formulated for evaluating numerator coefficients of the approximate model that leads to minimum ISE, and also for evaluating ISE. The efficacy of the proposed method is shown by illustrating three typical numerical examples employed from the literature, and the results are compared with many familiar reduction methods in terms of the ISE and relative ISE values pertaining to impulse input. Furthermore, time and frequency responses of the original system and the respective approximate model are plotted.