SDE Monte Carlo Distribution Estimates

This notebook is a practical SDE case study for pyna. It uses pyna’s SDE classes for the model boundary and single-path geometry, then uses vectorized NumPy ensembles for distribution estimates.

The notebook is intentionally pre-executed locally. The saved outputs are rendered by Sphinx/nbsphinx on GitHub Pages, so the documentation workflow does not need to spend CI CPU time on repeated Monte Carlo sampling.

%matplotlib inline
from matplotlib_inline.backend_inline import set_matplotlib_formats
set_matplotlib_formats('svg')

import numpy as np
import matplotlib.pyplot as plt

from pyna.dynamics import BrownianMotion, GeometricBrownianMotion, ItoSDE

plt.rcParams.update({
    "figure.figsize": (10, 4),
    "axes.grid": True,
    "axes.spines.top": False,
    "axes.spines.right": False,
})

rng = np.random.default_rng(20260701)
print("rng seed = 20260701")

rng seed = 20260701

Brownian motion: terminal distribution

A single Brownian realization is a Trajectory. A large ensemble estimates a probability distribution. For Brownian motion, the terminal distribution is analytic, so it is a useful regression check for Monte Carlo code.

bm = BrownianMotion(dim=1, diffusion=1.0, drift=[0.15])
single_path = bm.euler_maruyama([0.0], (0.0, 2.0), dt=0.002, rng=7)

n_paths = 60_000
n_steps = 600
T = 2.0
dt = T / n_steps

increments = np.sqrt(dt) * rng.normal(size=(n_paths, n_steps))
terminal = bm.drift_vector[0] * T + increments.sum(axis=1)

times = np.linspace(0.0, T, n_steps + 1)
fan_count = 40
fan_paths = np.empty((fan_count, n_steps + 1))
fan_paths[:, 0] = 0.0
fan_paths[:, 1:] = bm.drift_vector[0] * times[1:] + np.cumsum(increments[:fan_count], axis=1)

analytic_mean = bm.mean([0.0], T)[0]
analytic_var = bm.variance(T)[0]
empirical_mean = float(np.mean(terminal))
empirical_var = float(np.var(terminal, ddof=1))

x_grid = np.linspace(np.quantile(terminal, 0.001), np.quantile(terminal, 0.999), 400)
pdf = np.exp(-0.5 * (x_grid - analytic_mean) ** 2 / analytic_var) / np.sqrt(2 * np.pi * analytic_var)

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(times, fan_paths.T, color="#2f6f9f", alpha=0.08, linewidth=0.8)
axes[0].plot(single_path.t, single_path.y[:, 0], color="#b23a48", linewidth=1.6, label="one pyna Trajectory")
axes[0].set_title("Brownian sample paths")
axes[0].set_xlabel("t")
axes[0].set_ylabel("X(t)")
axes[0].legend(loc="upper left")

axes[1].hist(terminal, bins=90, density=True, color="#7aa95c", alpha=0.75, label="Monte Carlo")
axes[1].plot(x_grid, pdf, color="#1d3557", linewidth=2.0, label="analytic normal")
axes[1].set_title("Terminal distribution X(T)")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"paths = {n_paths:,}, steps = {n_steps}, T = {T}")
print(f"mean: empirical {empirical_mean:.5f}, analytic {analytic_mean:.5f}, abs error {abs(empirical_mean - analytic_mean):.5f}")
print(f"var : empirical {empirical_var:.5f}, analytic {analytic_var:.5f}, abs error {abs(empirical_var - analytic_var):.5f}")

paths = 60,000, steps = 600, T = 2.0
mean: empirical 0.29974, analytic 0.30000, abs error 0.00026
var : empirical 2.00725, analytic 2.00000, abs error 0.00725

Ornstein-Uhlenbeck: approach to a stationary law

The OU process is a compact teaching model for mean reversion. It also demonstrates a custom ItoSDE with user-defined drift and diffusion. The Monte Carlo ensemble checks the finite-time normal law and the stationary variance.

theta = 1.8
mean_level = 1.5
sigma = 0.35
x0 = -1.0
T = 3.0
n_steps = 1200
dt = T / n_steps
n_paths = 80_000

ou = ItoSDE(
    drift=lambda x, t: theta * (mean_level - x),
    diffusion=lambda x, t: np.array([sigma]),
    dim=1,
    brownian_dim=1,
    coordinate_names=("x",),
    label="Ornstein-Uhlenbeck",
)
ou_path = ou.euler_maruyama([x0], (0.0, T), dt=dt, rng=11)

x = np.full(n_paths, x0, dtype=float)
sqrt_dt_sigma = sigma * np.sqrt(dt)
for _ in range(n_steps):
    x += theta * (mean_level - x) * dt + sqrt_dt_sigma * rng.normal(size=n_paths)

