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.,   1.,   6.,   9.,   9.,  22.,  35.,  87., 112., 186., 276.,
       408., 506., 692., 727., 803., 852., 764., 651., 591., 477., 364.,
       258., 147., 105.,  55.,  19.,  16.,   8.,   5.]), array([-0.43577188, -0.34742472, -0.25907756, -0.17073039, -0.08238323,
        0.00596393,  0.09431109,  0.18265826,  0.27100542,  0.35935258,
        0.44769975,  0.53604691,  0.62439407,  0.71274123,  0.8010884 ,
        0.88943556,  0.97778272,  1.06612989,  1.15447705,  1.24282421,
        1.33117137,  1.41951854,  1.5078657 ,  1.59621286,  1.68456002,
        1.77290719,  1.86125435,  1.94960151,  2.03794868,  2.12629584,
        2.214643  ]), <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.0011
>>> print(f"simulation variance: {np.var(XF) :0.5}")
simulation variance: 0.12365
>>> 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