Thanks Chris. That sounds like an elegant way to do it.

I have been doing some research into the paper you posted as well as other papers by the same authors. In particular I implemented the diff eq. given in Eq. (14) of this follow up: [1111.2995] Tsallis distributions and 1/f noise from nonlinear stochastic differential equations, i.e.,

dx = \sigma^2 ( \eta - \lambda/2) (x^2 + x_0^2)^{\eta - 1} x dt + \sigma (x^2 + x_0^2)^{\eta/2} dW

Unlike the other equations in this paper and the one you linked to, this is the only equation that allows negative x (according to the authors). I would like to have a noise process that describes symmetric fluctuations around a mean value, so this equation seems to be the best choice, for that reason.

However, I have some challenges in the numerical implementation. Occasionally, this equation will produce solutions with sudden extremely large (positive or negative) values, which causes very slow integration. My current implementation is so slow that it is not practical; I can’t produce a decent number of trajectories in a reasonable time. I’m wondering if the situation can be improved by choosing better solver options, or perhaps modifying the equation to prevent very large (absolute) values, without modifying the noise spectrum. Here is a minimal code example, with the solver options I am using:

```
λ = 3 # 1/f case
σ = 1.0
η = 5/2 # somewhat arbitrary choice, same as in the article
x0 = 1.0
f(x, p, t) = σ^2 * (η - λ/2) * (x0^2 + x^2)^(η - 1) * x
g(x, p, t) = σ * (x0^2 + x^2)^(η / 2)
prob = SDEProblem(f, g, 0.0, (0.0, 1.0))
sol = solve(prob, SOSRI(), maxiters=1e9, dtmin=1e-16);
```

In the article I linked to above, they don’t give a lot of details about the numerical solver, but in [1003.1155] Modeling scaled processes and 1/f^b noise by the nonlinear stochastic differential equations [Eq. (36)], the use an x-dependent time-step \Delta t = \kappa^2 / (x_0 + x)^l, with \kappa \ll 1 and l parameters. I also implemented such a function using a `StepsizeLimiter`

(let me know if this looks right):

```
κ = 0.1
l = 2η - 2
dtFE(x, p, t) = abs(κ^2 / (x0 + x)^l) # function to set step size with StepsizeLimiter
step_size_limiter = StepsizeLimiter(dtFE; safety_factor=1, max_step=true)
sol = solve(prob, EM(), maxiters=1e9, callback=step_size_limiter, adative=false, dt=1e-5);
```

seems to work, but I regularly get a warning `Warning: Instability detected. Aborting`

. Not sure if I could simply discard the aborted trajectories, as one possible work around.