crflynn / stochastic

Generate realizations of stochastic processes in python.

stochastic

A python package for generating realizations of stochastic processes.

Installation

The stochastic package is available on pypi and can be installed using pip

pip install stochastic

Dependencies

Stochastic uses numpy for many calculations and scipy for sampling specific random variables.

Processes

This package offers a number of common discrete-time, continuous-time, and noise process objects for generating realizations of stochastic processes as numpy arrays.

The diffusion processes are approximated using the Euler–Maruyama method.

Here are the currently supported processes and their class references within the package.

• stochastic.processes

  • continuous

    • BesselProcess
    • BrownianBridge
    • BrownianExcursion
    • BrownianMeander
    • BrownianMotion
    • CauchyProcess
    • FractionalBrownianMotion
    • GammaProcess
    • GeometricBrownianMotion
    • InverseGaussianProcess
    • MixedPoissonProcess
    • MultifractionalBrownianMotion
    • PoissonProcess
    • SquaredBesselProcess
    • VarianceGammaProcess
    • WienerProcess

  • diffusion

    • DiffusionProcess (generalized)
    • ConstantElasticityVarianceProcess
    • CoxIngersollRossProcess
    • ExtendedVasicekProcess
    • OrnsteinUhlenbeckProcess
    • VasicekProcess

  • discrete

    • BernoulliProcess
    • ChineseRestaurantProcess
    • DirichletProcess
    • MarkovChain
    • MoranProcess
    • RandomWalk

  • noise

    • BlueNoise
    • BrownianNoise
    • ColoredNoise
    • PinkNoise
    • RedNoise
    • VioletNoise
    • WhiteNoise
    • FractionalGaussianNoise
    • GaussianNoise

Usage patterns

Sampling

To use stochastic, import the process you want and instantiate with the required parameters. Every process class has a sample method for generating realizations. The sample methods accept a parameter n for the quantity of steps in the realization, but others (Poisson, for instance) may take additional parameters. Parameters can be accessed as attributes of the instance.

from stochastic.processes.discrete import BernoulliProcess


bp = BernoulliProcess(p=0.6)
s = bp.sample(16)
success_probability = bp.p

Continuous processes provide a default parameter, t, which indicates the maximum time of the process realizations. The default value is 1. The sample method will generate n equally spaced increments on the interval [0, t].

Sampling at specific times

Some continuous processes also provide a sample_at() method, in which a sequence of time values can be passed at which the object will generate a realization. This method ignores the parameter, t, specified on instantiation.

from stochastic.processes.continuous import BrownianMotion


bm = BrownianMotion(drift=1, scale=1, t=1)
times = [0, 3, 10, 11, 11.2, 20]
s = bm.sample_at(times)

Sample times

Continuous processes also provide a method times() which generates the time values (using numpy.linspace) corresponding to a realization of n steps. This is particularly useful for plotting your samples.

import matplotlib.pyplot as plt
from stochastic.processes.continuous import FractionalBrownianMotion


fbm = FractionalBrownianMotion(hurst=0.7, t=1)
s = fbm.sample(32)
times = fbm.times(32)

plt.plot(times, s)
plt.show()

Specifying an algorithm

Some processes provide an optional parameter algorithm, in which one can specify which algorithm to use to generate the realization using the sample() or sample_at() methods. See the documentation for process-specific implementations.

from stochastic.processes.noise import FractionalGaussianNoise


fgn = FractionalGaussianNoise(hurst=0.6, t=1)
s = fgn.sample(32, algorithm='hosking')

stochastic's Issues

Geometric Brownian Motion process seemingly can't be seeded / results can't be reproduced

Hi there,

I created an MVP of a piece of code that has been causing issues, I don't think it's intended behaviour unless I'm misunderstanding some difference between GBM and BM processes. Whereas the BrownianMotion process will always return the same result when seeded, the GeometricBrownianMotion process does not. I have tried a couple of variations of the below code to test that it's not some bug in the order of code execution or seeding of the RNG. I did dig into the stochastic code a bit, but didn't find any obvious errors - maybe a maintainer of this codebase has a better idea?

MVP:

from stochastic import processes
import numpy as np


rng = np.random.default_rng(1)

gbm_process = processes.continuous.GeometricBrownianMotion(rng=rng)
gbm_sample = gbm_process.sample(1)

bm_process = processes.continuous.BrownianMotion(rng=rng)
bm_sample = bm_process.sample(1)

print(f"Geometric Brownian Motion sample: {gbm_sample}")
print(f"Brownian Motion sample: {bm_sample}")

Output:

