Estimations de distribution Monte Carlo pour SDE (EDS)

Ce notebook est une étude de cas pratique de SDE (équations différentielles stochastiques, EDS) pour pyna. Il utilise les classes SDE de pyna pour la frontière du modèle et la géométrie d’une trajectoire unique, puis des ensembles NumPy vectorisés pour les estimations de distribution.

Le notebook est volontairement préexécuté localement. Les sorties sauvegardées sont rendues par Sphinx/nbsphinx sur GitHub Pages, de sorte que le workflow de documentation n’a pas besoin de consacrer du temps CPU CI à des échantillonnages Monte Carlo répétés.

%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("graine rng = 20260701")

graine rng = 20260701

Mouvement brownien : distribution terminale

Une réalisation brownienne unique est une Trajectory. Un grand ensemble estime une distribution de probabilité. Pour le mouvement brownien, la distribution terminale est analytique ; c’est donc un contrôle de régression utile pour le code Monte Carlo.

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_trajectoires = 60_000
n_pas = 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_trajectoires = 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)

analytique_mean = bm.mean([0.0], T)[0]
analytique_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 - analytique_mean) ** 2 / analytique_var) / np.sqrt(2 * np.pi * analytique_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="une Trajectory pyna")
axes[0].set_title("trajectoires browniennes exemples")
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="normale analytique")
axes[1].set_title("Distribution terminale X(T)")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"trajectoires = {n_paths:,}, pas = {n_steps}, T = {T}")
print(f"moyenne : empirique {empirical_mean:.5f}, analytique {analytique_mean:.5f}, erreur absolue {abs(empirical_mean - analytique_mean):.5f}")
print(f"var : empirique {empirical_var:.5f}, analytique {analytique_var:.5f}, erreur absolue {abs(empirical_var - analytique_var):.5f}")

trajectoires = 60,000, pas = 600, T = 2.0
moyenne : empirique 0.29974, analytique 0.30000, erreur absolue 0.00026
var : empirique 2.00725, analytique 2.00000, erreur absolue 0.00725

Ornstein-Uhlenbeck : approche d’une loi stationnaire

Le processus OU est un modèle pédagogique compact de retour vers la moyenne. Il montre aussi une ItoSDE personnalisée avec dérive et diffusion définies par l’utilisateur. L’ensemble Monte Carlo vérifie la loi normale à temps fini et la variance stationnaire.

theta = 1.8
mean_level = 1.5
sigma = 0.35
x0 = -1.0
T = 3.0
n_pas = 1200
dt = T / n_steps
n_trajectoires = 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)

analytique_mean = mean_level + (x0 - mean_level) * np.exp(-theta * T)
analytique_var = sigma**2 / (2 * theta) * (1 - np.exp(-2 * theta * T))
stationnaire_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 - analytique_mean) ** 2 / analytique_var) / np.sqrt(2 * np.pi * analytique_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="niveau moyen")
axes[0].set_title("Une trajectoire OU")
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="normale à temps fini")
axes[1].set_title("loi terminale OU")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"trajectoires = {n_paths:,}, pas = {n_steps}, T = {T}")
print(f"moyenne : empirique {empirical_mean:.5f}, analytique {analytique_mean:.5f}, erreur absolue {abs(empirical_mean - analytique_mean):.5f}")
print(f"var : empirique {empirical_var:.5f}, temps fini {analytique_var:.5f}, stationnaire {stationnaire_var:.5f}")

trajectoires = 80,000, pas = 1200, T = 3.0
moyenne : empirique 1.48901, analytique 1.48871, erreur absolue 0.00030
var : empirique 0.03427, temps fini 0.03403, stationnaire 0.03403

Mouvement brownien géométrique : comportement à long terme de type action

Le mouvement brownien géométrique est un modèle pédagogique stylisé pour le bruit multiplicatif. Il est utile pour expliquer la croissance logarithmique, mais ce n’est ni une recommandation de trading ni un simulateur de marché réaliste.

gbm = GeometricBrownianMotion(mu=[0.08], sigma=[0.20])
S0 = 100.0
T = 10.0
n_trajectoires = 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

analytique_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="prix initial")
axes[0].set_title("valeurs terminales GBM")
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("croissance logarithmique annualisée")
axes[1].set_xlabel("log(S(T)/S0) / T")
axes[1].legend()
plt.tight_layout()

print(f"trajectoires = {n_paths:,}, horizon = {T:g} ans")
print(f"croissance log attendue = {log_growth_rate:.4f} par an")
print(f"moyenne terminale : empirique {empirical_mean:.2f}, analytique {analytique_mean:.2f}")
print(f"P[S(T) < S0] = {prob_loss:.3f}")
print("quantiles terminaux 5/25/50/75/95 % =", np.round(quantiles, 2))

trajectoires = 250,000, horizon = 10 ans
croissance log attendue = 0.0600 par an
moyenne terminale : empirique 222.55, analytique 222.55
P[S(T) < S0] = 0.171
quantiles terminaux 5/25/50/75/95 % = [ 64.56 119.13 182.27 279.15 516.16]

Ce que cela valide

• Les modèles SDE de pyna produisent des objets Trajectory à trajectoire unique reproductibles.

• Les lois terminales empiriques brownienne et OU correspondent à leurs moyennes et variances analytiques dans l’erreur Monte Carlo.

• Le GBM met en évidence la différence entre croissance moyenne arithmétique et croissance logarithmique à long terme.

Pour un Monte Carlo de production, gardez la frontière du modèle dans pyna et déplacez le noyau d’ensemble vers une implémentation vectorisée, parallèle ou accélérée. Ne promouvez un ensemble en Cycle, Tube ou IslandChain que si l’assertion géométrique est explicite et numériquement justifiée.