The SdePy Package¶

Getting Started¶

SdePy¶

The SdePy package provides tools to state and numerically integrate Ito Stochastic Differential Equations (SDEs), including equations with time-dependent parameters, time-dependent correlations, and stochastic jumps, and to compute with, and extract statistics from, their realized paths.

Several preset processes are provided, including lognormal, Ornstein-Uhlenbeck, Hull-White n-factor, Heston, and jump-diffusion processes.

Computations are fully vectorized across paths, via NumPy and SciPy, making live sessions with 100000 paths reasonably fluent on single cpu hardware.

This package came out of practical need, so expect a flexible tool that gets real-life things done. On the other hand, not every part of it is clean and polished, so expect rough edges, and the occasional bug (please report!).

Developers are committed to the stability of the public API, here again out of practical need to safeguard dependencies.

Start here¶

License¶

BSD 3-Clause License

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

This package reuses the compatibly licensed files listed below.

File: sdepy/doc/_templates/autosummary/class.rst License: 3-clause BSD

For details, see sdepy/doc/_templates/autosummary/LICENSE.txt

Quick Guide¶

Install and import¶

Install using pip install sdepy, or copy the package source code in a directory in your Python path.

Import as

>>> import sdepy
>>> import numpy as np
>>> import matplotlib.pyplot as plt  # optional, if plots are needed
>>> plt.rcParams['figure.figsize'] = (11., 5.5)
>>> plt.rcParams['lines.linewidth'] = 1.

How to state an SDE¶

Here follows a bare-bone definition of a Stochastic Differential Equation (SDE), in this case a Ornstein-Uhlenbeck process:

>>> @sdepy.integrate
... def my_process(t, x, theta=1., k=1., sigma=1.):
...     return {'dt': k*(theta - x), 'dw': sigma}

This represents the SDE dX = k*(theta - X)*dt + sigma*dW(t), where theta, k and sigma are parameters and dW(t) are Wiener process increments. A further 'dn' or 'dj' entry in the returned dictionary would allow for Poisson or compound Poisson jumps.

A number of preset processes are provided, including lognormal processes, Hull-White n-factor processes, Heston processes, and jump-diffusion processes.

How to integrate an SDE¶

Now my_process is a class, a subclass of the cooperating sdepy.SDE and sdepy.integrator classes:

>>> issubclass(my_process, sdepy.integrator), issubclass(my_process, sdepy.SDE)
(True, True)

It is to be instantiated with a number of parameters, including the SDE parameters theta, k and sigma; its instances are callable, given a timeline they will integrate and return the process along it. Decorating my_process with sdepy.kunfc allows for more concise handling of parameters:

>>> myp = sdepy.kfunc(my_process)

It is best explained by examples, involving my_process, myp and

>>> coarse_timeline = (0., 0.25, 0.5, 0.75, 1.0)
>>> timeline = np.linspace(0., 1., 500)

Scalar process in 100000 paths, with default parameters, computed at 5 time points (coarse_timeline), using 100 steps in between:

>>> x = my_process(x0=1, paths=100*1000,
...                steps=100)(coarse_timeline)
>>> x.shape
(5, 100000)

The same scalar process computed on a fine-grained timeline (timeline) and 1000 paths, using one integration step for each point in the timeline (no steps parameter):
```
>>> x = my_process(x0=1, paths=1000,
...                steps=100)(timeline)
>>> x.shape
(500, 1000)
```
A plot of a few paths may be used to inspect the integration result:
```
>>> gr = plt.plot(timeline, x[:, :30])
>>> plt.show()  
```

Vector process with three components and correlated Wiener increments (same other parameters as above):

>>> corr = ((1, .2, -.3), (.2, 1, .1), (-.3, .1, 1))
>>> x = my_process(x0=1, vshape=3, corr=corr,
...                paths=1000)(timeline)
>>> x.shape
(500, 3, 1000)

Vector process as above, with 10000 paths and time-dependent parameters and correlations:

>>> sigma = lambda t: 0.1 + t
>>> theta = lambda t: 2-t
>>> k = lambda t: 2/(t+1)
>>> c02 = lambda t: -0.1*np.cos(3*t)
>>> c12 = lambda t: 0.1*np.sign(0.5-t)
>>> corr = lambda t: ((     1,    -.2, c02(t)),
...                   (   -.2,      1, c12(t)),
...                   (c02(t), c12(t),      1))
>>> x = my_process(x0=1, vshape=3, corr=corr,
...                theta=theta, k=k, sigma=sigma, paths=10*1000)(timeline)
>>> x.shape
(500, 3, 10000)

This plot illustrates the correlations among the components of x increments, as a function of time and as compared to corr(t):

>>> dx = np.diff(x, axis=0)
>>> for i in range(3):
...     for j in range(i + 1, 3):
...         gr = plt.plot(
...             timeline, corr(timeline)[i][j] + 0*timeline,
...             timeline[1:], [np.cov(z)[i, j]/(z[i].std()*z[j].std())
...                            for z in dx]
...             )
>>> plt.show()  

A 1000 paths scalar process with path-dependent initial conditions and parameters, integrated backwards (i0=-1):

>>> x0, sigma = np.zeros(1000), np.zeros(1000)
>>> x0[::2], x0[1::2] = 0., 2.
>>> sigma[::2], sigma[1::2] = 0.5, 0.1
>>> x = my_process(x0=x0, sigma=sigma, paths=1000,
...                theta=1, k=-2,
...                i0=-1)(timeline)
>>> x.shape
(500, 1000)

When integrating backwards, the inital conditions are applied at the final point in the given timeline:

>>> assert (x[-1, :] == x0).all()
>>> gr = plt.plot(timeline, x[:, :30])
>>> gr = plt.plot(timeline, np.full_like(timeline, 1), 'k--')
>>> plt.show()  

Note the negative value of k, with mean reversion towards theta=1 occurring backwards in time.

A scalar process computed on a 10 x 15 grid of parameters sigma and k (note that the shape of the initial conditions and of each parameter should be broadcastable to the values of the process across paths, i.e. to shape vshape + (paths,)):

