.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/statphys/stochastic_processes/plot_ornstein_uhlenbeck_process.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_statphys_stochastic_processes_plot_ornstein_uhlenbeck_process.py: The Ornstein-Uhlenbeck process: Gaussian noise with a memory ================================================================ In 1930 George Uhlenbeck and Leonard Ornstein solved Langevin's equation for the velocity of a free Brownian particle exactly. The velocity is not white noise. It is a process pulled back toward its mean at a rate :math:`\theta` while being kicked with strength :math:`\sigma`, .. math:: dX = \theta(\mu - X)\,dt + \sigma\,dW. From a fixed start its mean relaxes as :math:`e^{-\theta t}` and its variance fills in as :math:`\frac{\sigma^2}{2\theta}(1 - e^{-2\theta t})`. Once stationary, it forgets its past exponentially, :math:`\langle \delta X(0)\,\delta X(\tau)\rangle = \frac{\sigma^2}{2\theta}e^{-\theta|\tau|}`. Doob (1942) later showed it is the only process that is at once stationary, Gaussian, and Markov. :class:`~physicskit.statphys.OrnsteinUhlenbeck` samples it with the exact Gaussian transition, or by Euler-Maruyama integration of the SDE. This example checks every one of these formulas, and shows the small, predictable bias that time stepping adds to the stationary variance. .. GENERATED FROM PYTHON SOURCE LINES 26-34 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.statphys import OrnsteinUhlenbeck theta, mu, sigma = 1.0, 0.0, 1.0 .. GENERATED FROM PYTHON SOURCE LINES 35-37 Relaxation from a fixed start ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 37-42 .. code-block:: Python ou = OrnsteinUhlenbeck(theta=theta, mu=mu, sigma=sigma, n_paths=20000, seed=1930) t, X = ou.run(t_max=5.0, dt=0.01, n_frames=250, x0=3.0) print(f"stationary variance sigma^2 / 2 theta = {ou.stationary_variance:.4f}; sampled at t = 5: {X[-1].var():.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none stationary variance sigma^2 / 2 theta = 0.5000; sampled at t = 5: 0.4924 .. GENERATED FROM PYTHON SOURCE LINES 43-48 Stationary memory ----------------- Starting from the stationary distribution, the autocovariance at lag :math:`\tau` decays as :math:`e^{-\theta\tau}`: the process remembers its value for a time :math:`1/\theta`. .. GENERATED FROM PYTHON SOURCE LINES 48-53 .. code-block:: Python t_s, X_s = ou.run(t_max=5.0, dt=0.01, n_frames=100) lags = np.arange(0, 60) acov = np.array([np.mean((X_s[0] - mu) * (X_s[lag] - mu)) for lag in lags]) .. GENERATED FROM PYTHON SOURCE LINES 54-59 Time-step bias of Euler-Maruyama -------------------------------- The exact update has no time-step error. Euler-Maruyama overshoots the stationary variance by the factor :math:`1/(1 - \theta\Delta t/2)`, which is negligible for :math:`\theta\Delta t \ll 1`. .. GENERATED FROM PYTHON SOURCE LINES 59-97 .. code-block:: Python dts = np.array([0.02, 0.1, 0.25, 0.5, 1.0]) em_var = [] for dt in dts: em = OrnsteinUhlenbeck(theta=theta, mu=mu, sigma=sigma, n_paths=40000, seed=int(1000 * dt)) _, Xe = em.run(t_max=8.0, dt=dt, n_frames=4, x0=0.0, method="euler_maruyama") em_var.append(Xe[-1].var()) print(f"dt = {dt:4.2f}: Euler-Maruyama variance {em_var[-1]:.4f}, predicted {sigma**2 / (theta * (2 - theta * dt)):.4f}") fig, axes = plt.subplots(1, 3, figsize=(15, 4.3)) for i in range(6): axes[0].plot(t, X[:, i], lw=0.6) axes[0].plot(t, X.mean(axis=1), "k", lw=2, label="ensemble mean") axes[0].plot(t, ou.mean(t, 3.0), "w--", lw=1) sd = np.sqrt(ou.variance(t)) axes[0].fill_between(t, ou.mean(t, 3.0) - sd, ou.mean(t, 3.0) + sd, color="gray", alpha=0.3, label=r"exact mean $\pm$ sd") axes[0].set_xlabel("t") axes[0].set_ylabel("X") axes[0].set_title("relaxing from X = 3") axes[0].legend(fontsize=8) axes[1].plot(t_s[lags], acov, "o", ms=3, label="sampled") axes[1].plot(t_s[lags], ou.autocovariance(t_s[lags]), "k--", label=r"$\frac{\sigma^2}{2\theta}e^{-\theta\tau}$") axes[1].set_xlabel(r"lag $\tau$") axes[1].set_title("stationary autocovariance") axes[1].legend() axes[2].plot(dts, em_var, "o", label="Euler-Maruyama") fine = np.linspace(0, 1.05, 100) axes[2].plot(fine, sigma**2 / (theta * (2 - theta * fine)), "k--", lw=1, label=r"$\sigma^2/\theta(2-\theta\Delta t)$") axes[2].axhline(ou.stationary_variance, color="gray", ls=":", label="exact") axes[2].set_xlabel(r"$\theta\,\Delta t$") axes[2].set_ylabel("stationary variance") axes[2].set_title("time-step bias") axes[2].legend() fig.tight_layout() plt.show() .. image-sg:: /api/gallery/statphys/stochastic_processes/images/sphx_glr_plot_ornstein_uhlenbeck_process_001.png :alt: relaxing from X = 3, stationary autocovariance, time-step bias :srcset: /api/gallery/statphys/stochastic_processes/images/sphx_glr_plot_ornstein_uhlenbeck_process_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none dt = 0.02: Euler-Maruyama variance 0.5093, predicted 0.5051 dt = 0.10: Euler-Maruyama variance 0.5225, predicted 0.5263 dt = 0.25: Euler-Maruyama variance 0.5714, predicted 0.5714 dt = 0.50: Euler-Maruyama variance 0.6619, predicted 0.6667 dt = 1.00: Euler-Maruyama variance 0.9916, predicted 1.0000 .. GENERATED FROM PYTHON SOURCE LINES 98-102 Check ----- Stationary variance sigma^2 / 2 theta, autocovariance e^(-theta tau), and the Euler-Maruyama bias sigma^2 / theta (2 - theta dt). .. GENERATED FROM PYTHON SOURCE LINES 102-105 .. code-block:: Python assert abs(X[-1].var() / ou.stationary_variance - 1) < 0.03 assert np.max(np.abs(acov - ou.autocovariance(t_s[lags]))) < 0.02 np.testing.assert_allclose(em_var, sigma**2 / (theta * (2 - theta * dts)), rtol=0.03) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.493 seconds) .. _sphx_glr_download_api_gallery_statphys_stochastic_processes_plot_ornstein_uhlenbeck_process.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_ornstein_uhlenbeck_process.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_ornstein_uhlenbeck_process.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_ornstein_uhlenbeck_process.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_