In this paper, the problem of approximating a linear time-invariant stable system by a finite weighted sum of given exponentials is considered. System identification schemes using Laguerre models are extended and generalized to Kautz models, which correspond to representations using several different possible complex exponentials. In particular, linear regression methods to estimate this sort of model from measured data are analyzed. The advantages of the proposed approach are the simplicity of the resulting identification scheme and the capability of modeling resonant systems using few parameters. The subsequent analysis is based on the result that the corresponding linear regression normal equations have a block Toeplitz structure. Several results on transfer function estimation are extended to discrete Kautz models, for example, asymptotic frequency domain variance expressions.