initial commit

.gitignore

@@ -0,0 +1,50 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+*.manifest
+*.spec
+
+# PyBuilder
+target/
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.cache
+coverage.xml
+
+# Translations
+*.mo
+*.pot
+
+# Sphinx documentation
+_build
+
+*.log

MANIFEST.in

@@ -0,0 +1,11 @@
+include MANIFEST.in
+include *.txt
+include LICENSE
+include README.rst
+include tox.ini
+recursive-include sdeint *
+recursive-include examples *
+recursive-include docs *
+prune docs/_build
+prune */__pycache__
+global-exclude *.pyc *.pyo *.pyd *~ *.swp

README.rst

@@ -0,0 +1,52 @@
+sdeint
+======
+
+| Numerical integration of Ito or Stratonovich SDEs.
+
+Overview
+--------
+sdeint is a collection of numerical algorithms for integrating Ito and Stratonovich stochastic ordinary differential equations (SODEs). It has simple functions that can be used similarly to ``scipy.integrate.odeint()`` or MATLAB's ``ode45``.
+
+There already exist some python and MATLAB packages providing Euler-Maruyama and Milstein algorithms, and a couple of others. So why am I bothering to make another package?  
+
+It is because there has been 25 years of further research with better methods but for some reason I can't find any open source reference implementations. Not even for those methods published by Kloeden and Platen way back in 1992. So I will aim to put some improved methods here.
+
+This is prototype code in python, so not aiming for speed. Later can always rewrite these with loops in C when speed is needed.
+
+Bug reports are very welcome!
+
+functions
+---------
+
+| ``itoint(f, G, y0, tspan)`` for Ito equation dy = f dt + G dW
+| ``stratint(f, G, y0, tspan)`` for Stratonovich equation dy = f dt + G∘dW
+
+These will choose an algorithm for you. Or you can use a specific algorithm directly:
+
+specific algorithms:
+--------------------
+So far have these algorithms as a starting point.
+    ``itoEuler()``: the Euler Maruyama algorithm for Ito equations
+    ``stratHeun()``: the Stratonovich Heun algorithm for Stratonovich equations
+
+
+TODO
+----
+- Write tests (using systems that can be solved exactly)
+
+- Add more recent strong stochastic Runge-Kutta algorithms.
+  Perhaps starting with Rößler (2010) or Burrage and Burrage (1996)
+
+- Add an implicit strong algorithm useful for stiff equations, perhaps the one
+  from Kloeden and Platen (1999) section 12.4.
+
+- Currently prioritizing those algorithms that work for very general d-dimensional systems with arbitrary noise coefficient matrix, and which are derivative free. Eventually will add special case algorithms that give a speed increase for systems with certain symmetries. That is, 1-dimensional systems, systems with scalar noise, diagonal noise or commutative noise, etc.
+
+- Eventually implement the main loops in C for speed.
+
+- Some time in the dim future, implement support for stochastic delay differential equations (SDDEs).
+
+See also:
+---------
+
+``nsim``: Framework that uses this ``sdeint`` library to do massive parallel simulations of SDE systems (using multiple CPUs or a cluster) and provides some tools to analyze the resulting timeseries. https://github.com/mattja/nsim

requirements.txt

@@ -0,0 +1 @@
+numpy>=1.6

sdeint/__init__.py

@@ -0,0 +1,6 @@
+from __future__ import absolute_import
+
+from .wiener import dW, Ikp, Jkp, Iwiktorsson, Jwiktorsson
+from .integrate import itoint, stratint, itoEuler, stratHeun
+
+__version__ = '0.1.0-dev'

sdeint/integrate.py