>>> sigma = np.linspace(0., 1., 10).reshape(10, 1, 1)
>>> k = np.linspace(1., 2., 15).reshape(1, 15, 1)
>>> x = my_process(x0=1, theta=2, k=k, sigma=sigma, vshape=(10, 15),
...                paths=10*1000)(coarse_timeline)
>>> x.shape
(5, 10, 15, 10000)

A plot of the final average process values against k illustrates a faster reversion to theta=2 as k increases, as well as the independence of the process mean from sigma:

>>> for i in range(10):
...     gr = plt.plot(k[0, :, 0], x[-1, i, :, :].mean(axis=-1))
>>> lb = plt.xlabel('k'), plt.ylabel('x(t=2).mean()')
>>> plt.show()  

In the example above, set steps>=100 to go from inaccurate and fast, to meaningful and slow.

Interactive modification of process and integration parameters using the sdepy.kfunc decorator myp = sdepy.kfunc(my_process).

The sdepy.kfunc decorated version of my_process is a subclass of sdepy.integrator and sdepy.SDE, as my_process is, and fully replicates its functionality and interface:
```
>>> issubclass(myp, sdepy.integrator), issubclass(myp, sdepy.SDE)
(True, True)
```
In addition, and in contrast to my_process, myp instances accept either an integration timeline, or a modified value of some integration or SDE parameters, or both, as illustrated below:
```
>>> p = myp(x0=1, sigma=1, paths=1000)

>>> x = p(timeline)
>>> x1, x2 = p(timeline, sigma=0.5), p(timeline, sigma=1.5)
>>> q = p(paths=100, vshape=(3,), k=2)
>>> y = q(timeline, sigma=0.5)
```
x is the result of integrating p along timeline (no difference here from a my_process instance); x1, x2 are obtained by integration along timeline by setting sigma to the given values, and keeping other parameters as stated when p was instantiated; q is another myp instance with updated default values for paths, vshape and k, and all else set as in p; and finally, y was obtained by integrating q along timeline, with its own parameters, save for sigma that was modified to 0.5.

Moreover, for sdepy.kfunc classes, instantiation and computation may happen contextually:
```
>>> x = myp(timeline, x0=1, sigma=1, paths=1000)
```
is equivalent to:
```
>>> x = my_process(x0=1, sigma=1, paths=1000)(timeline)
```
sdepy.kfunc-decorated classes allow to state some central values of parameters for a given problem, and to explore the effects of variatons in some of them via a concise interface, that keeps the modified parameters in focus and all the rest in the background.

To inspect the parameters stored in a sdepy.kfunc instance, use the read-only params attribute:
```
>>> q.params  
{
   'paths': 100,
   'vshape': (3,),
   'x0': array(1),
   'sigma': array(1),
   'k': array(2),
   ...,
}
```
To test if an object is a kfunc, use sdepy.iskfunc():
```
>>> sdepy.iskfunc(myp), sdepy.iskfunc(p), sdepy.iskfunc(my_process)
(True, True, False)
```
The examples that follow illustrate, among other things, the use of myp as a sdepy.kfunc class.

Processes generated using integration results as stochasticity sources (mind using consistent vshape and paths, and synchronizing timelines):

>>> my_dw = sdepy.integrate(lambda t, x: {'dw': 1})(vshape=1, paths=1000)(timeline)
>>> p = myp(dw=my_dw, vshape=3, paths=1000,
...         x0=1, sigma=((1,), (2,), (3,)))
>>> x = p(timeline)
>>> x.shape
(500, 3, 1000)

Now, x1, x2, x3 = = x[:, 0], x[:, 1], x[:, 2] have different sigma, but share the same dw increments, as can be seen plotting a path:

>>> k = 0  # path to be plotted
>>> gr = plt.plot(timeline, x[:, :, k])
>>> plt.show()  

If more integrations steps are needed between points in the output timeline, use steps to keep the integration timeline consistent with the one of my_dw:

>>> x = p(coarse_timeline, steps=timeline)
>>> x.shape
(5, 3, 1000)

Using stochasticity sources with memory (mind using consistent vshape and paths):

>>> my_dw = sdepy.true_wiener_source(paths=1000)
>>> p = myp(x0=1, theta=1, k=1, sigma=1, dw=my_dw, paths=1000)

my_dw, as a sdepy.true_wiener_source instance has memory of, and generates new Wiener process increments consistent with, its formerly realized values. As a consequence, processes defined invoking p share the same underlying Wiener process increments:

>>> t1 = np.linspace(0., 1.,  30)
>>> t2 = np.linspace(0., 1., 100)
>>> t3 = t = np.linspace(0., 1., 300)
>>> x1, x2, x3 = p(t1), p(t2), p(t3)
>>> y1, y2, y3 = p(t, theta=1.5), p(t, theta=1.75), p(t, theta=2)

x1, x2, x3 illustrate SDE integration convergence as time steps become smaller, and y1, y2, y3 illustrate how theta affects paths, all else being equal:

>>> i = 0 # path to be plotted
>>> gr = plt.plot(t, x1(t)[:, i], t, x2(t)[:, i], t, x3(t)[:, i])
>>> plt.show() 
>>> gr = plt.plot(t, y1[:, i], t, y2[:, i], t, y3[:, i])
>>> plt.show()  

How to handle the integration output¶

SDE integrators return instances of sdepy.process, a subclass of np.ndarray with a timeline stored in the t attribute (note the shape of x, repeatedly used in the examples below):

>>> coarse_timeline = (0., 0.25, 0.5, 0.75, 1.0)
>>> timeline = np.linspace(0., 1., 101)
>>> x = my_process(x0=1, vshape=3, paths=1000)(timeline)
>>> x.shape
(101, 3, 1000)

x is a sdepy.process instance:

>>> type(x)
<class 'sdepy.infrastructure.process'>

and is based on the given timeline:

>>> np.isclose(timeline, x.t).all()
True

