A new path integration (PI) method is studied to improve the efficiency of computing probability density of nonlinear oscillators subject to combined harmonic and random excitation. The new PI method utilizes Fourier series to obtain a spectral presentation of the short time transition probability density (STPD) in the time domain and uses the method of linear least squares to determine the Fourier coefficients. It also utilizes a tensor product spline interpolation method to obtain an accurate representation of the STPD in the state space. The new PI method is applied to the monostable and bistable Duffing oscillators to predict the response statistics, including the time average of asymptotic mean squares and the spectral amplification factor. Specifically, the spectral amplification is used to characterize stochastic resonance of the bistable oscillator. The predictions show good agreement with Monte Carlo simulations and available data in the literature. The new PI method is also used to investigate the influence of noise intensity on stochastic P-bifurcation of the bistable oscillator. Finally, a case study shows that the new PI method reduces the computational time by 1–2 orders of magnitude in comparison with a traditional PI method.