The object of the study is a fast simulation method for generation and analysis of the plastic response of a randomly excited MDOF oscillatro with several potential elements with elasto-plastic constitutive behavior. The oscillator is statically determinate with linear damping. The external excitation is modeled as a Gaussian with noise process or a filtered Gaussian with noise process. Since the excitation is a radom process, the response of the oscillator is also a random process. An exact analytical solution for the distribution properties of the reponse does not exist on the basis of the governing equations of motion. One of the possible approaches to obtain statistical information about the response of the oscillator is to use direct numerical integration of the equations of motion for each simulated realization of the excitation process. However, this approach is consuming much time. An alternative faster method is based on the so-called Slepian model process concept. As the direct simulation procedure the Slepian process based algorithm is of Monte Carlo simulation type. Instead of simulating the individual response realizations the focus is solely on those parts of the response that are interesting for the plastic displacement process. These parts are those that are in the vicinity of the points of crossing of the response from the elastic domain to the plastic domain. The principle of the method is to use the known response properties at the crossing points of the stationary Gaussian white noise excited linear oscillator obtained from the elasto-plastic oscillator by totally removing the plastic domain. Thus the key to the applicability of the method is that the oscillator has a linear domain within which the response stays for a sufficiently long time to make the random response behave approximately as a stationary Gaussian process. This requires that the standard deviation of the stationary response is not too large as compared to the plastic yield limits. The Slepian model process for the behavior of the linear response is then simply the conditional mean (linear regression) of the process conditional on the position at the crossing level and a given entrance velocity to the plastic domain plus a Gaussian residual process. This description is combined with a model for the distribution properties of the entrance velocity to the plastic domain. Having this, the initial conditions are set for the movement of the elasto-plastic oscillator in the plastic domain. The simulated data are used to obtain the distribution of the plastic displacement increments as well as the distribution of the time intervals between the occurrences of these increments. From these distributions approximate distribution properties of development of the plastic displacements with respect to the time can be obtained in a reasonably simple way. The results are presented in terms of probability density and distribution functions, and the results obtained by Slepian simulations are compared with the results obtained by direct numerical integration simulation. The comparative analysis shows very good agreement. The study is exemplified for plane multi-story shear frames with rigid traverses where several of the connecting columns have finite symmetrical or non-symmetrical yield limits. In principle the method is the same for any statically determinate structural MDOF system with uniaxial elasto-plastic constitutive behavior of the elements of the structure. Subject to an empirically based assumption that there is a localization tendency of close to simultaneously yielding elements, the Slepian method is successfully extended in simplified form to a manageable approximate method applicable to cases with a very large number of degrees of freedom.