Whenever possible, a process will store references, not copies, of timeline and values. In fact:

>>> timeline is x.t
True

The first axis is reserved for the timeline, the last for paths, and axes in the middle match the shape of process values:

>>> x.shape == x.t.shape + x.vshape + (x.paths,)
True

Calling processes interpolates in time:

>>> y = x(coarse_timeline)
>>> y.shape
(5, 3, 1000)

The result is always an array, not a process:

>>> type(y)
<class 'numpy.ndarray'>

Indexing works as usual, and returns NumPy arrays:

>>> type(x[0])
<class 'numpy.ndarray'>

All array methods are unchanged (no overriding), and return NumPy arrays as well:

>>> type(x.mean(axis=0))
<class 'numpy.ndarray'>

You can slice processes along time, values and paths with special indexing.

Time indexing:

>>> y = x['t', ::2]
>>> y.shape
(51, 3, 1000)

Values indexing:

>>> y = x['v', 0]
>>> y.shape
(101, 1000)

Paths indexing:

>>> y = x['p', :10]
>>> y.shape
(101, 3, 10)

The output of a special indexing operation is a process:

>>> isinstance(y, sdepy.process)
True

Smart indexing is allowed. To select paths that cross x=0 at some point and for some component, use:

>>> i_negative = x.min(axis=(0, 1)) < 0
>>> y = x['p', i_negative]
>>> y.shape == (101, 3, i_negative.sum())
True

You can do algebra with processes that either share the same timeline, or are constant (a process with a one-point timeline is assumed to be constant), and either have the same number of paths, or are deterministic (with one path):

>>> x_const = x['t', 0]  # a constant process
>>> x_one_path = x['p', 0]  # a process with one path

>>> y = np.exp(x) - x_const
>>> z = np.maximum(x, x_one_path)

>>> isinstance(y, sdepy.process), isinstance(z, sdepy.process)
(True, True)

When integrating SDEs, the SDE parameters and/or stochasticity sources accept processes as valid values (mind using deterministic processes, or synchronizing the number of paths, and make sure that the shape of values do broadcast together). To use a realization of my_process as the volatility of a 3-component lognormal process, do as follows:

>>> stochastic_vol = my_process(x0=1, paths=10*1000)(timeline)
>>> stochastic_vol_x = sdepy.lognorm_process(x0=1, vshape=3, paths=10*1000,
...     mu=0, sigma=stochastic_vol)(timeline)

Processes have specialized methods, and may be analyzed, and their statistics cumulated across multiple runs, using the sdepy.montecarlo class. Some examples follow:

Cumulative probability distribution function at t=0.5 of the process values of x across paths:
```
>>> cdf = x.cdf(0.5, x=np.linspace(-2, 2, 100))  # an array
```

Characteristic function at t=0.5 of the same distribution:

>>> chf = x.chf(0.5, u=np.linspace(-2, 2, 100))  # an array

Standard deviation across paths:

>>> std = x.pstd()  # a one-path process
>>> std.shape
(101, 3, 1)

Maximum value reached along the timeline:

>>> xmax = x.tmax()  # a constant process
>>> xmax.shape
(1, 3, 1000)

A linearly interpolated, or Gaussian kernel estimate (default) of the probability distribution function (pdf) and its cumulated values (cdf) across paths, at a given time point, may be obtained using the montecarlo class:

>>> y = x(1)[0]  # 0-th component of x at time t=1
>>> a = sdepy.montecarlo(y, bins=30)
>>> ygrid = np.linspace(y.min(), y.max(), 200)
>>> gr = plt.plot(ygrid, a.pdf(ygrid), ygrid, a.cdf(ygrid))
>>> gr = plt.plot(ygrid, a.pdf(ygrid, method='interp', kind='nearest'))
>>> plt.show()  # doctest: +SKIP

A sdepy.montecarlo instance can be used to cumulate the results of multiple simulations, across multiple components of process values:

>>> p = my_process(x0=1, vshape=3, paths=10*1000)
>>> a = sdepy.montecarlo(bins=100)  # empty montecarlo instance
>>> for _ in range(10):
...     x = p(timeline)  # run simulation
...     a.update(x(1))  # cumulate x values at t=1
>>> a.paths
100000
>>> gr = plt.plot(ygrid, a[0].pdf(ygrid), ygrid, a[0].cdf(ygrid))
>>> gr = plt.plot(ygrid, a[0].pdf(ygrid, method='interp', kind='nearest'))
>>> plt.show()  

Example - Stochastic Runge-Kutta¶

Minimal implementation of a basic stochastic Runge-Kutta integration scheme, as a subclass of sdepy.integrator (the A and dZ methods below are the standardized way in which equations are exposed to integrators):

>>> from numpy import sqrt
>>> class my_integrator(sdepy.integrator):
...     def next(self):
...         t, new_t = self.itervars['sw']
...         x, new_x = self.itervars['xw']
...         dt = new_t - t
...         A, dZ = self.A(t, x), self.dZ(t, dt)
...         a, b, dw = A['dt'], A['dw'], dZ['dw']
...         b1 = self.A(t, x + a*dt + b*sqrt(dt))['dw']
...         new_x[...] = x + a*dt + b*dw + (b1 - b)/2 * (dw**2 - dt)/sqrt(dt)

SDE of a lognormal process, as a subclass of sdepy.SDE, and classes that integrate it with the default integration method (euler) and via my_integrator (rk):

>>> class my_SDE(sdepy.SDE):
...     def sde(self, t, x):
...         return {'dt': 0, 'dw': x}

>>> class euler(my_SDE, sdepy.integrator):
...     pass

>>> class rk(my_SDE, my_integrator):
...     pass

Comparison of integration errors, as the integration from t=0 to t=1 is carried out with an increasing number of steps, against the integration result of sdepy.lognorm_process, which returns an exact result irrespective of the number and size of the integration steps (this happens since, by implementation, it integrates the linear SDE for log(x)):

