Übersetzung möglicherweise veraltet Die englische Quellseite ist neuer; diese deutsche Seite ist möglicherweise noch nicht vollständig synchronisiert.
Englische Quelle öffnen

Monte-Carlo-Verteilungsschätzungen für SDE#

Dieses Notebook ist eine praktische SDE-Fallstudie für pyna. Es verwendet die SDE-Klassen von pyna für Modellgrenzen und Einzelpfadgeometrie und nutzt anschließend vektorisierte NumPy-Ensembles für Verteilungsschätzungen.

Das Notebook wird absichtlich lokal vorab ausgeführt. Die gespeicherten Ausgaben werden von Sphinx/nbsphinx auf GitHub Pages gerendert, sodass der Dokumentations-Workflow keine CI-CPU-Zeit für wiederholtes Monte-Carlo-Sampling aufwenden muss.

[1]:

%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

Brownsche Bewegung: Endverteilung#

Eine einzelne brownsche Realisierung ist eine Trajectory. Ein großes Ensemble schätzt eine Wahrscheinlichkeitsverteilung. Für brownsche Bewegung ist die Endverteilung analytisch bekannt und damit ein nützlicher Regressionstest für Monte-Carlo-Code.

[2]:

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_Pfade = 60_000
n_Schritte = 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_Pfade = 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)

analytisch_mean = bm.mean([0.0], T)[0]
analytisch_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 - analytisch_mean) ** 2 / analytisch_var) / np.sqrt(2 * np.pi * analytisch_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="eine pyna-Trajectory")
axes[0].set_title("Brown'sche Beispielpfade")
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="analytische Normalverteilung")
axes[1].set_title("Endverteilung X(T)")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"Pfade = {n_paths:,}, Schritte = {n_steps}, T = {T}")
print(f"Mittelwert: empirisch {empirical_mean:.5f}, analytisch {analytisch_mean:.5f}, absoluter Fehler {abs(empirical_mean - analytisch_mean):.5f}")
print(f"Varianz: empirisch {empirical_var:.5f}, analytisch {analytisch_var:.5f}, absoluter Fehler {abs(empirical_var - analytisch_var):.5f}")

Pfade = 60,000, Schritte = 600, T = 2.0
Mittelwert: empirisch 0.29974, analytisch 0.30000, absoluter Fehler 0.00026
Varianz: empirisch 2.00725, analytisch 2.00000, absoluter Fehler 0.00725
../../../../_images/notebooks_i18n_de_tutorials_sde_monte_carlo_distribution_3_1.svg

Ornstein-Uhlenbeck: Annäherung an ein stationäres Gesetz#

Der OU-Prozess ist ein kompaktes Lehrmodell für Mittelwert-Rückkehr. Er zeigt außerdem eine benutzerdefinierte ItoSDE mit eigener Drift und Diffusion. Das Monte-Carlo-Ensemble prüft das endzeitliche Normalgesetz und die stationäre Varianz.

[3]:

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

analytisch_mean = mean_level + (x0 - mean_level) * np.exp(-theta * T)
analytisch_var = sigma**2 / (2 * theta) * (1 - np.exp(-2 * theta * T))
stationär_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 - analytisch_mean) ** 2 / analytisch_var) / np.sqrt(2 * np.pi * analytisch_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="Mittelwertniveau")
axes[0].set_title("Ein OU-Pfad")
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="Normalverteilung bei endlicher Zeit")
axes[1].set_title("OU-Endgesetz")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"Pfade = {n_paths:,}, Schritte = {n_steps}, T = {T}")
print(f"Mittelwert: empirisch {empirical_mean:.5f}, analytisch {analytisch_mean:.5f}, absoluter Fehler {abs(empirical_mean - analytisch_mean):.5f}")
print(f"Varianz: empirisch {empirical_var:.5f}, endliche Zeit {analytisch_var:.5f}, stationär {stationär_var:.5f}")

Pfade = 80,000, Schritte = 1200, T = 3.0
Mittelwert: empirisch 1.48901, analytisch 1.48871, absoluter Fehler 0.00030
Varianz: empirisch 0.03427, endliche Zeit 0.03403, stationär 0.03403
../../../../_images/notebooks_i18n_de_tutorials_sde_monte_carlo_distribution_5_1.svg

Geometrische brownsche Bewegung: aktienähnliches Langzeitverhalten#

Geometrische brownsche Bewegung ist ein stilisiertes Lehrmodell für multiplikatives Rauschen. Sie eignet sich zur Erklärung logarithmischen Wachstums, ist aber weder eine Handelsempfehlung noch ein realistischer Marktsimulator.

[4]:

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

analytisch_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="Anfangspreis")
axes[0].set_title("GBM-Endwerte")
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("Annualisiertes Log-Wachstum")
axes[1].set_xlabel("log(S(T)/S0) / T")
axes[1].legend()
plt.tight_layout()

print(f"Pfade = {n_paths:,}, Horizont = {T:g} Jahre")
print(f"erwartetes Log-Wachstum = {log_growth_rate:.4f} pro Jahr")
print(f"Endmittelwert: empirisch {empirical_mean:.2f}, analytisch {analytisch_mean:.2f}")
print(f"P[S(T) < S0] = {prob_loss:.3f}")
print("Endquantile 5/25/50/75/95% =", np.round(quantiles, 2))

Pfade = 250,000, Horizont = 10 Jahre
erwartetes Log-Wachstum = 0.0600 pro Jahr
Endmittelwert: empirisch 222.55, analytisch 222.55
P[S(T) < S0] = 0.171
Endquantile 5/25/50/75/95% = [ 64.56 119.13 182.27 279.15 516.16]
../../../../_images/notebooks_i18n_de_tutorials_sde_monte_carlo_distribution_7_1.svg

Was hier validiert wird#

  • pyna-SDE-Modelle erzeugen reproduzierbare Einzelpfad-Trajectory-Objekte.

  • Die empirischen Endgesetze von brownscher Bewegung und OU stimmen innerhalb des Monte-Carlo-Fehlers mit den analytischen Mittelwerten und Varianzen überein.

  • GBM macht den Unterschied zwischen arithmetischem Mittelwachstum und langfristigem logarithmischem Wachstum sichtbar.

Für produktive Monte-Carlo-Simulationen sollte die Modellgrenze in pyna bleiben, während der Ensemble-Kernel in eine vektorisierte, parallele oder beschleunigerbasierte Implementierung ausgelagert wird. Ein Ensemble sollte nur dann zu Cycle, Tube oder IslandChain hochgestuft werden, wenn die geometrische Aussage explizit und numerisch begründet ist.