The paper presents a general algorithm for parameter identification in the transfer function matrix of a multi-input multi-output system. The method, whose basis is in the theory of distributions, employs the concept of generalized functions in the manner originally conceived and established by Dirac and Schwartz respectively, while the basis of most of the methods of system identification is to treat the process signals as functions in the ordinary sense with its differentiability limitations. The present method employs exponentially weighted series of the generalized time derivatives of the impulse distribution, known as the Poisson moment functional (PMF) expansion with certain advantages such as infinite differentiability and noise immunity. Some important aspects of identifiability and suitability of input signals in the context of the PMF approach are discussed. The algorithm is capable of treating the practically important case of data on an arbitrary but active period of time implying unknown initial conditions to be accounted for.