@@ -0,0 +1,192 @@
+# Copyright 2015 Matthew J. Aburn
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version. See <http://www.gnu.org/licenses/>.
+
+"""
+Numerical integration algorithms for Ito and Stratonovich stochastic ordinary
+differential equations.
+
+Usage:
+    itoint(f, G, y0, tspan)  for Ito equation dy = f dt + G dW
+    stratint(f, G, y0, tspan)  for Stratonovich equation dy = f dt + G \circ dW
+
+    y0 is the initial value
+    tspan is an array of time values (currently these must be equally spaced)
+    f is the deterministic part of the system (can be a scalar or  dx1  vector)
+    G is the stochastic part of the system (can be a scalar or  d x m matrix)
+
+sdeint will choose an algorithm for you. Or you can choose one explicitly:
+
+Algorithms implemented so far:
+    itoEuler: the Euler Maruyama algorithm for Ito equations.
+    stratHeun: the Stratonovich Heun algorithm for Stratonovich equations.
+"""
+
+import numpy as np
+from .wiener import deltaW
+
+
+class Error(Exception):
+    pass
+
+
+class SDEValueError(Error):
+    """Thrown if integration arguments fail some basic sanity checks"""
+    pass
+
+
+def _check_args(f, G, y0, tspan):
+    """Do some common validations. Find dimension d, number of Wiener noises m.
+    """
+    if not np.isclose(min(np.diff(tspan)), max(np.diff(tspan))):
+        raise SDEValueError('Currently time steps must be equally spaced.')
+    # Be flexible to allow scalar equations. convert them to a 1D vector system
+    if isinstance(y0, numbers.Number):
+        if isinstance(y0, numbers.Integral):
+            numtype = np.float64
+        else:
+            numtype = type(y0)
+        y0_orig = y0
+        y0 = np.array([y0], dtype=numtype)
+        def make_vector_fn(fn):
+            def newfn(y, t):
+                return np.array([fn(y[0], t)], dtype=numtype)
+            newfn.__name__ = fn.__name__
+            return newfn
+        def make_matrix_fn(fn):
+            def newfn(y, t):
+                return np.array([[fn(y[0], t)]], dtype=numtype)
+            newfn.__name__ = fn.__name__
+            return newfn
+        if isinstance(f(y0_orig, 0.0), numbers.Number):
+            f = make_vector_fn(f)
+        if isinstance(G(y0_orig, 0.0), numbers.Number):
+            G = make_matrix_fn(G)
+    # determine dimension d of the system
+    d = len(y0)
+    if len(f(y0, tspan[0])) != d or len(G(y0, tspan[0])) != d:
+        raise SDEValueError('y0, f and G have incompatible shapes.')
+    # determine number of independent Wiener processes m
+    m = G(y0, tspan[0]).shape[1]
+    return (d, m, f, G, y0, tspan)
+
+
+def itoint(f, G, y0, tspan):
+    """ Numerically integrate Ito equation  dy = f dt + G dW
+    """
+    # In future versions we can automatically choose here the most suitable
+    # Ito algorithm based on properties of the system and noise.
+    (d, m, f, G, y0, tspan) = _check_args(f, G, y0, tspan)
+    chosenAlgorithm = itoEuler
+    return chosenAlgorithm(f, G, y0, tspan)
+
+
+def stratint(f, G, y0, tspan):
+    """ Numerically integrate Stratonovich equation  dy = f dt + G \circ dW
+    """
+    # In future versions we can automatically choose here the most suitable
+    # Stratonovich algorithm based on properties of the system and noise.
+    (d, m, f, G, y0, tspan) = _check_args(f, G, y0, tspan)
+    chosenAlgorithm = stratHeun
+    return chosenAlgorithm(f, G, y0, tspan)
+
+
+def itoEuler(f, G, y0, tspan):
+    """Use the Euler-Maruyama algorithm to integrate the Ito equation
+    dy = f(y,t)dt + G(y,t) dW(t)
+
+    where y is the d-dimensional state vector, f is a vector-valued function,
+    G is an d x m matrix-valued function giving the noise coefficients and
+    dW(t) = (dW_1, dW_2, ... dW_m) is a vector of independent Wiener increments
+
+    Args:
+      f: callable(y, t) returning (d,) array
+         Vector-valued function to define the deterministic part of the system
+      G: callable(y, t) returning (d,m) array
+         Matrix-valued function to define the noise coefficients of the system
+      y0: array of shape (d,) giving the initial state vector y(t==0)
+      tspan (array): The sequence of time points for which to solve for y.
+        These must be equally spaced, e.g. np.arange(0,10,0.005)
+        tspan[0] is the intial time corresponding to the initial state y0.
+
+    Returns:
+      y: array, with shape (len(tspan), len(y0))
+         With the initial value y0 in the first row
+
+    Raises:
+      SDEValueError
+
+    See also:
+      G. Maruyama (1955) Continuous Markov processes and stochastic equations
+      Kloeden and Platen (1999) Numerical Solution of Differential Equations
+    """
+    (d, m, f, G, y0, tspan) = _check_args(f, G, y0, tspan)
+    n = len(tspan)
+    dt = (tspan[n-1] - tspan[0])/(n - 1)
+    # allocate space for result
+    y = np.zeros((n, d), dtype=type(y0[0]))
+    # pre-generate all Wiener increments (for m independent Wiener processes):
+    dWs = deltaW(n - 1, m, dt)
+    y[0] = y0;
+    for i in range(1, n):
+        t1 = tspan[i - 1]
+        t2 = tspan[i]
+        y1 = y[i - 1]
+        dW = dWs[i - 1]
+        y[i] = y1 + f(y1, t1)*dt + G(y1, t1).dot(dW)
+    return y
+
+
+def stratHeun(f, G, y0, tspan):
+    """Use the Stratonovich Heun algorithm to integrate Stratonovich equation
+    dy = f(y,t)dt + G(y,t) \circ dW(t)
+
+    where y is the d-dimensional state vector, f is a vector-valued function,
+    G is an d x m matrix-valued function giving the noise coefficients and
+    dW(t) = (dW_1, dW_2, ... dW_m) is a vector of independent Wiener increments
+
+    Args:
+      f: callable(y, t) returning (d,) array
+         Vector-valued function to define the deterministic part of the system
+      G: callable(y, t) returning (d,m) array
+         Matrix-valued function to define the noise coefficients of the system
+      y0: array of shape (d,) giving the initial state vector y(t==0)
+      tspan (array): The sequence of time points for which to solve for y.
+        These must be equally spaced, e.g. np.arange(0,10,0.005)
+        tspan[0] is the intial time corresponding to the initial state y0.
+
+    Returns:
+      y: array, with shape (len(tspan), len(y0))
+         With the initial value y0 in the first row
+
+    Raises:
+      SDEValueError
+
+    See also:
+      W. Rumelin (1982) Numerical Treatment of Stochastic Differential
+         Equations
+      R. Mannella (2002) Integration of Stochastic Differential Equations
+         on a Computer
+      K. Burrage, P. M. Burrage and T. Tian (2004) Numerical methods for strong
+         solutions of stochastic differential equations: an overview
+    """
+    (d, m, f, G, y0, tspan) = _check_args(f, G, y0, tspan)
+    n = len(tspan)
+    dt = (tspan[n-1] - tspan[0])/(n - 1)
+    # allocate space for result
+    y = np.zeros((n, d), dtype=type(y0[0]))
+    # pre-generate all Wiener increments (for m independent Wiener processes):
+    dWs = deltaW(n - 1, m, dt)
+    y[0] = y0;
+    for i in range(1, n):
+        t1 = tspan[i - 1]
+        t2 = tspan[i]
+        y1 = y[i - 1]
+        dW = dWs[i - 1]
+        ybar = y1 + f(y1, t1)*dt + G(y1, t1).dot(dW)
+        y[i] = (y1 + 0.5*(f(y1, t1) + f(ybar, t2))*dt +
+                0.5*(G(y1, t1) + G(ybar, t2)).dot(dW))
+    return y