>>> args = dict(dw=sdepy.true_wiener_source(paths=100),
...             paths=100, x0=10)
>>> timeline = (0, 1)
>>> steps = np.array((2, 3, 5, 10, 20, 30, 50, 100,
...                   200, 300, 500, 1000, 2000, 3000))

>>> # exact integration results at t=1
>>> exact = sdepy.lognorm_process(mu=0, sigma=1, **args)(timeline)[-1].mean()

>>> # errors of approximate integration results at t=1
>>> errors = np.abs(np.array([
...     [euler(**args, steps=s)(timeline)[-1].mean()/exact - 1,
...      rk(**args, steps=s)(timeline)[-1].mean()/exact - 1]
...     for s in steps]))

>>> # plots
>>> ax = plt.axes(label=0); ax.set_xscale('log'); ax.set_yscale('log')
>>> gr = ax.plot(steps, errors)
>>> plt.show()  
>>> print(f'euler error: {errors[-1,0]:.2e}\n'
...       f'   rk error: {errors[-1,1]:.2e}')  
euler error: 1.70e-03
   rk error: 8.80e-06

Example - Fokker-Planck Equation¶

Monte Carlo integration of partial differential equations, illustrated in the simplest example of the heat equation diff(u, t) - k*diff(u, x, 2) == 0, for the function u(x, t), i.e. the Fokker-Planck equation for the SDE dX(t) = sqrt(2*k)*dW(t). Initial conditions at t=t0, two examples:

u(x, t0) = 1 for lb < x < hb and 0 otherwise,
u(x, t0) = sin(x).

Setup:

>>> from numpy import exp, sin
>>> from scipy.special import erf
>>> from scipy.integrate import quad

>>> k = .5
>>> x0, x1 = 0, 10;
>>> t0, t1 = 0, 1
>>> lb, hb = 4, 6

Exact green function and solutions for initial conditions 1. and 2., to be checked against results:

>>> def green_exact(y, s, x, t):
...     return exp(-(x - y)**2/(4*k*(t - s)))/sqrt(4*np.pi*k*(t - s))

>>> def u1_exact(x, t):
...     return (erf((x - lb)/2/sqrt(k*(t - t0))) - erf((x - hb)/2/sqrt(k*(t - t0))))/2

>>> def u2_exact(x, t):
...     return exp(-k*(t - t0))*sin(x)

Realization of the needed stochastic process, by backward integration from a grid of final values of x at t=t1, using the preset wiener_process class (the steps keyword is added as a reminder of the setup needed for less-than-trivial equations, it does not actually make a difference here):

>>> xgrid = np.linspace(x0, x1, 51)
>>> tgrid = np.linspace(t0, t1, 5)
>>> xp = sdepy.wiener_process(
...     paths=10000, steps=100,
...     sigma=sqrt(2*k),
...     vshape=xgrid.shape, x0=xgrid[..., np.newaxis],
...     i0=-1,
...     )(timeline=tgrid)

Computation of the green function and of the solutions u(x, t1) via Monte Carlo integration (note the liberal use of scipy.integrate.quad below, enabled by the smoothness of the Gaussian kernel estimate a[i, j].pdf):

>>> a = sdepy.montecarlo(xp, bins=100)

>>> def green(y, i, j):
...     """green function from (y=y, s=tgrid[i]) to (x=xgrid[j], t=t1)"""
...     return a[i, j].pdf(y)

>>> u1, u2 = np.empty(51), np.empty(51)
>>> for j in range(51):
...     u1[j] = quad(lambda y: green(y, 0, j), lb, hb)[0]
...     u2[j] = quad(lambda y: sin(y)*green(y, 0, j), -np.inf, np.inf)[0]

Comparison against exact values:

>>> y = np.linspace(x0, x1, 500)
>>> for i, j in ((1, 20), (2, 30), (3, 40)):
...     gr = plt.plot(y, green(y, i, j),
...                   y, green_exact(y, tgrid[i], xgrid[j], t1), ':')
>>> plt.show()  

>>> gr = plt.plot(xgrid, u1, y, u1_exact(y, t1), ':')
>>> gr = plt.plot(xgrid, u2, y, u2_exact(y, t1), ':')
>>> plt.show()  

>>> print(f'u1 error: {np.abs(u1 - u1_exact(xgrid, t1)).mean():.2e}\n'
...       f'u2 error: {np.abs(u2 - u2_exact(xgrid, t1)).mean():.2e}'
... )  
u1 error: 2.49e-03
u2 error: 5.51e-03

Example - Basket Lookback Option¶

Take a basket of 4 financial securities, with risk-neutral probabilities following lognormal processes in the Black-Scholes framework. Correlations, dividend yields and term structure of volatility (will be linearly interpolated) are given below:

>>> corr = [
...     [1,    0.50, 0.37, 0.35],
...     [0.50,    1, 0.47, 0.46],
...     [0.37, 0.47,    1, 0.19],
...     [0.35, 0.46,  0.19,   1]]

>>> dividend_yield = sdepy.process(c=(0.20, 4.40, 0., 4.80))/100
>>> riskfree = 0  # to keep it simple

>>> vol_timepoints = (0.1, 0.2, 0.5, 1, 2, 3)
>>> vol = np.array([
...     [0.40, 0.38, 0.30, 0.28, 0.27, 0.27],
...     [0.31, 0.29, 0.22, 0.16, 0.18, 0.21],
...     [0.24, 0.22, 0.19, 0.19, 0.21, 0.22],
...     [0.35, 0.31, 0.21, 0.18, 0.19, 0.19]])
>>> sigma = sdepy.process(t=vol_timepoints, v=vol.T)
>>> sigma.shape
(6, 4, 1)

The prices of the securities at the end of each quarter for the next 2 years, simulated across 50000 independent paths and their antithetics (sdepy.odd_wiener_source is used), are:

>>> maturity = 2
>>> timeline = np.linspace(0, maturity, 4*maturity + 1)
>>> p = sdepy.lognorm_process(
...     x0=100, corr=corr, dw=sdepy.odd_wiener_source,
...     mu=(riskfree - dividend_yield),
...     sigma=sigma,
...     vshape=4, paths=100*1000, steps=maturity*250)
>>> x = p(timeline)
>>> x.shape
(9, 4, 100000)

