A statistical linearization method is proposed for determining the response of a single-degree-of-freedom Bouc-Wen system subjected to combined colored noise and harmonic loads. The proposed method is based on the assumption that the system response can be decomposed into the sum of deterministic harmonic and zero-mean random components. Specifically，the equation of motion is decomposed into two sets of nonlinear differential equations governing deterministic response and stochastic response，respectively. The harmonic balance method is used to solve the equation of motion with deterministic excitation，whereas the statistical linearization method is utilized to obtain the variance of the stochastic response. These treatments lead to a set of coupled algebraic equations in terms of the Fourier coefficients of the deterministic response and the stochastic response variance. Standard numerical schemes such as Newton’s iteration method are adopted to solve the preceding non-linear algebraic equations. Pertinent numerical examples demonstrate the applicability and accuracy of the proposed method.