模拟布朗运动与几何布朗运动

1. 模拟标准布朗运动

1.1 模拟原理

我们通常用 WtW_tWt 或 BtB_tBt 来表示标准布朗运动（Standard Brownian Motion），它是随机过程/随机代数的重要基础。布朗运动满足

B0=0B_0 = 0B0=0；
Bt−BsB_t-B_sBt−Bs 独立于 BrB_rBr，t≥s≥rt \geq s \geq rt≥s≥r；
Bt−BsB_t - B_sBt−Bs 服从均值为 000 方差为 t−st-st−s 的正态分布；
函数 t↦Btt \mapsto B_tt↦Bt 是连续的。

我们利用 Bt−Bs∼N(0,t−s)B_t - B_s \sim N(0, t-s)Bt−Bs∼N(0,t−s)这一点，通过对其离散化，来模拟标准布朗运动。假设在 t∈[0,1]t\in [0,1]t∈[0,1]内均匀取NNN个点，令 Δt=1/N\Delta t = 1/NΔt=1/N，创建这样一个时间网格（grid）：
0,Δt,2Δt,3Δt,⋯,(N−1)Δt,1,0, \Delta t, 2\Delta t, 3\Delta t,\cdots, (N-1)\Delta t, 1, 0,Δt,2Δt,3Δt,⋯,(N−1)Δt,1,那么对 k=0,1,⋯,N−1k = 0, 1, \cdots, N-1k=0,1,⋯,N−1，都有 B(k+1)Δt−BkΔt∼N(0,Δt)B_{(k+1)\Delta t} - B_{k\Delta t} \sim N(0, \Delta t)B(k+1)Δt−BkΔt∼N(0,Δt)，即
B(k+1)Δt−BkΔtΔt∼N(0,1).\frac{B_{(k+1)\Delta t} - B_{k\Delta t}}{\sqrt{\Delta t}} \sim N(0, 1). ΔtB(k+1)Δt−BkΔt∼N(0,1).假设我们可以生成服从标准正态分布的随机数，将上述关系离散化，即可得到
B(k+1)Δt=BkΔt+ΔtNk,B_{(k+1)\Delta t} = B_{k\Delta t} + \sqrt{\Delta t} N_k, B(k+1)Δt=BkΔt+ΔtNk,其中 NkN_kNk 即为在生成 B(k+1)ΔtB_{(k+1)\Delta t}B(k+1)Δt 时随机生成的服从标准正态分布的随机数。重复利用上式，即可模拟出标准布朗运动的路径。这种模拟方法即为欧拉法（Euler method）。

1.2 模拟代码

根据上述模拟原理，写出对应Python代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as pltdef SBM(T=1, N=100, seed=None):'''Simulate a standard Brownian Motion, using formula$B_{(k+1)\Delta t} = B_{k\Delta t} + \sqrt{\Delta t} N_k$,where N_k is random number drawn from standard normal distributionInput:------T, time interval is [0, T]N, number of sample point in [0, 1], \Delta t = 1 / Nseed: random seedOutput:------points of standard Brownian MotionExample:------>>> y = SBM(T=1, N=100, seed=13)'''if seed:np.random.seed(seed)Normal = lambda : np.random.randn() # generate random number distributed by standard normal distributiondelta_t = 1 / Ns_delta_t = np.sqrt(delta_t)res = np.zeros(shape=(T * N + 1, ))for i in range(T * N):res[i+1] = res[i] + s_delta_t * Normal()return res

如，每单位区间内撒 N=200N=200N=200 个点，模拟区间 t∈[0,2]t \in [0, 2]t∈[0,2] 上的标准布朗运动：

T = 2
N = 200plt.figure(figsize=(20, 12))
plt.plot(np.arange(T*N+1)/N, SBM(T, N, seed=13))
plt.xlabel('t')
plt.ylabel('$B_t$')
plt.title('Simulate standard Brownian Motion, $t\in [0, 2], \Delta t = \\frac{1}{200}$')
plt.show()

模拟图：

1.3 理论验证

下面我们验证上述模拟方式是合理的。

我们知道
P{max⁡0≤t≤1Bt>12}=2P{B1>12}=2[1−Φ(12)],\mathbb P\left\{\max_{0 \leq t \leq 1} B_t > \frac{1}{2}\right\} = 2 \mathbb P\left\{ B_1 > \frac{1}{2} \right\} = 2\left[ 1 - \Phi\left(\frac{1}{2}\right) \right], P{0≤t≤1maxBt>21}=2P{B1>21}=2[1−Φ(21)],具体值可用如下代码计算：
```
import scipy as spPhi = lambda x: sp.stats.norm.cdf(x)
print(2 * (1 - Phi(1/2))
```
0.6170750774519738

