Stochastic differential equations: Python+Numpy vs. Cython. FIGHT!!

Cython

# Same with Cython:
import numpy as np
cimport cython
cimport numpy as np

from libc.stdint cimport uint32_t, int32_t
from libc.math cimport sqrt
from libc.math cimport fabs
from libc.math cimport pow


ctypedef float     real_t
ctypedef uint32_t  uint_t
ctypedef int32_t   int_t

cdef:
    real_t f1 = .1
    real_t g1 = .01
    real_t g2 = .1

cdef real_t f(real_t X) nogil:
    return -f1*X
cdef real_t g(real_t X) nogil:
    return sqrt((g1 + g2*X*X))
cdef real_t scale(real_t X) nogil:
    return 1/fabs(g1 + g2*X*X)

@cython.boundscheck(False)
@cython.wraparound(False)
def simulation_c(real_t L = 0, uint_t N = 100000, real_t dt = .001, real_t init = .1):
    
    cdef:
        uint_t t
        np.ndarray[real_t, ndim=1] dW = np.asarray(np.random.randn(N)*np.sqrt(dt), dtype = np.float32)
        np.ndarray[real_t, ndim=1] X = np.zeros(N, dtype = np.float32)
        np.ndarray[real_t, ndim=1] T = np.zeros(N, dtype = np.float32)
    
    X[0] = init
    
    for t in range(N - 1):
        X[t+1] = X[t] + f(X[t])*dt*pow(scale(X[t]), L) +  g(X[t])*dW[t]*sqrt(pow(scale(X[t]), L))
        T[t+1] = T[t] + dt*pow(scale(X[t]), L)
    
    return X, T

Pure Python (+Numpy)

#The numpy arrays will not improve the speed here much, since there's still a loop over the arrays...
def simulation(L = 0, N = 100000, dt = 1E-3, init = .1):
    """Simulate a stochastic differential equation.
    """

    #Set up some parameters:
    f1 = .1
    g1 = .01
    g2 = .1
    
    dW = np.random.randn(N)*np.sqrt(dt)
    X = np.zeros(N)
    T = np.zeros(N)
    X[0] = init
    
    def f(X):
        return -f1*X
    
    def g(X):
        return np.sqrt((g1 + g2*X**2))
    
    def scale(X):
        return 1/np.abs(g1 + g2*X**2)
    
    for t in xrange(0,N - 1):
        X[t+1] = X[t] + f(X[t])*dt*scale(X[t]) +  g(X[t])*dW[t]*np.sqrt(scale(X[t]))
        T[t+1] = T[t] + dt*scale(X[t])
    
    return X, T