Pure Python

gistfile1.py

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