In this thesis a novel approximate/analytical approach based on the concepts of stochastic averaging and of statistical linearization is developed for the response determination of nonlinear/hysteretic multi-degree-of-freedom (MDOF) systems subject to evolutionary stochastic excitation. The significant advantage of the approach relates to the fact that it is readily applicable for excitations possessing even non-separable evolutionary power spectra (EPS) circumventing ad hoc pre-filtering and pre-processing excitation treatments associated with existing alternative schemes of linearization. Further, the approach can be used, in a rather straightforward manner, in conjunction with recently developed design spectrum based analyses for obtaining peak response estimates without resorting to numerical integration of the nonlinear equations of motion. Furthermore, a novel approximate/analytical Wiener path integral based solution (PIS) is developed and a numerical PIS approach is extended to determine the response and first-passage probability density functions (PDFs) of nonlinear/hysteretic systems subject to evolutionary stochastic excitation. Applications include the versatile Preisach hysteretic model, recently applied in modeling systems equipped with smart material (shape memory alloys) devices used for seismic hazard risk mitigation. The approach is also applied to determine the capsizing probability of a ship, whose rolling dynamics is captured by a softening Duffing oscillator. Finally, novel harmonic wavelets based joint time-frequency response analysis and identification approaches are developed capable of determining the time-varying frequency content of non-stationary complex stochastic phenomena encountered in engineering applications. Specifically, a harmonic wavelets based statistical linearization approach is developed to determine the EPS of the response of nonlinear/hysteretic systems subject to stochastic excitation. In a similar context, an identification approach for nonlinear time-variant systems based on the localization properties of the harmonic wavelet transform is also developed. It can be construed as a generalization of the well established reverse multiple-input/single-output (MISO) spectral identification approach to account for non-stationary inputs and time-varying system parameters. Several linear and nonlinear time-variant systems are used to demonstrate the reliability of the approach.