我们通过变换 seed 来生成 100010001000 次模拟，并计算对应概率：
```
times = 1000
ls = []
for i in range(times):ls.append(SBM(T=1, N=500, seed=i))cnt1 = 0
for bm in ls:if max(bm > 1/2):cnt1 += 1print(cnt1 / 1000)
```
0.615

可以看到，模拟概率与理论概率十分相近。

我们知道
P{B1>0,B2<0}=18,\mathbb P \{B_1 > 0, B_2 < 0\} = \frac{1}{8}, P{B1>0,B2<0}=81,同理，生成 100010001000 次模拟并计算概率：

times = 1000
ls = []
for i in range(times):ls.append(SBM(T=2, N=1000, seed=i))cnt2 = 0
for bm in ls:if bm[1000] > 0 and bm[1000*2] < 0:cnt2 += 1print(cnt2 / 1000)

0.129

2. 模拟布朗运动

2.1 模拟原理

一般的布朗运动 XtX_tXt 满足
dXt=m(t,Xt)dt+σ(t,Xt)dBt,X0=0.dX_t = m(t, X_t) dt + \sigma(t, X_t) dB_t, \quad X_0 = 0. dXt=m(t,Xt)dt+σ(t,Xt)dBt,X0=0.这里我们简便起见，假设 m(t,Xt)=m,σ(t,Xt)=σm(t,X_t) = m, \sigma(t, X_t) = \sigmam(t,Xt)=m,σ(t,Xt)=σ。同样地，在单位区间 [0,1][0, 1][0,1] 内均匀撒 NNN 个点，取 Δt=1/N\Delta t = 1/NΔt=1/N，利用性质
P{Xt+Δt−Xt=σΔt∣Xt}=12[1+mσΔt],P{Xt+Δt−Xt=−σΔt∣Xt}=12[1−mσΔt],\begin{aligned} &P\{X_{t+\Delta t}-X_t=\sigma \sqrt{\Delta t} \mid X_t\}=\frac{1}{2}\left[1+\frac{m}{\sigma} \sqrt{\Delta t}\right], \\ &P\{X_{t+\Delta t}-X_t=-\sigma \sqrt{\Delta t} \mid X_t\}=\frac{1}{2}\left[1-\frac{m}{\sigma} \sqrt{\Delta t}\right], \end{aligned} P{Xt+Δt−Xt=σΔt∣Xt}=21[1+σmΔt],P{Xt+Δt−Xt=−σΔt∣Xt}=21[1−σmΔt],采用二项抽样法（Binomial Sampling Schemes）来模拟。即，给定 XkΔtX_{k\Delta t}XkΔt，它在未来的 Δt\Delta tΔt 时间内以概率为 12[1+mσΔt]\frac{1}{2}\left[1+\frac{m}{\sigma} \sqrt{\Delta t}\right]21[1+σmΔt] 向上走 σΔt\sigma \sqrt{\Delta t}σΔt 个单位，以概率为 12[1−mσΔt]\frac{1}{2}\left[1-\frac{m}{\sigma} \sqrt{\Delta t}\right]21[1−σmΔt] 向下走 σΔt\sigma \sqrt{\Delta t}σΔt 个单位。代入上述性质即得到对 k=0,1,⋯,N−1k=0, 1, \cdots, N-1k=0,1,⋯,N−1，有
P{X(k+1)Δt=XkΔt+σΔt∣XkΔt}=12[1+mσΔt],P{X(k+1)Δt=XkΔt−σΔt∣XkΔt}=12[1−mσΔt].\begin{aligned} &P\left\{X_{(k+1)\Delta t}=X_{k\Delta t} + \sigma \sqrt{\Delta t} \mid X_{k \Delta t}\right\}=\frac{1}{2}\left[1+\frac{m}{\sigma} \sqrt{\Delta t}\right], \\ &P\left\{X_{(k+1)\Delta t}=X_{k\Delta t} - \sigma \sqrt{\Delta t} \mid X_{k \Delta t}\right\}=\frac{1}{2}\left[1-\frac{m}{\sigma} \sqrt{\Delta t}\right]. \end{aligned} P{X(k+1)Δt=XkΔt+σΔt∣XkΔt}=21[1+σmΔt],P{X(k+1)Δt=XkΔt−σΔt∣XkΔt}=21[1−σmΔt].

