pyNSDE - solving nonlinear stochastic differential equation

Here we have implemented numerical solution of nonlinear SDE,

$$ d x = \left(\eta - \frac{\lambda}{2} \right) x^{2 \eta - 1} d t + x^\eta d W . $$

By "solution" we mean obtaining a sample trajectory of the process. Notably, this SDE is written in Ito sense.

Instead of solving the above NSDE directly, we solve SDE of the Bessel process,

$$ d y = \frac{\lambda - \eta}{2 \left( \eta - 1 \right)} \cdot \frac{d t}{y} + d W , $$

using Euler-Maruyama method with variable time step. Notably, this SDE is also written in Ito sense.

Sample trajectories of the Bessel process, $y$, are then transformed into sample trajectories of the NSDE above, $x$. Exact form of transformation depends on the parameter $\eta$ value. For $\eta > 1 $,

$$ x = \left[ \left( \eta - 1 \right) y \right]^\frac{1}{1-\eta} $$

is used. While, for $\eta < 1$,

$$ x = \left[ \left( 1 - \eta \right) y \right]^\frac{1}{1-\eta} $$

is employed instead. $\eta = 1$ case should be avoided, as Bessel process diverges in this case.

Usage

Under the hood this Python module uses C program through ctypes built-in module. Thus you'll need to compile C code first (makefile is provided for your convenience). Note that C code depends on GNU Scientific Library.

After compiling the C code, you can use this library to generate time series that exhibit pink or $1/f$ noise. Depending on the model and simulation parameters, $1/f$ noise can be observed in an arbitrarily broad range of frequencies. Below follows an example with simulation results of the calculation with mostly default parameter values.

import numpy as np
import matplotlib.pyplot as plt

from pyNSDE import generate_series

from stats.pdf import make_log_pdf
from stats.psd import make_seg_log_psd

# simulation
series = generate_series(1048576, 1e-3, seed=123)

# calculating PDF / PSD
pdf = make_log_pdf(series)
psd = make_seg_log_psd(series, fs=1e3)

# creating simple visualization
plt.figure(figsize=(12,3))
plt.subplot(131)
plt.xlabel('t')
plt.ylabel('x(t)')
plt.plot(series[::256], 'r-')
plt.subplot(132)
plt.loglog()
plt.xlabel('x')
plt.ylabel('p(x)')
plt.plot(pdf[:, 0], pdf[:, 1], 'r-')
plt.plot(pdf[:, 0], 2*(pdf[:, 0]**-3), 'k--')
plt.subplot(133)
plt.loglog()
plt.xlabel('f')
plt.ylabel('S(f)')
plt.plot(psd[:, 0], psd[:, 1], 'r-')
plt.plot(psd[20:, 0], 1.5*(psd[20:, 0]**-1), 'k--')
plt.tight_layout()
plt.show()

In this code snippet stats library was cloned from https://github.com/akononovicius/python-stats.

Related repositories

Earlier implementation of a less flexible program solving the same nonlinear stochastic differential equation is available at https://github.com/JuliusRuseckas/numerical-sde-variable-step.

https://github.com/akononovicius/python-stats library might be useful when analyzing simulated time series.

Related research

Couple of scientific review papers specific to the SDE being solved:

• B. Kaulakys and J. Ruseckas, Stochastic nonlinear differential equation generating 1/f noise, Phys. Rev. E 70, 020101 (2004). doi: 10.1103/PhysRevE.70.020101. arXiv:cond-mat/0408507 [cond-mat.stat-mech].
• B. Kaulakys, J. Ruseckas, V. Gontis and M. Alaburda, Nonlinear stochastic models of 1/f noise and power-law distributions, Physica A 365, 217-221 (2006). doi: 10.1016/j.physa.2006.01.017. arXiv:cond-mat/0509626 [cond-mat.stat-mech].

Recent review with variety of applications related to the SDE being solved:

• R. Kazakevičius, A. Kononovicius, B. Kaulakys, V. Gontis. Understanding the nature of the long-range memory phenomenon in socio-economic systems. Entropy 23: 1125 (2021). doi: 10.3390/e23091125. arXiv:2108.02506 [physics.soc-ph].