analytic_mean = mean_level + (x0 - mean_level) * np.exp(-theta * T)
analytic_var = sigma**2 / (2 * theta) * (1 - np.exp(-2 * theta * T))
stationary_var = sigma**2 / (2 * theta)
empirical_mean = float(np.mean(x))
empirical_var = float(np.var(x, ddof=1))

grid = np.linspace(np.quantile(x, 0.001), np.quantile(x, 0.999), 400)
pdf = np.exp(-0.5 * (grid - analytic_mean) ** 2 / analytic_var) / np.sqrt(2 * np.pi * analytic_var)

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(ou_path.t, ou_path.y[:, 0], color="#b23a48", linewidth=1.4)
axes[0].axhline(mean_level, color="#1d3557", linestyle="--", linewidth=1.2, label="mean level")
axes[0].set_title("One OU path")
axes[0].set_xlabel("t")
axes[0].set_ylabel("X(t)")
axes[0].legend()

axes[1].hist(x, bins=90, density=True, color="#e0a458", alpha=0.76, label="Monte Carlo")
axes[1].plot(grid, pdf, color="#1d3557", linewidth=2.0, label="finite-time normal")
axes[1].set_title("OU terminal law")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"paths = {n_paths:,}, steps = {n_steps}, T = {T}")
print(f"mean: empirical {empirical_mean:.5f}, analytic {analytic_mean:.5f}, abs error {abs(empirical_mean - analytic_mean):.5f}")
print(f"var : empirical {empirical_var:.5f}, finite-time {analytic_var:.5f}, stationary {stationary_var:.5f}")

paths = 80,000, steps = 1200, T = 3.0
mean: empirical 1.48901, analytic 1.48871, abs error 0.00030
var : empirical 0.03427, finite-time 0.03403, stationary 0.03403

Geometric Brownian motion: stock-like long-run behavior

Geometric Brownian motion is a stylized teaching model for multiplicative noise. It is useful for explaining log growth, but it is not a trading recommendation or a realistic market simulator.

gbm = GeometricBrownianMotion(mu=[0.08], sigma=[0.20])
S0 = 100.0
T = 10.0
n_paths = 250_000

z = rng.normal(size=n_paths)
log_growth_rate = gbm.expected_log_growth()[0]
log_terminal = np.log(S0) + log_growth_rate * T + gbm.sigma[0] * np.sqrt(T) * z
terminal = np.exp(log_terminal)
annualized_log_returns = np.log(terminal / S0) / T

analytic_mean = gbm.mean([S0], T)[0]
empirical_mean = float(np.mean(terminal))
prob_loss = float(np.mean(terminal < S0))
quantiles = np.quantile(terminal, [0.05, 0.25, 0.50, 0.75, 0.95])

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].hist(terminal, bins=120, density=True, color="#6b8e9f", alpha=0.75)
axes[0].axvline(S0, color="#b23a48", linestyle="--", linewidth=1.4, label="initial price")
axes[0].set_title("GBM terminal values")
axes[0].set_xlabel("S(T)")
axes[0].set_xlim(0, np.quantile(terminal, 0.995))
axes[0].legend()

axes[1].hist(annualized_log_returns, bins=100, density=True, color="#7aa95c", alpha=0.75)
axes[1].axvline(log_growth_rate, color="#1d3557", linewidth=2.0, label="E[d log S]/dt")
axes[1].set_title("Annualized log growth")
axes[1].set_xlabel("log(S(T)/S0) / T")
axes[1].legend()
plt.tight_layout()

print(f"paths = {n_paths:,}, horizon = {T:g} years")
print(f"expected log growth = {log_growth_rate:.4f} per year")
print(f"mean terminal: empirical {empirical_mean:.2f}, analytic {analytic_mean:.2f}")
print(f"P[S(T) < S0] = {prob_loss:.3f}")
print("terminal quantiles 5/25/50/75/95% =", np.round(quantiles, 2))

paths = 250,000, horizon = 10 years
expected log growth = 0.0600 per year
mean terminal: empirical 222.55, analytic 222.55
P[S(T) < S0] = 0.171
terminal quantiles 5/25/50/75/95% = [ 64.56 119.13 182.27 279.15 516.16]

What this validates

• pyna SDE models produce reproducible single-path Trajectory objects.

• Brownian and OU empirical terminal laws match their analytic means and variances within Monte Carlo error.

• GBM exposes the difference between arithmetic mean growth and long-time log growth.

For production Monte Carlo, keep the model boundary in pyna and move the ensemble kernel to a vectorized, parallel or accelerator-backed implementation. Do not promote an ensemble to Cycle, Tube or IslandChain unless the geometry claim is explicit and numerically justified.