2.2 模拟代码

基于上述二项抽样，给出Python代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as pltdef BM_b(m=1, sigma=1, T=1, N=100, seed=None):r'''Simulate a Brownian Motion, using binomial schemes:\begin{aligned}&P\{X_{t+\Delta t}-X_t=\sigma \sqrt{\Delta t} \mid X_t\}=\frac{1}{2}\left[1+\frac{m}{\sigma} \sqrt{\Delta t}\right], \\&P\{X_{t+\Delta t}-X_t=-\sigma \sqrt{\Delta t} \mid X_t\}=\frac{1}{2}\left[1-\frac{m}{\sigma} \sqrt{\Delta t}\right].\end{aligned}The Brownian motion satisfies dynamic$dX_t = m dt + \sigma dB_t$,with drift m, variance \sigma^2 constant, X_0 = 0Input:------m, driftsigma, square root of variance parameterT, time interval is [0, T]N, number of sample point in [0, 1], \Delta t = 1 / Nseed: random seedOutput:------points of Brownian MotionExample:------>>> y = BM_b(m=1, sigma=1, T=1, N=100, seed=13)'''if seed:np.random.seed(seed)dt = 1 / Ns_dt = np.sqrt(dt)P_up = 0.5 * (1 + m / sigma * s_dt) # probability of goes upP_down = 1 - P_up# vector of depending on whether each step goes up or downomega = np.random.choice([1, -1], size=N * T, p=[P_up, P_down]) res = np.cumsum(omega) * sigma * s_dtreturn res # note that this series doesn't include X_0

如，每单位区间内撒 N=1000N=1000N=1000 个点，模拟区间 t∈[0,2]t \in [0, 2]t∈[0,2] 上 m=−3/10,σ=1.5m =-3/\sqrt{10}, \sigma = 1.5m=−3/10,σ=1.5 的布朗运动：

m = -3 / np.sqrt(10)
sigma = 1.5
T = 2
N = 1000Xt = BM_b(m=m, sigma=sigma, T=T, N=N, seed=13)
t = (np.arange(Xt.shape[0]) + 1) / 1000plt.figure(figsize=(20, 12))
plt.plot(t, Xt, label='$X_t$')
plt.xlabel('t', fontsize=20)
#plt.ylabel('$B_t$', fontsize=20)
plt.legend(fontsize=20)
plt.title('Simulate Brownian Motion Using Binomial Schemes, $t\in [0, 2], \Delta t = \\frac{1}{1000}$', fontsize=20)
plt.show()

模拟图：

2.3 理论验证

下面我们验证上述对drift与方差项均为常数的布朗运动模拟是合理的。

对满足
dXt=mdt+σdBtdX_t = mdt + \sigma dB_t dXt=mdt+σdBt的布朗运动，我们有 Xt∼N(mt,σ2t)X_t \sim N(mt, \sigma^2t)Xt∼N(mt,σ2t)。因此，对常数 aaa，我们有
P{Xt>a}=P{Xt−mtσt>a−mtσt}=Φ(mt−aσt).\begin{aligned} &\mathbb P \left\{ X_t > a \right\} \\ = &\mathbb P \left\{ \frac{X_t - mt}{\sigma \sqrt{t}} > \frac{a - mt}{\sigma \sqrt{t}} \right\} \\ = &\Phi \left( \frac{mt - a}{\sigma \sqrt{t}} \right). \end{aligned} ==P{Xt>a}P{σtXt−mt>σta−mt}Φ(σtmt−a).对在2.2节中模拟的布朗运动而言（即 m=−3/10,σ=1.5m =-3/\sqrt{10}, \sigma = 1.5m=−3/10,σ=1.5），我们有
P{X2>0}=Φ(−25),P{X2>1}=Φ(−610−11.52),\mathbb P\{ X_2 > 0 \} = \Phi \left( \frac{-2}{\sqrt{5}} \right), \mathbb P\{ X_2 > 1 \} = \Phi \left( \frac{\frac{-6}{\sqrt{10}} - 1}{1.5 \sqrt{2}} \right), P{X2>0}=Φ(5−2),P{X2>1}=Φ(1.5210−6−1),他们分别可以被计算为：

from scipy.stats import normprint('{:.4f}'.format(norm.cdf(-2 / np.sqrt(5))))
print('{:.4f}'.format(norm.cdf((-6 / np.sqrt(10) - 1) / (1.5 * np.sqrt(2)))))

