Abstract: Determination of stochastic response of a nonlinear dynamic systme endowed with fractional-order derivative element has always been the difficulty of the random vibration research.A multi-harmonic balance method is presented for determining the stochastic response of a nonlinear Duffing oscillator endowed with fractional derivative elements and subject to stochastic excitation.Specifically, the fractional-order nonlinear differential equation is transferred into a set of nonlinear algebra equations in terms of unknown Fourier coefficients of response by employing the multi-harmonic balance method.Next, Newton's iterative approach is adopted for determining the unknown response Fourier coefficient.Thus, the response time history can be obtained by the inverse Fourier transform.Further, the Fourier coefficients of sample responses or the response power spectral density can be calculated by repeated use of the proposed method from the Fourier coefficients of excitations sampled from its power spectral density.A closed form solution of the Fourier coefficients of the cubic term is derived to avoid the so-called aliasing effect due to system non-linearity and enhance compuational accuracy and efficiency.Finally, a pertinent numerical example demonstrates the accuracy and reliability of the proposed method for determing response power spectral density of a Duffing oscillator with different levels of nonlinearity.