Efficient iteration of arbitrary non-linear dynamical system (deterministic or stochastic)

For a while I have been wondering what the most efficient way to simulate an arbitrary non-linear (deterministic of stochastic) dynamical system in Python. I end up doing this a lot either teaching or research. I am convinced that there must be a simple and efficient way of doing this.

At the pub this evening I came up with the following...

def iterate(F, X, T, **params):
    """Iterate a non-linear map F starting from some initial condition X for T periods."""
    t = 0
    while t < T:
        yield X
        X = F(X, **params)
        t += 1

...a test case using the Tinkerbell Map...

def tinker_bell_map(X, a, b, c, d):
    return [X[0]**2 - X[1]**2 + a * X[0] + b * X[1], 2 * X[0] * X[1] + c * X[0] + d * X[1]]

...yields...

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 10, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 26 µs per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 100, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 254 µs per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 1000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 2.36 ms per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 10000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 19.6 ms per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 100000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 192 ms per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 1000000, a=0.9, b=-0.6013, c=2.0, d=0.5)] 
1 loops, best of 3: 2.02 s per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 10000000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 20.5 s per loop

...I have tried several other test cases for deterministic and stochastic systems and the above works like a charm. While I think the above is pretty good, I wonder if it could be made even faster using Numba? I tried doing

from numba import jit

@jit
def iterate(F, X, T, **params):
    """Iterate a non-linear map F starting from some initial condition X for T periods."""
    t = 0
    while t < T:
        yield X
        X = F(X, **params)
        t += 1

But this didn't lead to any real speedups. Thoughts on how to speed this up in Numba? Thoughts on a better way to accomplish the above? Thoughts on whether something like this would be useful for Quantecon?

The issue names no repository files, tests, or entry points; begin by reviewing the proposed iterate function and its Numba experiment. Done would require an agreed implementation and measurable speedup, or a documented decision about whether this belongs in QuantEcon.