0.1855
0.0860

我们通过模拟 100001000010000 次，计算出对应概率，并做对比：

m = -3 / np.sqrt(10)
sigma = 1.5cnt1 = 0
cnt2 = 0
nums = 10000for i in range(nums):res = BM_b(m=m, sigma=sigma, T=T, N=N, seed=i)if res[-1] > 0:cnt1 += 1if res[-1] > 1:cnt2 += 1print(cnt1 / nums)
print(cnt2 / nums)

0.1799
0.0871

3. 模拟几何布朗运动

3.1 模拟原理

几何布朗运动（Geometric Brownian Motion）满足
dXt=Xt(mdt+σdBt),X0=1,dX_t = X_t\left(mdt + \sigma dB_t \right), \quad X_0=1, dXt=Xt(mdt+σdBt),X0=1,其对应的解为
Xt=X0e(m−σ22)t+σBt,X_t = X_0 e^{\left( m - \frac{\sigma^2}{2} \right)t + \sigma B_t}, Xt=X0e(m−2σ2)t+σBt,可以通过伊藤微分直接验证。

我们使用欧拉法（Euler method）来对几何布朗运动做模拟。生成单位区间 [0,1][0, 1][0,1] 内 NNN 个均匀的点，取 Δt=1/N\Delta t = 1/NΔt=1/N，对 k=0,1,⋯,N−1k=0, 1, \cdots, N-1k=0,1,⋯,N−1，利用迭代式
X(k+1)Δt=XkΔt(1+mΔt+σΔtNk)X_{(k+1)\Delta t} = X_{k\Delta t} \left( 1 + m \Delta t + \sigma \sqrt{\Delta t} N_k \right) X(k+1)Δt=XkΔt(1+mΔt+σΔtNk)做模拟。

3.2 模拟代码与验证

基于上述欧拉法模拟，给出对应Python代码：

def GBM(T=1, N=100, m=1, sigma=1, seed=None, X_0=1):'''Simulate a Geometric Brownian Motion, using formula$X_{(k+1)\Delta t} = X_{k\Delta t} \left( 1 + m \Delta t + \sigma \sqrt{\Delta t} N_k \right)$,where N_k is random number drawn from standard normal distributionA Geometric Brownian Motion satisfies dynamic$dX_t = X_t(m dt + \sigma dB_t)$Input:------T, time interval is [0, T]N, number of sample point in [0, 1], \Delta t = 1 / Nm, sigma, parameters in dynamicseed: random seedX_0: start point of GBMOutput:------points of Geometric Brownian MotionExample:------y = GBM(T=1, N=1000, m=1, sigma=1, seed=13)'''if seed:np.random.seed(seed)Normal = lambda : np.random.randn() # generate random number distributed by standard normal distributiondelta_t = 1 / Ns_delta_t = np.sqrt(delta_t)res = np.zeros(shape=(T * N + 1, ))res[0] = X_0for i in range(T * N):res[i+1] = res[i] * (1 + m * delta_t + sigma * s_delta_t * Normal())return res

对于 m=1/2,σ=1m=1/2, \sigma=1m=1/2,σ=1 的几何布朗运动 XtX_tXt，即其满足
dXt=12Xtdt+XtdBt,X0=1,dX_t = \frac{1}{2} X_t dt + X_t dB_t, \quad X_0 = 1, dXt=21Xtdt+XtdBt,X0=1,可以计算得到 Xt=eBtX_t = e^{B_t}Xt=eBt。因此，我们指定相同的随机数种子 seed，通过两种方式做模拟：一种即直接通过几何布朗运动做模拟；一种即先生成标准布朗运动，再做指数运算。两种模拟方式理应给出相同的路径图。

Xt = GBM(T=1, N=1000, m=1/2, sigma=1, seed=11)
eBt = np.exp(SBM(T=1, N=1000, seed=11))
t = np.arange(Xt.shape[0]) / 1000plt.figure(figsize=(20, 12))
plt.plot(t, eBt, label='$e^{B_t}$')
plt.plot(t, Xt, label='$X_t$')
plt.xlabel('t', fontsize=20)
#plt.ylabel('$B_t$', fontsize=20)
plt.legend(fontsize=20)
plt.title('Simulate Geometric Brownian Motion, $t\in [0, 1], \Delta t = \\frac{1}{1000}$', fontsize=20)
plt.show()

可以看到，抛除浮点数运算误差，两种不同的模拟方式确实给出了近乎相同的结果。