A bistable Duffing oscillator subjected to additive and multiplicative Ornstein–Uhlenbeck (OU) colored excitations is examined. It is modeled through a set of four first-order stochastic differential equations by representing the OU excitations as filtered Gaussian white noise excitations. Enlargement in the state-space vector leads to four-dimensional (4D) Fokker–Planck–Kolmogorov (FPK) equation. The exponential-polynomial closure (EPC) method, proposed previously for the case of white noise excitations, is further improved and developed to solve colored noise case, resulting in much more polynomial terms included in the approximate solution. Numerical results show that approximate solutions from the EPC method compare well with the predictions obtained via Monte Carlo simulation (MCS) method. Investigation is also carried out to examine the influence of intensity level on the probability distribution solutions of system responses.