CODEWITHSUNDEEP
pythonnumpystatistics
physicsGuy
physicsGuy

Reputation: 3675

Stochastic integration with python

I want to numerically solve integrals that contain white noise.

Mathematically white noise can be described by a variable X(t), which is a random variable with a time average, Avg[X(t)] = 0 and the correlation function, Avg[X(t), X(t')] = delta_distribution(t-t').

A simple example would be to calculate the integral over X(t) from t=0 to t=1. On average this is of course zero, but what I need are different realizations of this integral.

The problem is that this does not work with numpy.integrate.quad().

Are there any packages for python that deal with stochastic integrals?

Upvotes: 8

Views: 3978

Answers (1)

overfull hbox
overfull hbox

Reputation: 861

This is a good starting point for numerical SDE methods: http://math.gmu.edu/~tsauer/pre/sde.pdf.

Here is a simple numpy solver for the stochastic differential equation dX_t = a(t,X_t)dt + b(t,X_t)dW_t which I wrote for a class project last year. It is based on the forward euler method for regular differential equations, and in practice is fairly widely used when solving SDEs.

def euler_maruyama(a,b,x0,t):
    N = len(t)
    x = np.zeros((N,len(x0)))
    x[0] = x0
    for i in range(N-1):
        dt = t[i+1]-t[i]
        dWt = np.random.normal(0,dt)
        x[i+1] = x[i] + a(t[i],x[i])*dt + b(t[i],x[i])*dWt
    return x

Essentially, at each timestep, the deterministic part of the function is integrated using forward Euler, and the stochastic part is integrated by generating a normal random variable dWt with mean 0 and variance dt and integrating the stochastic part with respect to this.

The reason we generate dWt like this is based on the definition of Brownian motions. In particular, if $W$ is a Brownian motion, then $(W_t-W_s)$ is normally distributed with mean 0 and variance $t-s$. So dWt is a discritization of the change in $W$ over a small time interval.

This is a the docstring from the function above:

Parameters
----------
a : callable a(t,X_t),
    t is scalar time and X_t is vector position
b : callable b(t,X_t),
    where t is scalar time and X_t is vector position
x0 : ndarray
    the initial position
t : ndarray
    list of times at which to evaluate trajectory

Returns
-------
x : ndarray
    positions of trajectory at each time in t

Upvotes: 4

Related Questions