A call option knocks in if any of the securities reaches a price below 80 at any quarter (starting from 100), and pays the lookback maximum attained by the basket (equally weighted), minus 105, if positive. Its price is:

>>> x_worst = x.min(axis=1)
>>> x_basket = x.mean(axis=1)
>>> down_and_in_paths = (x_worst.min(axis=0) < 80)
>>> lookback_x_basket = x_basket.max(axis=0)
>>> payoff = np.maximum(0, lookback_x_basket - 105)
>>> payoff[np.logical_not(down_and_in_paths)] = 0
>>> a = sdepy.montecarlo(payoff, use='even')
>>> print(a)  
4.997 +/- 0.027

Release Notes¶

This section reports notable changes and additions in SdePy releases.

Not all releases are covered. For further information, and a complete list of releases, please refer to the SdePy repository on GitHub.

SdePy 1.2.0¶

New in this release

SdePy was upgraded to the current NumPy pseudo-random numbers generation framework, based on the PCG-64 algorithm. The following additions were made:

SDEs and stochasticity source classes now accept upon instantiation a new optional parameter rng, to locally set the random number generator used by each instance. Such rng is expected to expose the interface used by numpy.random.Generator and numpy.random.RandomState instances.
In case no rng parameter is given, random number generation is delegated to a global numpy.random.default_rng() object, created upon import: absent any user intervention, unmodified pre-existing code will rely on PCG-64 random number generation after upgrading to current NumPy and SdePy releases.
To access the stably reproducible NumPy legacy random generation, it is recommended to instantiate sources and SDEs with the rng=numpy.random.RandomState(SEED) parameter. As a result, each object is served the random numbers stream formerly obtained from the global NumPy random state after a numpy.random.seed(SEED) call.
Compatibility with legacy NumPy versions was maintained, if now deprecated: in such case, SdePy silently falls back on the legacy global NumPy random state.

The present SdePy version safeguards functional compatibility with code written for previous versions, note however that, if such code is run with current NumPy and SdePy versions, it will break pseudo-random numbers reproducibility:

numpy.random.seed calls no longer affect SdePy objects.
To reproduce the default output of former SdePy versions, numpy.random.seed(SEED) statements may be replaced with the following assignment:
```
sdepy.infrastructure.default_rng = numpy.random.RandomState(SEED)
```
Any other use of such global variable should rather be avoided, in favor of the rng keyword.

Improvements

The testing suite was updated, and its interface sdepy.test() extended, for use beyond NumPy legacy random generation.
GitHub Actions were updated accordingly, to perform CI tests using both legacy and current NumPy random generation.

Changes

Python 3.5 is no longer supported.

SdePy 1.1.0¶

New in this release

The sdepy.process class acquired new methods vmin, vmax, vmean, vvar, vstd to perform summary operations across values, for each time point and path.
A piecewise constant process constructor was added as the piecewise function; it replaces, and improves upon, the private _piecewise_constant_process (undocumented, not part of the API), now deprecated.

Bug-fixes

An incompatibility issue of the sdepy.process class with NumPy versions >= 1.16.0 was solved.
A bug in the out parameter of process methods pmin, pmax, tmin, tmax was fixed. No backward incompatible changes were made against the stated API, note however that code relying on the out parameter of the mentioned process methods might need fixing as well.

API Documentation¶

Overview¶

This package provides tools to state and numerically integrate Ito Stochastic Differential Equations (SDEs), including equations with time-dependent parameters, time-dependent correlations, and stochastic jumps, and to compute with, and extract statistics from, their realized paths.

Package contents:

A set of tools to ease computations with stochastic processes, as obtained from numerical integration of the corresponding SDE, is provided via the process and montecarlo classes (see Infrastructure):

The process class, a subclass of numpy.ndarray representing a sequence of values in time, realized in one or several paths. Algebraic manipulations and ufunc computations are supported for instances that share the same timeline, or are constant, and comply with numpy broadcasting rules. Interpolation along the timeline is supported via callability of process instances. Process-specific functionalities, such as averaging and indexing along time or across paths, are delegated to process-specific methods, attributes and properties (no overriding of numpy.ndarray operations).

The montecarlo class, as an aid to cumulate the results of several Monte Carlo simulations of a given stochastic variable, and to extract summary estimates for its probability distribution function and statistics.

Numerical realizations of the differentials commonly found as stochasticity sources in SDEs, are provided via the source class and its subclasses, with or without memory of formerly invoked realizations (see Stochasticity Sources).

A general framework for stochastic step by step simulations, and for numerical SDE integration, is provided via the paths_generator class, and its cooperating subclasses integrator, SDE and SDEs (see SDE Integration Framework). The full API allows for extensive customization of preprocessing, post-processing, stochasticity sources instantiation and handling, integration algorithms etc. The integrate decorator provides a simple and concise interface to handle standard use cases, via Euler-Maruyama integration.

Several preset stochastic processes are provided, including lognormal, Ornstein-Uhlenbeck, Hull-White n-factor, Heston, and jump-diffusion processes (see Stochastic Processes). Each process consists of a process generator class, a subclass of integrator and SDE, named with a _process suffix, and a definition of the underlying SDE, a subclass of SDE or SDEs, named with a _SDE suffix.

Several analytical results relating to the preset stochastic processes are made available, as a general reference and for testing purposes (see Analytical Results). They are limited to the case of constant process parameters, and with some further limitations on the parameters’ domains. Function arguments are consistent with those of the corresponding processes. Suffixes _pdf, _cdf and _chf stand respectively for probability distribution function, cumulative probability distribution function, and characteristic function. Black-Scholes formulae for the valuation of call and put options have been included (with prefix bs).

