4

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

在此处输入图像描述

该代码是此 Wikipedia 文章中代码的压缩版本。

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？

4

1 回答 1

8

SDE 的参数估计是一个研究级别的领域，因此相当重要。关于这个主题的整本书都有。随意查看这些以获取更多详细信息。

但是对于这种情况，这是一种简单的方法。首先，注意GBM 的对数是一个仿射变换的维纳过程（即线性Ito 漂移-扩散过程）。所以

d ln(S_t) = (mu - sigma^2 / 2) dt + sigma 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
于 2019-01-13T21:20:50.630 回答