sdeint/wiener.py

@@ -0,0 +1,68 @@
+# Copyright 2015 Matthew J. Aburn
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version. See <http://www.gnu.org/licenses/>.
+
+"""
+Simulation of standard multiple stochastic integrals, both Ito and Stratonovich
+I_{ij}(t) = \int_{0}^{t}\int_{0}^{s} dW_i(u) dW_j(s)  (Ito)
+J_{ij}(t) = \int_{0}^{t}\int_{0}^{s} \circ dW_i(u) \circ dW_j(s)  (Stratonovich)
+
+These multiple integrals I and J are important building blocks that will be
+used by most of the higher-order algorithms that integrate multi-dimensional
+SODEs.
+
+We first implement the method of Kloeden, Platen and Wright (1992) to
+approximate the integrals by the first n terms from the series expansion of a
+Brownian bridge process. By default using n=5.
+
+Finally we implement the method of Wiktorsson (2001) which improves on the
+previous method by also approximating the tail-sum distribution by a
+multivariate normal distribution.
+"""
+
+import numpy as np
+
+
+def deltaW(n, m, delta_t):
+    """pre-generate sequence of Wiener increments for each of m independent
+    Wiener processes W_j(t) j=1..m over n time steps with constant
+    time step size delta_t.
+    
+    Returns:
+      array of shape (n, m), where the [k, j] element has the value
+        W_j((k+1)delta_t) - W_j(k*delta_t)
+    """
+    sqrth = np.sqrt(h)
+    return np.random.normal(0.0, sqrth, (n, m))
+
+
+def _kprod(A, B):
+    """Kroenecker tensor product of two matrices.
+    If A is m x n matrix and B is p x q then _kprod(A,B) is a mp x nq matrix.
+    """
+    pass
+
+
+def _vec(A):
+    """Stack columns of matrix A on top of each other to give a long vector.
+    If A is m x n matrix then _vec(A) will be mn x 1
+    """
+
+
+def Ikp(h, m, n=5):
+    pass
+
+
+def Jkp(h, m, n=5):
+    pass
+
+
+def Iwiktorsson(h, m, n=5):
+    pass
+
+
+def Jwiktorsson(h, m, n=5):
+    pass