Stochastic simulations

Stochastic simulations are supported via the Euler-Maruyama method. Random normal variables \(W_t\sim N(0, dt)\) can be included in the vector field function, and each simulation instance will have an independent stream of random variables.

As an example, consider the Ornstein-Uhlenbeck process with \(\theta = 1\):

\[ dx_t = (\mu - x_t)dt + \sigma W_t \]

To implement this system, we heuristically write it in Langevin form:

\[ \frac{dx}{dt} = \mu - x + \sigma\eta(t)\]

where \(\eta(t)\) represents white noise. When \(\sigma=0\), the system is a first-order ODE with steady state x=\(\mu\). Increasing \(\sigma\) produces solutions with a steady state distribution \(x\sim N(\mu, \frac{\sigma^2}{2})\).

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import clode
>>> from typing import List
>>> def wiener_process(
...         t: float,
...         y: List[float],
...         p: List[float],
...         dy: List[float],
...         aux: List[float],
...         wiener: List[float],
... ) -> None:
...     x: float = y[0]
...     mu: float = p[0]
...     sigma: float = p[1]
...     weiner_variable: float = wiener[0]
...     dx: float = mu - x + sigma * weiner_variable
...     dy[0] = dx
... 
>>> variables = {"x": 0.0}
>>> parameters = {"mu": 1.0, "sigma": 0.5}
>>> t_span = (0.0, 1000.0)
>>> integrator = clode.Simulator(
...     rhs_equation=wiener_process,
...     variables=variables,
...     parameters=parameters,
...     num_noise=1,
...     stepper=clode.Stepper.stochastic_euler,
...     single_precision=True,
...     t_span=t_span,
...     dt=0.001,
... )
>>> nPts = 8192
>>> integrator.set_repeat_ensemble(nPts)
>>> integrator.transient()
>>> XF = integrator.get_final_state()
>>> plt.hist(XF, 30)
(array([  1.,   2.,   5.,  18.,  35.,  62., 117., 180., 269., 422., 576.,
       736., 762., 813., 876., 810., 717., 569., 432., 326., 193., 126.,
        69.,  31.,  25.,   6.,   6.,   5.,   1.,   2.]), array([-0.3228741 , -0.22877825, -0.1346824 , -0.04058655,  0.0535093 ,
        0.14760516,  0.24170101,  0.33579686,  0.42989271,  0.52398856,
        0.61808441,  0.71218026,  0.80627611,  0.90037196,  0.99446781,
        1.08856367,  1.18265952,  1.27675537,  1.37085122,  1.46494707,
        1.55904292,  1.65313877,  1.74723462,  1.84133047,  1.93542632,
        2.02952218,  2.12361803,  2.21771388,  2.31180973,  2.40590558,
        2.50000143]), <BarContainer object of 30 artists>)
>>> plt.xlabel("x")
Text(0.5, 0, 'x')
>>> plt.show()
>>> print(f"mean xf: {np.mean(XF) :0.5}")
mean xf: 1.0036
>>> print(f"simulation variance: {np.var(XF) :0.5}")
simulation variance: 0.12492
>>> print(f"expected variance: {parameters['sigma'] ** 2 / 2 :0.5}")
expected variance: 0.125

Output:

mean xf: 1.0041
simulation variance: 0.12629
expected variance: 0.125