As an aid to interactive and notebook sessions, shortcuts are provided for stochasticity sources and preset processes (see Shortcuts). Shortcuts have been wrapped as “kfuncs”, objects with managed keyword arguments that simplify interactive workflow when frequent parameters tuning operations are needed (see kfunc decorator documentation). Analytical results are wrapped as kfuncs as well.

For all sources and processes, values can take any shape, scalar or multidimensional. Correlated multivariate stochasticity sources are supported. Poisson jumps are supported, and may be compounded with any random variable supported by scipy.stats. Time-varying process parameters (correlations, intensity of Poisson processes, volatilities etc.) are allowed whenever applicable. process instances act as valid stochasticity source realizations (as does any callable object complying with a source protocol), and may be passed as a source specification when computing the realization of a given process.

Computations are fully vectorized across paths, providing an efficient infrastructure for simulating a large number of process realizations. Less so, for large number of time steps: integrating 100 time steps across one million paths takes seconds, one million time steps across 100 paths takes minutes.

Infrastructure¶

`process`([t, x, v, c, dtype])	Array representation of a process (a subclass of numpy.ndarray).
`piecewise`([t, x, v, dtype, mode])	Return a process that interpolates to a piecewise constant function.
`montecarlo`([sample, axis, bins, range, use, …])	Summary statistics of Monte Carlo simulations.

Stochasticity Sources¶

`source`(*[, paths, vshape, dtype, rng])	Base class for stochasticity sources.
`wiener_source`(*[, paths, vshape, dtype, …])	dw, a source of standard Wiener process (Brownian motion) increments.
`poisson_source`(*[, paths, vshape, dtype, …])	dn, a source of Poisson process increments.
`cpoisson_source`(*[, paths, vshape, dtype, …])	dj, a source of compound Poisson process increments (jumps).
`odd_wiener_source`(*[, paths, vshape, dtype, rng])	dw, a source of standard Wiener process (Brownian motion) increments with antithetic paths exposing opposite increments (averages exactly to 0 across paths).
`even_poisson_source`(*[, paths, vshape, …])	dn, a source of Poisson process increments with antithetic paths exposing identical increments.
`even_cpoisson_source`(*[, paths, vshape, …])	dj, a source of compound Poisson process increments (jumps) with antithetic paths exposing identical increments.
`true_source`(*[, paths, vshape, dtype, rng, …])	Base class for stochasticity sources with memory.
`true_wiener_source`(*[, paths, vshape, …])	dw, source of standard Wiener process (brownian motion) increments with memory.
`true_poisson_source`(*[, paths, vshape, …])	dn, a source of Poisson process increments with memory.
`true_cpoisson_source`(*[, paths, vshape, …])	dj, a source of compound Poisson process increments (jumps) with memory.
`norm_rv`([a, b])	Normal distribution with mean a and standard deviation b, possibly time-dependent.
`uniform_rv`([a, b])	Uniform distribution between a and b, possibly time-dependent.
`exp_rv`([a])	Exponential distribution with scale a, possibly time-dependent.
`double_exp_rv`([a, b, pa])	Double exponential distribution, with scale a with probability pa, and -b with probability (1 - pa), possibly time-dependent.
`rvmap`(f, y)	Map f to random variates of distribution y, possibly time-dependent.

SDE Integration Framework¶

`paths_generator`(*[, paths, xshape, wshape, …])	Step by step generation of stochastic simulations across multiple paths, intended for subclassing.
`integrator`(*[, paths, xshape, wshape, …])	Step by step numerical integration of Ito Stochastic Differential Equations (SDEs), intended for subclassing.
`SDE`(*[, paths, vshape, dtype, rng, steps, …])	Class representation of a user defined Stochastic Differential Equation (SDE), intended for subclassing.
`SDEs`(*[, paths, vshape, dtype, rng, steps, …])	Class representation of a user defined system of Stochastic Differential Equations (SDEs), intended for subclassing.
`integrate`([sde, q, sources, log, addaxis])	Decorator for Ito Stochastic Differential Equation (SDE) integration.

Stochastic Processes¶

`wiener_process`([paths, vshape, dtype, rng, …])	Wiener process (Brownian motion) with drift.
`lognorm_process`([paths, vshape, dtype, rng, …])	Lognormal process.
`ornstein_uhlenbeck_process`([paths, vshape, …])	Ornstein-Uhlenbeck process (mean-reverting Brownian motion).
`hull_white_process`([paths, vshape, dtype, …])	F-factors Hull-White process (sum of F correlated mean-reverting Brownian motions).
`hull_white_1factor_process`([paths, vshape, …])	1-factor Hull-White process (F=1 Hull-White process with F-index collapsed to a scalar).
`cox_ingersoll_ross_process`([paths, vshape, …])	Cox-Ingersoll-Ross mean reverting process.
`full_heston_process`([paths, vshape, dtype, …])	Heston stochastic volatility process (returns both process and volatility).
`heston_process`([paths, vshape, dtype, rng, …])	Heston stochastic volatility process (stores and returns process only).
`jumpdiff_process`([paths, vshape, dtype, …])	Jump-diffusion process (lognormal process with compound Poisson logarithmic jumps).
`merton_jumpdiff_process`([paths, vshape, …])	Merton jump-diffusion process (jump-diffusion process with normal jump size distribution).
`kou_jumpdiff_process`([paths, vshape, dtype, …])	Double exponential (Kou) jump-diffusion process (jump-diffusion process with double exponential jump size distribution).
`wiener_SDE`(*[, paths, vshape, dtype, rng, …])	SDE for a Wiener process (Brownian motion) with drift.
`lognorm_SDE`(*[, paths, vshape, dtype, rng, …])	SDE for a lognormal process with drift.
`ornstein_uhlenbeck_SDE`(*[, paths, vshape, …])	SDE for an Ornstein-Uhlenbeck process.
`hull_white_SDE`(*[, paths, vshape, dtype, …])	SDE for an F-factors Hull White process.
`cox_ingersoll_ross_SDE`(*[, paths, vshape, …])	SDE for a Cox-Ingersoll-Ross mean reverting process.
`full_heston_SDE`(*[, paths, vshape, dtype, …])	SDE for a Heston stochastic volatility process.
`heston_SDE`(*[, paths, vshape, dtype, rng, …])	SDE for a Heston stochastic volatility process.
`jumpdiff_SDE`(*[, paths, vshape, dtype, rng, …])	SDE for a jump-diffusion process (lognormal process with compound Poisson logarithmic jumps).
`merton_jumpdiff_SDE`(*[, paths, vshape, …])	SDE for a Merton jump-diffusion process.
`kou_jumpdiff_SDE`(*[, paths, vshape, dtype, …])	SDE for a double exponential (Kou) jump-diffusion process.

