Оценки распределений SDE методом Монте-Карло#

Этот ноутбук - практический пример SDE в pyna. Он использует классы SDE из pyna для описания модели и геометрии одиночной траектории, а затем применяет векторизованные ансамбли NumPy для оценок распределений.

Ноутбук намеренно предварительно выполнен локально. Сохраненные выводы отображаются Sphinx/nbsphinx на GitHub Pages, поэтому рабочий процесс документации не тратит CPU-время CI на повторную выборку Монте-Карло.

[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")

seed rng = 20260701

Броуновское движение: распределение в конечный момент времени#

Одна реализация броуновского движения является Trajectory. Большой ансамбль оценивает вероятностное распределение. Для броуновского движения распределение в конечный момент времени известно аналитически, поэтому это полезная регрессионная проверка кода Монте-Карло.

[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_путей = 60_000
n_шагов = 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_путей = 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("Выборочные траектории броуновского движения")
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("Терминальное распределение X(T)")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"путей = {n_paths:,}, шагов = {n_steps}, T = {T}")
print(f"среднее: эмпирическое {empirical_mean:.5f}, аналитическое {analytic_mean:.5f}, абс. ошибка {abs(empirical_mean - analytic_mean):.5f}")
print(f"дисп.   : эмпирическая {empirical_var:.5f}, аналитическое {analytic_var:.5f}, абс. ошибка {abs(empirical_var - analytic_var):.5f}")

путей = 60,000, шагов = 600, T = 2.0
среднее: эмпирическое 0.29974, аналитическое 0.30000, абс. ошибка 0.00026
дисп.   : эмпирическая 2.00725, аналитическая 2.00000, абс. ошибка 0.00725

../../../../_images/notebooks_i18n_ru_tutorials_sde_monte_carlo_distribution_3_1.svg

Орнштейн-Уленбек: приближение к стационарному закону#

Процесс OU - компактная учебная модель возврата к среднему. Он также демонстрирует пользовательский ItoSDE с заданными drift и diffusion. Ансамбль Монте-Карло проверяет нормальный закон за конечное время и стационарную дисперсию.

[3]:

theta = 1.8
mean_level = 1.5
sigma = 0.35
x0 = -1.0
T = 3.0
n_шагов = 1200
dt = T / n_steps
n_путей = 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="уровень среднего")
axes[0].set_title("Одна траектория 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="finite-time normal")
axes[1].set_title("Терминальный закон OU")
axes[1].set_xlabel("X(T)")
axes[1].legend()
plt.tight_layout()

print(f"путей = {n_paths:,}, шагов = {n_steps}, T = {T}")
print(f"среднее: эмпирическое {empirical_mean:.5f}, аналитическое {analytic_mean:.5f}, абс. ошибка {abs(empirical_mean - analytic_mean):.5f}")
print(f"дисп.   : эмпирическая {empirical_var:.5f}, за конечное время {analytic_var:.5f}, стационарная {stationary_var:.5f}")

путей = 80,000, шагов = 1200, T = 3.0
среднее: эмпирическое 1.48901, аналитическое 1.48871, абс. ошибка 0.00030
дисп.   : эмпирическая 0.03427, за конечное время 0.03403, стационарная 0.03403

../../../../_images/notebooks_i18n_ru_tutorials_sde_monte_carlo_distribution_5_1.svg

Геометрическое броуновское движение: долгосрочное поведение, похожее на цену актива#

Геометрическое броуновское движение - стилизованная учебная модель мультипликативного шума. Она полезна для объяснения логарифмического роста, но не является торговой рекомендацией или реалистичным симулятором рынка.

[4]:

gbm = GeometricBrownianMotion(mu=[0.08], sigma=[0.20])
S0 = 100.0
T = 10.0
n_путей = 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")
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("Годовой логарифмический рост")
axes[1].set_xlabel("log(S(T)/S0) / T")
axes[1].legend()
plt.tight_layout()

print(f"путей = {n_paths:,}, горизонт = {T:g} лет")
print(f"ожидаемый логарифмический рост = {log_growth_rate:.4f} в год")
print(f"терминальное среднее: эмпирическое {empirical_mean:.2f}, аналитическое {analytic_mean:.2f}")
print(f"P[S(T) < S0] = {prob_loss:.3f}")
print("терминальные квантили 5/25/50/75/95% =", np.round(quantiles, 2))

путей = 250,000, горизонт = 10 лет
ожидаемый логарифмический рост = 0.0600 в год
терминальное среднее: эмпирическое 222.55, аналитическое 222.55
P[S(T) < S0] = 0.171
терминальные квантили 5/25/50/75/95% = [ 64.56 119.13 182.27 279.15 516.16]

../../../../_images/notebooks_i18n_ru_tutorials_sde_monte_carlo_distribution_7_1.svg

Что это валидирует#

Модели SDE в pyna создают воспроизводимые объекты одиночных путей Trajectory.
Эмпирические законы в конечный момент времени для броуновского движения и OU совпадают с аналитическими средними и дисперсиями в пределах ошибки Монте-Карло.
GBM показывает различие между ростом арифметического среднего и долгосрочным логарифмическим ростом.

Для производственного Монте-Карло держите описание модели в pyna, а ядро ансамбля переносите в векторизованную, параллельную или аппаратно ускоренную реализацию. Не преобразуйте ансамбль в Cycle, Tube или IslandChain, если геометрическое утверждение не является явным и численно обоснованным.