An algorithm for the computation of the stochastic non-stationary non-linear response of large FE-models is presented. In stochastic analysis, the response is described by the mean and covariance function and possibly by the probability distribution of the non-linear response. To estimate the covariance matrix, the method of equivalent statistical linearization is applied for linearizing all non-linear elements. The large fully populated and symmetric covariance matrix of dimension P 2n is described by the so called Karhunen-Lo e eve expansion representation, which allows one to employ a feasible description. In this context, an efficient procedure is suggested to determine the minimal number of Karhunen-Lo e eve vectors necessary to assure a sufficiently accurate representation. This method in fact allows one to employ any available deterministic integration scheme to compute the Karhunen-Lo e eve vectors. To increase the computational efficiency, modal coordinates are used to represent the linear sub-system. Special attention is given to the effects of mode truncation which are generally not negligible for large FE-systems. Moreover, the suggested approach has no limitation w.r.t. the size of the model in terms of degrees of freedom (DOFs) or the number of non-linear elements. The feasibility of the proposed procedure is demonstrated in the numerical examples where the methodology is applied to an office building with approximately 25, 50 and 100 thousand DOFs containing hysteretic damping elements.