Analytical Results¶

`wiener_mean`(t, *[, x0, mu, sigma])	Mean of values at time t of a Wiener process (as per the wiener_process class) with time-independent parameters.
`wiener_var`(t, *[, x0, mu, sigma])	Variance of values at time t of a Wiener process (as per the wiener_process class) with time-independent parameters.
`wiener_std`(t, *[, x0, mu, sigma])	Standard deviation of values at time t of a Wiener process (as per the wiener_process class) with time-independent parameters.
`wiener_pdf`(t, x, *[, x0, mu, sigma])	Probability distribution function of values at time t of a Wiener process (as per the wiener_process class) with time-independent parameters, evaluated at x.
`wiener_cdf`(t, x, *[, x0, mu, sigma])	Cumulative probability distribution function of values at time t of a Wiener process (as per the wiener_process class) with time-independent parameters, evaluated at x.
`wiener_chf`(t, u, *[, x0, mu, sigma])	Characteristic function of the probability distribution of values at time t of a Wiener process (as per the wiener_process class) with time-independent parameters, evaluated at u.
`lognorm_mean`(t, *[, x0, mu, sigma])	Mean of values at time t of a lognormal process (as per the lognorm_process class) with time-independent parameters.
`lognorm_var`(t, *[, x0, mu, sigma])	Variance of values at time t of a lognormal process (as per the lognorm_process class) with time-independent parameters.
`lognorm_std`(t, *[, x0, mu, sigma])	Standard deviation of values at time t of a lognormal process (as per the lognorm_process class) with time-independent parameters.
`lognorm_pdf`(t, x, *[, x0, mu, sigma])	Probability distribution function of values at time t of a lognormal process (as per the lognorm_process class) with time-independent parameters, evaluated at x.
`lognorm_cdf`(t, x, *[, x0, mu, sigma])	Cumulative probability distribution function of values at time t of a lognormal process (as per the lognorm_process class) with time-independent parameters, evaluated at x.
`lognorm_log_chf`(t, u, *[, x0, mu, sigma])	Characteristic function of the probability distribution of values at time t of the logarithm of a lognormal process (as per the lognorm_process class) with time-independent parameters, evaluated at u.
`oruh_mean`(t, *[, x0, theta, k, sigma])	Mean of values at time t of an Ornstein-Uhlenbeck process (as per the ornstein_uhlenbeck_process class) with time-independent parameters.
`oruh_var`(t, *[, x0, theta, k, sigma])	Variance of values at time t of an Ornstein-Uhlenbeck process (as per the ornstein_uhlenbeck_process class) with time-independent parameters.
`oruh_std`(t, *[, x0, theta, k, sigma])	Standard deviation of values at time t of an Ornstein-Uhlenbeck process (as per the ornstein_uhlenbeck_process class) with time-independent parameters.
`oruh_pdf`(t, x, *[, x0, theta, k, sigma])	Probability distribution function of values at time t of an Ornstein-Uhlenbeck process (as per the ornstein_uhlenbeck_process class) with time-independent parameters, evaluated at x.
`oruh_cdf`(t, x, *[, x0, theta, k, sigma])	Cumulative probability distribution function of values at time t of an Ornstein-Uhlenbeck process (as per the ornstein_uhlenbeck_process class) with time-independent parameters, evaluated at x.
`hw2f_mean`(t, *[, x0, theta, k, sigma, rho])	Mean of values at time t of a Hull-White 2-factors process (as per the hull_white_process class) with time-independent parameters.
`hw2f_var`(t, *[, x0, theta, k, sigma, rho])	Variance of values at time t of a Hull-White 2-factors process (as per the hull_white_process class) with time-independent parameters.
`hw2f_std`(t, *[, x0, theta, k, sigma, rho])	Standard deviation of values at time t of a Hull-White 2-factors process (as per the hull_white_process class) with time-independent parameters.
`hw2f_pdf`(t, x, *[, x0, theta, k, sigma, rho])	Probability distribution function of values at time t of a Hull-White 2-factors process (as per the hull_white_process class) with time-independent parameters, evaluated at x.
`hw2f_cdf`(**args)	Cumulative probability distribution function of values at time t of a Hull-White 2-factors process (as per the hull_white_process class) with time-independent parameters, evaluated at x.
`cir_mean`(t, *[, x0, theta, k, xi])	Mean of values at time t of a Cox-Ingersoll-Ross process (as per the cox_ingersoll_ross_process class) with time-independent parameters.
`cir_var`(t, *[, x0, theta, k, xi])	Variance of values at time t of a Cox-Ingersoll-Ross process (as per the cox_ingersoll_ross_process class) with time-independent parameters.
`cir_std`(t, *[, x0, theta, k, xi])	Standard deviation of values at time t of a Cox-Ingersoll-Ross process (as per the cox_ingersoll_ross_process class) with time-independent parameters.
`cir_pdf`(t, x, *[, x0, theta, k, xi])	Probability distribution function of values at time t of a Cox-Ingersoll-Ross process (as per the cox_ingersoll_ross_process class) with time-independent parameters, evaluated at x.
`heston_log_mean`(t, *[, x0, mu, sigma, y0, …])	Mean of the logarithm of values at time t of a Heston process (as per the full_heston_process class) with time-independent parameters.
`heston_log_var`(**args)	Variance of the logarithm of values at time t of a Heston process (as per the full_heston_process class) with time-independent parameters.
`heston_log_std`(t, *[, x0, mu, sigma, y0, …])	Standard deviation of the logarithm of values at time t of a Heston process (as per the full_heston_process class) with time-independent parameters.
`heston_log_pdf`(t, logx, *[, x0, mu, sigma, …])	Probability distribution function of values at time t of the logarithm of a Heston process, (as per the full_heston_process class) with time-independent parameters, evaluated at logx.
`heston_log_chf`(t, u, *[, x0, mu, sigma, y0, …])	Characteristic function of the probability distribution of values at time t of the logarithm of a Heston process (as per the full_heston_process class) , with time-independent parameters, evaluated at u.
`mjd_log_pdf`(t, logx, *[, x0, mu, sigma, …])	Probability distribution function of values at time t of the logarithm of a Merton jump-diffusion process (as per the merton_jumpdiff_process class), with time-independent parameters, evaluated at logx.
`mjd_log_chf`(t, u, *[, x0, mu, sigma, lam, a, b])	Characteristic function of the probability distribution of values at time t of the logarithm of a Merton jump-diffusion process (as per the merton_jumpdiff_process class), with time-independent parameters, evaluated at u.
`kou_mean`(t, *[, x0, mu, sigma, lam, a, b, pa])	Mean of values at time t of a double exponential (Kou) jump-diffusion process (as per the kou_jumpdiff_process class) with time-independent parameters.
`kou_log_pdf`(t, logx, *[, x0, mu, sigma, …])	Probability distribution function of values at time t of the logarithm of a double-exponential (Kou) jump-diffusion process (as per the kou_jumpdiff_process class), with time-independent parameters, evaluated at logx.
`kou_log_chf`(t, u, *[, x0, mu, sigma, lam, …])	Characteristic function of the probability distribution of values at time t of the logarithm of a Kou jump-diffusion process, (as per the kou_jumpdiff_process class) with time-independent parameters, evaluated at u.
`bsd1d2`(k, t, *[, x0, r, q, sigma])	Black-Scholes d1 and d2 coefficients.
`bscall`(k, t, *[, x0, r, q, sigma])	Black-Scholes call option value.
`bscall_delta`(k, t, *[, x0, r, q, sigma])	Black-Scholes call option delta.
`bsput`(k, t, *[, x0, r, q, sigma])	Black-Scholes put option value.
`bsput_delta`(k, t, *[, x0, r, q, sigma])	Black-Scholes put option delta.

