Based on Volterra series theory, a Volterra-pseudo excitation method (Volterra-PEM) is developed to compute the response power spectral density (PSD) for nonlinear systems under random excitation. Firstly, within the framework of Volterra series theory, an improved pseudo excitation method (PEM) is established for multi-dimensional power spectral density (MPSD) analysis of nonlinear systems. As a generalized PEM for linear random vibration analysis, the Volterra-PEM is used to analyse the response MPSD, which also has a very concise expression. Secondly, in the case of computation difficulties when multi-dimensional integrating from MPSD to PSD, the computational efficiency is improved by converting the multi-dimensional integral into a matrix operation. Finally, as numerical examples, the Volterra-PEM is used to estimate the response PSD for stationary random vibration of a nonlinear spring-damped oscillator and a non-ideal boundary beam with geometrical nonlinearity. Compared with Monte Carlo simulation (MCS), the results show that the proposed method can predict the response PSD for nonlinear systems quickly and accurately under appropriate excitation intensities.