Geometric Brownian Motion sample: [1.         0.62310484]
Brownian Motion sample: [0.         0.34558419]

Thanks!

Support RNG seeding

It would be nice if the RNG used could be seeded (e.g. for reproducibility). Even better would be if users could just pass a RNG instance e.g. to the class constructor.

I can probably contribute such a feature (e.g. adding a rng kwarg to all classes), but before doing so, I would like to know 1) whether you agree at least with the basic idea and 2) whether you would be willing to drop support for numpy<1.17, as numpy1.17 introduced a new RNG API (https://numpy.org/doc/1.18/reference/random/index.html) which also has better performance (in particular using Ziggurat instead of Box-Muller for normal variates) -- it should not be too difficult to have a branch in the code to support the older API too, but I just want to know whether I can skip that work ;-)

The argument volexp behaves different to other arguments in the Difussion process.

Shouldn't the argument to volexp be t?

stochastic/stochastic/processes/diffusion/diffusion.py

Line 122 in 90ffba8

+ self._vol(t) * initial ** self._volexp(initial) * gns[k]

New Continuous Process: Non-Homogeneous Poisson

This continuous stochastic process generates points distributed according to a density function or density matrix. Can be used to generate poisson processes whose rate function varies with time, space, or any other data space.

The user would need to supply a deterministic function, or matrix representing the density in the data space, as well as the boundaries of this data space. These can be n-dimensional. Then, points are generated using the thinning/acceptance-rejection algorithm. This necessitaets the generation of a maximum lambda value in the data space using the private method _gen_lmax, generating uniformly distributed points in the space with rate lmax (in unthinned), then rejecting these points with probability proportional to the density at each generated point (getting thinned).

I found the NonHomogeneousPoissonProcess too dissimilar to inherit from the PoissonProcess. I hope this pull request is closer to being acceptable than my previous one.

New Continuous Process: Multidimensional Spatial Point Process.

This process would be able to generate any number or events in an arbitrary number of dimensions, using either a callable function or density array. This class would diverge slightly from other processes, maybe changing the length parameter to bounds as such: ((0,1),(0,2),...). Then, the generated n would only be generated in these bounds. For now, I am thinking about generating points using the thinning algorithm.

Would such a process fit in this package?

Wrong implementation of _fgn_autocovariance

Apologies if this is not the right to communicate with the developers. By looking at the code for generating fractional gaussian noise using the Davies and Harte method, I see an issue with the calculation of the circulant matrix first row. As shown in the appendix of the 1987 paper, the first row contains the autocovariances c0, c1,...cn, cn-1, ...c1. However, the current procedure returns a sequence of cs which just gives different values for each i in [0, n]. What I am saying is also corroborated by the description in Wood and Chan (1994 and 1997) whose method in 1 dimension is identical to Davies and Harte.

If you do not agree this is a problem I would love to hear the explanation for my own understanding. Otherwise, an implementation such as

H2 = 2 * hurst
k = np.hstack([np.arange(0, n // 2 + 1), np.arange(n // 2 - 1, 0, -1)])
return 0.5 * (abs(k - 1) ** H2 - 2 * k ** H2 + (k + 1) ** H2)

Does I believe exhibit the correct behavior.

Added to conda

Hi,

this is more of an FYI than an issue but I just wanted to let the maintainers know that I added stochastic to conda-forge so that it can be downloaded with conda.

The repo is here: https://github.com/conda-forge/stochastic-feedstock

Let me know if anyone of the maintainers here would like to be added as maintainers there or if I should make a PR here to add a badge or install instructions. Otherwise, feel free to just close the issue.

Thanks! 😃

BrownianMotion sample_at losing the drift and scale

Hello, I believe that _sample_brownian_motion_at is wrong as it calls it calls _sample_gaussian_noise_at and doesn't adjust for drift and scale. As a consequence, sample_at samples a Guassian Process. Please let me know if you would like me to work on it.

Interface for iterative sampling

It would be interesting to have an interface for dynamic iterative sampling...

Instead of doing something like:

process.sample(255)

for a more dynamic process (like a simulation of some kind), I would do something like:

while True:
    process.sample(1).
    plots.update().
    # etc

but I am guessing there are restrictions on how to do the sampling in an iterative way, depending on the kind of stochastic process we want to simulate.

I only know about Markov Chain specifying that the next sample depends only on the one before, so I am assuming other kinds of stochastic processes have various constraints about this... hence why I'm thinking about an extra interface for iterative sampling to support that use-case.

Thoughts ?

Add type annotations

Python has type annotations! Adding them would enable editors to identify mistakes and would allow static typechecking as a test (i.e. via mypy).

Sampling from `FractionalGaussianNoise` and `ColoredNoise` leads to arrays of different size

Hi,

sampling from FractionalGaussianNoise and ColoredNoiseas follows:

import stochastic.processes.noise as sn
import numpy as np

rng = np.random.default_rng(seed=42)

n = 10000

br_noise = sn.ColoredNoise(beta=2.0, rng=rng).sample(n)
fr_noise =  sn.FractionalGaussianNoise(hurst=0.5, rng=rng).sample(n)

Results in arrays of unequal length:

br_noise.shape
# Out: (10001,)

fr_noise.shape
# Out: (10000,)

Is this the intended behaviour?

I believe both classes internally check that n is a positive integer. However, ColoredNoise then does n = n + 1

stochastic/stochastic/processes/noise/colored_noise.py

Line 68 in 90ffba8

n = n + 1

In BernoulliProcess(p=P), bp.p is actually probability of "0" instead of "1"

Thank you for Stochastic!
I was preparing a class and realized that bp.p will give failure probability rather than success probability.
Very simple fix.
Thanks

Multi-Dimensional Process

What about multi-dimensional version of the processes ?

Add a more general diffusion process

Would it be possible to add a class with a simple general difussion? Something like $dX_t = b(t, X_t) dt + \sigma(t, X_t) dW_t$ where the functions $\sigma$ and $b$ passed as arguments to the class initialization are both functions of time and space.

I imagine in some circunstances (non-Lipschitz) the behaviour of the solution can't be ensured or easily controlled. But a simple warning in the docs should suffice and maybe some examples of how it fails would be academically interesting for the users.

Best.

Update build badge to gha

Hosking method for fgn regenerates the autocovariance matrix at wrong times

In FractionalGaussianNoise._hosking, there's the snippet

            if self._n != n or self._cov is None:
                self._cov = np.array([self._autocovariance(i)
                                      for i in range(n)])

which suggests that you were intending to save the calculation of self._cov for later calls. However, _hosking() never sets self._n, which remains None, so in fact that line will be run evertime fgn.sample(n, algorithm="hosking") is called, even with the same n (this can be checked by adding a print there). That specific point is just a missed optimization; however, that also means that an (admittedly slightly pathological) sequence of calls will cause a crash:

from stochastic.noise import FractionalGaussianNoise
fgn = FractionalGaussianNoise(.75)
fgn.sample(1000, algorithm="daviesharte")  # sets self._n to 1000
fgn.sample(100, algorithm="hosking")  # caches a size-100 covariance vector
fgn.sample(1000, algorithm="hosking")  # crashes: "index 100 is out of bounds for axis 0 with size 100"

Variance gamma process sample depends on the step number

Hello, I believe that _sample_variance_gamma_process generates realizations that depends on the number of the steps at fixed maximum time.
If I change the n parameter from 100 to 1000 I obtain that the excursion of the trajectories change from [-0.002 : 0.002] to [-0.0006 to 0.0006]. This gives a 2nd moment depending on the time step of the realization. It sounds bad.
The other parameters are:
VarianceGammaProcess(drift=0, variance=10., scale=0.002, t=1, rng=None)
Would it be possible to fix this?
Thanks

Generating many samples

Hi all, and thank you for the great package!

When one is interested in generating samples of stochastic processes for Monte Carlo methods and Machine Learning tasks it is essential to create many samples of the underlying process. The way the package is designed, this is very inefficient and doesn't exploit the power of numpy. (Using the package, you have to essentially iterate in a python loop for something like 1M+ times).
These are the tasks I am working on the most (they are particularly important in mathematical finance applications) and I was thinking about publishing some of the code that I have collected so far. But before doing so, I was wondering if there is interest in updating this project to be more efficient for generating many samples. In that case a would be happy to help (e.g. I have a Cython implementation of Hosking's method).

New Continuous Process: Mixed Poisson

This continuous stochastic process generates points distributed according to a rate generated from an arbitrary random distribution. I refer to the random distribution used to generate this rate is known as the information distribution.

The user would need to supply a random distribution in its parameters upon first call (numpy functions such as numpy.random.uniform). Then, a rate would be generated upon initialization. Upon each call of the sample method, a new random rate would be generated using the information distribution.

Since this is a Poisson Process, the suggested MixedPoissonProcess class is very similar to the PoissonProcess class. Additionnal methods are setters and getters for the information process and parameters. Upon each call of the setters, a new rate is generated. Another additional method is the genrate method, which generates a new random rate using the previously supplied information distribution and parameters.