File:Ornstein-Uhlenbeck-traces-sigma-theta.svg

#! /usr/bin/env python3
# -*- coding:utf8 -*-

import matplotlib.pyplot as plt
import numpy as np
from math import *

plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
np.random.seed(8)

sigma_theta = ((1., 1.), (4., 1.), (4., 16.))
# sigma: diffusion
# theta: mean reversion rate
mu = 0. # long-term mean
t = np.linspace(0, 5, 2001)
dt = t[1:] - t[:-1]

fig = plt.figure(figsize=(520 / 90.0, 340 / 90.0), dpi=72)

for st in sigma_theta:
    sigma, theta = st
    D = sigma**2 / 2 # diffusion constant
    sigmaX = sigma / sqrt(2 * theta) # standard deviation of X
    randnorm = np.random.normal(0, 1, len(t))
    X = np.empty_like(t)
    X[0] = sigmaX * randnorm[0]

    for i in range(1, len(t)):
        X[i] = X[i-1] + (mu - X[i-1]) * dt[i-1] * theta + sigma * sqrt(dt[i-1]) * randnorm[i]
    
    txt_sigx = f'{sigmaX:.1f}'.replace('0.7', r'\frac{1}{2}\sqrt{2}').replace('2.8', r'2\sqrt{2}')
    plt.plot(t, X, lw=1, label=rf'$\sigma={sigma:.0f}$  $\theta=${theta:2.0f}  $\sigma_X=' + txt_sigx + '$')

plt.grid(True)
plt.xlim(t[0], t[-1])
plt.ylim(-7.2, 7.2)
plt.xlabel('t')
plt.ylabel('X')
plt.legend(loc='upper right', framealpha=1)
plt.tight_layout()
plt.savefig('Ornstein-Uhlenbeck-traces-sigma-theta.svg')