1. 编程号首页

模拟布朗运动与几何布朗运动的区别_布朗运动有哪些例子

编程小号 • • 未分类

模拟布朗运动与几何布朗运动的区别_布朗运动有哪些例子模拟标准布朗运动，布朗运动，集合布朗运动_standardbrownianmotion

1. 模拟标准布朗运动

1.1 模拟原理

我们通常用 W t W_t Wt​ 或 B t B_t Bt​ 来表示标准布朗运动（Standard Brownian Motion），它是随机过程/随机代数的重要基础。布朗运动满足

  • B 0 = 0 B_0 = 0 B0​=0；
  • B t − B s B_t-B_s Bt​−Bs​ 独立于 B r B_r Br​， t ≥ s ≥ r t \geq s \geq r t≥s≥r；
  • B t − B s B_t – B_s Bt​−Bs​ 服从均值为 0 0 0 方差为 t − s t-s t−s 的正态分布；
  • 函数 t ↦ B t t \mapsto B_t t↦Bt​ 是连续的。

我们利用 B t − B s ∼ 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]内均匀取 N N N个点，令 Δ 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 − 1 k = 0, 1, \cdots, N-1 k=0,1,⋯,N−1，都有 B ( k + 1 ) Δ t − B k Δ 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 − B k Δ t Δ t ∼ N ( 0 , 1 ) . \frac{B_{(k+1)\Delta t} – B_{k\Delta t}}{\sqrt{\Delta t}} \sim N(0, 1). Δt
​B(k+1)Δt​−BkΔt​​∼N(0,1).假设我们可以生成服从标准正态分布的随机数，将上述关系离散化，即可得到
B ( k + 1 ) Δ t = B k Δ t + Δ t N k , B_{(k+1)\Delta t} = B_{k\Delta t} + \sqrt{\Delta t} N_k, B(k+1)Δt​=BkΔt​+Δt
​Nk​,其中 N k N_k Nk​ 即为在生成 B ( k + 1 ) Δ t B_{(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 plt



def 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 distribution Input: ------ T, time interval is [0, T] N, number of sample point in [0, 1], \Delta t = 1 / N seed: random seed Output: ------ points of standard Brownian Motion Example: ------ >>> 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 distribution

    delta_t = 1 / N
    s_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 = 200 N=200 N=200 个点，模拟区间 t ∈ [ 0 , 2 ] t \in [0, 2] t∈[0,2] 上的标准布朗运动：

T = 2
N = 200

plt.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()

模拟图：

SBM

1.3 理论验证

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

  • 我们知道
    P { max ⁡ 0 ≤ t ≤ 1 B t > 1 2 } = 2 P { B 1 > 1 2 } = 2 [ 1 − Φ ( 1 2 ) ] , \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≤1max​Bt​>21​}    =    2P{
    B1​>21​}    =    2[1−Φ(21​)],    具体值可用如下代码计算：

    import scipy as sp

Phi = lambda x: sp.stats.norm.cdf(x)
print(2 * (1 - Phi(1/2))

    0.6170750774519738

    我们通过变换 seed 来生成 1000 1000 1000 次模拟，并计算对应概率：

    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 += 1

print(cnt1 / 1000)

    0.615

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

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

    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 += 1

print(cnt2 / 1000)

    0.129

2. 模拟布朗运动

2.1 模拟原理

一般的布朗运动 X t X_t Xt​ 满足
d X t = m ( t , X t ) d t + σ ( t , X t ) d B t , X 0 = 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 , X t ) = m , σ ( t , X t ) = σ m(t,X_t) = m, \sigma(t, X_t) = \sigma m(t,Xt​)=m,σ(t,Xt​)=σ。同样地，在单位区间 [ 0 , 1 ] [0, 1] [0,1] 内均匀撒 N N N 个点，取 Δ t = 1 / N \Delta t = 1/N Δt=1/N，利用性质
P { X t + Δ t − X t = σ Δ t ∣ X t } = 1 2 [ 1 + m σ Δ t ] , P { X t + Δ t − X t = − σ Δ t ∣ X t } = 1 2 [ 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）来模拟。即，给定 X k Δ t X_{k\Delta t} XkΔt​，它在未来的 Δ t \Delta t Δt 时间内以概率为 1 2 [ 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
​ 个单位，以概率为 1 2 [ 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 − 1 k=0, 1, \cdots, N-1 k=0,1,⋯,N−1，有
P { X ( k + 1 ) Δ t = X k Δ t + σ Δ t ∣ X k Δ t } = 1 2 [ 1 + m σ Δ t ] , P { X ( k + 1 ) Δ t = X k Δ t − σ Δ t ∣ X k Δ t } = 1 2 [ 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 plt



def 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 = 0 Input: ------ m, drift sigma, square root of variance parameter T, time interval is [0, T] N, number of sample point in [0, 1], \Delta t = 1 / N seed: random seed Output: ------ points of Brownian Motion Example: ------ >>> y = BM_b(m=1, sigma=1, T=1, N=100, seed=13) '''

    if seed:
        np.random.seed(seed)

    dt = 1 / N
    s_dt = np.sqrt(dt)
    P_up = 0.5 * (1 + m / sigma * s_dt) # probability of goes up
    P_down = 1 - P_up

    # vector of depending on whether each step goes up or down
    omega = np.random.choice([1, -1], size=N * T, p=[P_up, P_down]) 
    res = np.cumsum(omega) * sigma * s_dt
    
    return res # note that this series doesn't include X_0

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

m = -3 / np.sqrt(10)
sigma = 1.5
T = 2
N = 1000

Xt = BM_b(m=m, sigma=sigma, T=T, N=N, seed=13)
t = (np.arange(Xt.shape[0]) + 1) / 1000

plt.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()

模拟图：

bm_b

2.3 理论验证

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

对满足
d X t = m d t + σ d B t dX_t = mdt + \sigma dB_t dXt​=mdt+σdBt​的布朗运动，我们有 X t ∼ N ( m t , σ 2 t ) X_t \sim N(mt, \sigma^2t) Xt​∼N(mt,σ2t)。因此，对常数 a a a，我们有
P { X t > a } = P { X t − m t σ t > a − m t σ t } = Φ ( m t − 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{
σt
​Xt​−mt​>σt
​a−mt​}Φ(σt
​mt−a​).​对在2.2节中模拟的布朗运动而言（即 m = − 3 / 10 , σ = 1.5 m =-3/\sqrt{10}, \sigma = 1.5 m=−3/10
​,σ=1.5），我们有
P { X 2 > 0 } = Φ ( − 2 5 ) , P { X 2 > 1 } = Φ ( − 6 10 − 1 1.5 2 ) , \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.52
​10
​−6​−1​),他们分别可以被计算为：

from scipy.stats import norm

print('{:.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

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

m = -3 / np.sqrt(10)
sigma = 1.5

cnt1 = 0
cnt2 = 0
nums = 10000


for i in range(nums):

    res = BM_b(m=m, sigma=sigma, T=T, N=N, seed=i)

    if res[-1] > 0:
        cnt1 += 1
    
    if res[-1] > 1:
        cnt2 += 1


print(cnt1 / nums)
print(cnt2 / nums)

0.1799
0.0871

3. 模拟几何布朗运动

3.1 模拟原理

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

我们使用欧拉法（Euler method）来对几何布朗运动做模拟。生成单位区间 [ 0 , 1 ] [0, 1] [0,1] 内 N N N 个均匀的点，取 Δ t = 1 / N \Delta t = 1/N Δt=1/N，对 k = 0 , 1 , ⋯   , N − 1 k=0, 1, \cdots, N-1 k=0,1,⋯,N−1，利用迭代式
X ( k + 1 ) Δ t = X k Δ t ( 1 + m Δ t + σ Δ t N k ) 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+σΔt
​Nk​)做模拟。

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 distribution A 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 / N m, sigma, parameters in dynamic seed: random seed X_0: start point of GBM Output: ------ points of Geometric Brownian Motion Example: ------ 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 distribution

    delta_t = 1 / N
    s_delta_t = np.sqrt(delta_t)

    res = np.zeros(shape=(T * N + 1, ))
    res[0] = X_0
    for i in range(T * N):
        res[i+1] = res[i] * (1 + m * delta_t + sigma * s_delta_t * Normal())
    
    return res

对于 m = 1 / 2 , σ = 1 m=1/2, \sigma=1 m=1/2,σ=1 的几何布朗运动 X t X_t Xt​，即其满足
d X t = 1 2 X t d t + X t d B t , X 0 = 1 , dX_t = \frac{1}{2} X_t dt + X_t dB_t, \quad X_0 = 1, dXt​=21​Xt​dt+Xt​dBt​,X0​=1,可以计算得到 X t = e B t X_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]) / 1000

plt.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()

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

gbm

今天的文章模拟布朗运动与几何布朗运动的区别_布朗运动有哪些例子分享到此就结束了，感谢您的阅读。

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。
如需转载请保留出处：https://bianchenghao.cn/86226.html

赞 (0)
编程小号编程小号
0

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注