One way of dealing with colored noise is to represent it by another SDE. For example, <a href="http://www.itpa.lt/~ruseckas/files/presentations/Turkey-2012.pdf">here is a presentation which presents some models</a> like:

$$dx = (1+\frac{1}{2x} - \frac{x}{2\times 10^3})x^{4} dt + x^{\frac{5}{2}} dW_t $$

has a solution with a $1/f$ spectrum. So now you just have to solve a system of nonlinear SDEs.

I wouldn't use the Strong Order 2 Method in practice because it takes so many more computations I don't think it actually gives a speedup. Rossler has some Runge-Kutta type methods which are easier to implement. But honestly try Euler-Maruyama / Milstein first and see if that's good enough for your needs before you do something fancy.

----------

Update: I came back to this as my own notes and noticed that the link was broken. Here's [another paper with an SDE model which gives $1/f^\beta$ spectra](http://ac.els-cdn.com/S0378437106000574/1-s2.0-S0378437106000574-main.pdf?_tid=104ca482-93c9-11e6-9927-00000aab0f26&acdnat=1476640425_a2513e09684a220f1dfd041e78d2c0cb).