Shortcuts¶

Stochasticity sources and preset processes may be addressed using the following shortcuts:

Full name	Shortcut
`wiener_source`	`dw`
`poisson_source`	`dn`
`cpoisson_source`	`dj`
`odd_wiener_source`	`odd_dw`
`even_poisson_source`	`even_dn`
`even_cpoisson_source`	`even_dj`
`true_wiener_source`	`true_dw`
`true_poisson_source`	`true_dn`
`true_cpoisson_source`	`true_dj`
`wiener_process`	`wiener`
`lognorm_process`	`lognorm`
`ornstein_uhlenbeck_process`	`oruh`
`hull_white_process`	`hwff`
`hull_white_1factor_process`	`hw1f`
`cox_ingersoll_ross_process`	`cir`
`full_heston_process`	`heston_xy`
`heston_process`	`heston`
`jumpdiff_process`	`jumpdiff`
`merton_jumpdiff_process`	`mjd`
`kou_jumpdiff_process`	`kou`

Shortcuts have been wrapped as “kfuncs”, objects with managed keyword arguments (see kfunc decorator documentation below).

Analytical results are named according to the shortcut of the corresponding process (e.g. lognorm_mean, lognorm_cdf etc. from the lognorm shortcut) and are wrapped as kfuncs as well.

`kfunc`([f, nvar])	Decorator to wrap classes or functions as objects with managed keyword arguments.
`iskfunc`(cls_or_object)	Tests if the given class or instance has been wrapped as a kfunc.

Testing¶

Tests have been set up within the numpy.testing framework, and require the pytest package to be installed. To launch tests, invoke sdepy.test() or sdepy.test('full').

The testing subpackage sdepy.tests was written in pursuit of the following goals:

Maximize case coverage, by exposing the package functions and methods to a plurality of different input shapes, values, data types etc., and of different combinations thereof, as may be encountered in practice.
Provide a quantitative validation of the algorithms, functions and processes covered in sdepy.
Keep dependencies of the test code on the adopted testing framework to a bare minimum.

Most often, a number of testing cases is declared as a list or lists of classes and inputs, a general testing procedure is set up, and the latter is iteratively applied to the former. Unfortunately, all this resulted in a thinly documented (if at all), hard to read, and hard to maintain testing code base - sorry about that.

The quantitative validation of the package, via tests marked as 'slow' and 'quant', is done in two steps:

To validate a sdepy release, tests are run with 100_000 or more paths. Numerical integration results for the mean, standard deviation, probability distribution, and/or characteristic function are compared against their exact values computed analytically from the process parameters. Comparisons are then plotted and visually inspected, and the occasional larger than usual deviation is manually checked to be statistically acceptable, i.e. only so few standard deviations off the mark.
Each time sdepy.test('full') is invoked, to keep testing times manageable and the testing procedure non invasive, tests are run with 100 paths, using NumPy legacy random generation with a fixed seed, the realized errors are then compared and checked against the expected errors, as distributed with the package and stored in the ./tests/cfr directory. Note that such default tests rely on the stable reproducibility of expected errors, across platforms and versions of Python, NumPy and SciPy.

In order to run tests using current NumPy random generation, with a set number of paths (200_000 in the example below), and inspect graphs and realized errors as saved in the current directory, use the following statement (non backward compatible changes to the testing interface may occur in the future, without prior warning):

sdepy.test('full', rng=numpy.random.default_rng(),
           paths=200_000, plot=True, outdir='.')

For further details, see sdepy.test documentation:

sdepy.test Invoke the sdepy testing suite (requires pytest>=3.8.1 to be installed).