Abstract A harmonic wavelets based technique is developed for determining the evolutionary power spectrum (EPS) matrix of the response of non-linear chain-like MDOF structural systems subject to a multi-component non-stationary stochastic excitation. Specifically, first a relationship between the evolutionary power spectrum matrices of the excitation and of the system response is derived by using a recently proposed locally stationary wavelet based representation of non-stationary stochastic processes. The relationship can be construed as a direct extension of the celebrated spectral input–output relationship of the linear stationary random vibration theory. Further, a harmonic wavelets based statistical linearization technique is proposed for the case of MDOF non-linear systems with chain-like architecture systems and hysteretic non-linearities. Numerical examples include MDOF non-linear systems comprising the versatile Bouc–Wen hysteretic model. Pertinent Monte Carlo simulations demonstrate the reliability of the technique.