如何在Python中解决/适应几何布朗运动过程？

Question

例如，下面的代码模拟几何布朗运动（GBM）过程，该过程满足以下随机微分方程：

这段代码是维基百科文章中的代码的简化版本。

import numpy as np
np.random.seed(1)

def gbm(mu=1, sigma = 0.6, x0=100, n=50, dt=0.1):
    step = np.exp( (mu - sigma**2 / 2) * dt ) * np.exp( sigma * np.random.normal(0, np.sqrt(dt), (1, n)))
    return x0 * step.cumprod()

series = gbm()

如何在Python中适应GBM过程？也就是说，如何估计mu和sigma并解决随机微分方程，给定时间序列series？

- Greg

---

我不太理解这里的物理问题，但是对于拟合参数，你可以尝试使用 scipy.optimize.curve_fit。 - Gerges

您可以使用许多过程的实现来计算其统计矩。这些矩将与μ和σ相关联，但我不确定如何。它们的名称相当暗示了如何做到这一点。 - kevinkayaks

你不能简单地取对数，进行线性拟合以获得mu-sigma^2/2和一些截距，然后减去线性拟合来估计sigma吗？ - Paul Panzer

你可能会对这个感兴趣：https://symfit.readthedocs.io/en/stable/fitting_types.html#ode-fitting这个链接使用了我编写的 symfit 包，使得在 Python 中处理这样的拟合过程变得更加容易。 - tBuLi

看这个方程，我有一种感觉，从你的时间序列（St和dSt）中构建Wt可能更容易，并将其设置为mu和sigma的函数。然后，您可以使用优化算法来拟合sigma和mu，以便Wt能够再现预期的统计分布。 - Tarifazo

Answer 1

SDE的参数估计是一个研究级领域，因此相对比较复杂。有关该主题的整本书存在于世上。如需更多详细信息，请随意查看那些书籍。

但这里有一个简单的方法来解决这个问题。首先，注意到对数GBM是一个仿射变换的维纳过程（即线性Ito漂移扩散过程）。所以

d ln(S_t) = (mu - sigma^2 / 2) dt + sigma dB_t

因此我们可以估计对数过程的参数并将其转换为适合原始过程的参数。例如，可以查看[1]、[2]、[3]和[4]等链接。

以下是一段脚本，它以两种简单的方式对漂移进行了估计（只是想看到它们之间的差异），并为扩散估计了一种方法（抱歉只有一种方法）。对数过程的漂移通过(X_T - X_0) / T 和增量最大似然估计法进行估计（请参见代码）。扩散参数以其定义的微小方差的形式被（有偏地）估计。

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

- user3658307