Как решить/сопоставить процесс геометрического броуновского движения в Python?

Оценка параметров для SDE является областью исследовательского уровня и, следовательно, довольно нетривиальной. На эту тему существуют целые книги. Не стесняйтесь изучать их для получения более подробной информации.

Но вот тривиальный подход для этого случая. Во-первых, обратите внимание, что логарифм GBM представляет собой аффинно преобразованный винеровский процесс (то есть линейный процесс дрейфа-диффузии Ито). Так

d ln (S_t) = (мю - сигма ^ 2 / 2) dt + сигма dB_t

Таким образом, мы можем оценить параметры логарифмического процесса и перевести их в соответствие с исходным процессом. Проверьте [1], [2], [3], [4], например.

Вот скрипт, который делает это двумя простыми способами для дрифта (просто хотел увидеть разницу) и только одним для диффузии (извините). Дрейф лог-процесса оценивается (X_T - X_0) / T и с помощью инкрементного MLE (см. код). Параметр диффузии оценивается (необъективно) с определением его как бесконечно малой дисперсии.

import numpy as np

np.random.seed(9713)

# Parameters
mu = 1.5
sigma = 0.9
x0 = 1.0
n = 1000
dt = 0.05

# Times
T = dt*n
ts = np.linspace(dt, T, n)

# Geometric Brownian motion generator
def gbm(mu, sigma, x0, n, dt):
    step = np.exp( (mu - sigma**2 / 2) * dt ) * np.exp( sigma * np.random.normal(0, np.sqrt(dt), (1, n)))
    return x0 * step.cumprod()

# Estimate mu just from the series end-points
# Note this is for a linear drift-diffusion process, i.e. the log of GBM
def simple_estimate_mu(series):
    return (series[-1] - x0) / T

# Use all the increments combined (maximum likelihood estimator)
# Note this is for a linear drift-diffusion process, i.e. the log of GBM
def incremental_estimate_mu(series):
    total = (1.0 / dt) * (ts**2).sum()
    return (1.0 / total) * (1.0 / dt) * ( ts * series ).sum()

# This just estimates the sigma by its definition as the infinitesimal variance (simple Monte Carlo)
# Note this is for a linear drift-diffusion process, i.e. the log of GBM
# One can do better than this of course (MLE?)
def estimate_sigma(series):
    return np.sqrt( ( np.diff(series)**2 ).sum() / (n * dt) )

# Estimator helper
all_estimates0 = lambda s: (simple_estimate_mu(s), incremental_estimate_mu(s), estimate_sigma(s))

# Since log-GBM is a linear Ito drift-diffusion process (scaled Wiener process with drift), we
# take the log of the realizations, compute mu and sigma, and then translate the mu and sigma
# to that of the GBM (instead of the log-GBM). (For sigma, nothing is required in this simple case).
def gbm_drift(log_mu, log_sigma):
    return log_mu + 0.5 * log_sigma**2

# Translates all the estimates from the log-series
def all_estimates(es):
    lmu1, lmu2, sigma = all_estimates0(es)
    return gbm_drift(lmu1, sigma), gbm_drift(lmu2, sigma), sigma

print('Real Mu:', mu)
print('Real Sigma:', sigma)

### Using one series ###
series = gbm(mu, sigma, x0, n, dt)
log_series = np.log(series)

print('Using 1 series: mu1 = %.2f, mu2 = %.2f, sigma = %.2f' % all_estimates(log_series) )

### Using K series ###
K = 10000
s = [ np.log(gbm(mu, sigma, x0, n, dt)) for i in range(K) ]
e = np.array( [ all_estimates(si) for si in s ] )
avgs = np.mean(e, axis=0)

print('Using %d series: mu1 = %.2f, mu2 = %.2f, sigma = %.2f' % (K, avgs[0], avgs[1], avgs[2]) )

Выход:

Real Mu: 1.5
Real Sigma: 0.9
Using 1 series: mu1 = 1.56, mu2 = 1.54, sigma = 0.96
Using 10000 series: mu1 = 1.51, mu2 = 1.53, sigma = 0.93