接上文：

在這篇文章中，我主要想展示一下weak order 2。0方法應用在幾何布朗運動上的效果。

一般來說，Euler方法的弱收斂階數是1，Milstein方法的弱收斂階數還達不到2。如果我們的目的是估計一個隨機微分方程的解在某個時刻的期望

E[X_T]

，那麼，使用更高階的數值方法就能得到更好的近似效果。

對於幾何布朗運動：

d X_{t}=a X_{t} d t+b X_{t} d W_{t}\\

它的顯式解是：

X_{t}=X_{0} \exp \left(\left(a-\frac{1}{2} b^{2}\right) t+b W_{t}\right)\\

因此

E[X_T]=X_0e^{aT}

。

假設我們每隔

\Delta t

取一系列離散的點

Y_0, Y_1, ..., Y_n

，根據我之前寫的那篇文章，weak order 2。0 的迭代格式是：

\begin{aligned} Y_{n+1}=& Y_{n}+a Y_{n} \Delta t+b Y_{n} \Delta W_{n}+\frac{1}{2} b^{2} Y_{n}\left(\left(\Delta W_{n}\right)^{2}-\Delta t\right) \\ &+a b Y_{n} \Delta Z_{n}+\frac{1}{2} a^{2} Y_{n} \Delta t^{2} \\ &+a b Y_{n}\left(\Delta W_{n} \Delta t-\Delta Z_{n}\right) \end{aligned}\\

其中

\Delta Z_{n}=I_{(1,0)}=\int_{\tau_{n}}^{\tau_{n+1}} \int_{\tau_{n}}^{s} d W_{z} d s=\Delta t \cdot \Delta W

。

Simplified weak order 2。0 的迭代格式是：

\begin{aligned} Y_{n+1}=& Y_{n}+a Y_{n} \Delta t+b Y_{n} \Delta \hat{W}+\frac{1}{2} b^{2} Y_{n}\left((\Delta \hat{W})^{2}-\Delta t\right) \\ &+a b Y_{n} \Delta \hat{W} \Delta t+\frac{1}{2} a^{2} Y_{n} \Delta t^{2} \end{aligned}\\

其中：

P(\Delta \hat{W}=\pm \sqrt{3 \Delta t})=\frac{1}{6}, \quad P(\Delta \hat{W}=0)=\frac{2}{3}

。

因為我們的最終目的是求

E[X_T]

，

本質上布朗運動就是隨機遊走的極限

，因此使用隨機遊走來近似能夠節省計算量。

現在我們要做的，就是驗證上面的格式的弱收斂階數是2。0，即：

\mu(\Delta t)=\left|E\left(g\left(X_{T}\right)\right)-E\left(g\left(Y_{T}^{\Delta}\right)\right)\right| \leq C_{g} \Delta t^{2}\\

從而，理論上來說， #FormatImgID_13# 與 #FormatImgID_14# 擬合出來大約是一條斜率為2的直線。

取

T = 1, a = 1.2, b = 0.5

，檢驗函式

g(x) = x

。

分別寫出兩種方法的路徑模擬函式：

function

path

=

PathWeakGBM

（

a， b， X0， T， delta_t

）

% a sample path for weak order 2。0 GBM

steps

=

0

：

delta_t

：

T

；

path

=

［

X0

］；

for

i

=

steps

（

1

：

end

-

1

）

Y0

=

path

（

end

）；

r

=

randn

（）

*

sqrt

（

delta_t

）；

Y1

=

Y0

+

a

*

Y0

*

delta_t

+

b

*

Y0

*

r

。。。

+

0。5

*

b

^

2

*

Y0

*

（

r

^

2

-

delta_t

）

。。。

+

0。5

*

a

^

2

*

Y0

*

delta_t

^

2

+

a

*

b

*

Y0

*

r

*

delta_t

；

path

=

［

path

Y1

］；

end

end

function

path

=

PathSimplifiedWeakGBM

（

a， b， X0， T， delta_t

）

% a sample path for simplified weak order 2。0 GBM

steps

=

0

：

delta_t

：

T

；

path

=

［

X0

］；

for

i

=

steps

（

1

：

end

-

1

）

Y0

=

path

（

end

）；

r

=

rv

（

1

，

delta_t

）；

Y1

=

Y0

+

a

*

Y0

*

delta_t

+

b

*

Y0

*

r

。。。

+

0。5

*

b

^

2

*

Y0

*

（

r

。^

2

-

delta_t

）

。。。

+

0。5

*

a

^

2

*

Y0

*

delta_t

^

2

+

a

*

b

*

Y0

*

r

*

delta_t

；

path

=

［

path

Y1

］；

end

end

取

N=16000000

，即模擬出非常多的路徑，來計算期望值

E[X_T]

。

a

=

1。2

；

b

=

0。5

；

T

=

1

；

X0

=

1

；

EXT

=

X0

*

exp

（

a

*

T

）；

function

order

=

WeakGBM

（

a， b， X0， T

）

rng

（

2020

）

EXT

=

X0

*

exp

（

a

*

T

）；

N

=

16000000

；

delta

=

2

。^

［

-

1

-

2

-

3

-

4

-

5

］；

err

=

［］；

for

delta_t

=

delta

y

=

X0

*

ones

（

1

，

N

）；

i

=

0

；

while

i

*

delta_t

<

T

r

=

randn

（

1

，

N

）

*

sqrt

（

delta_t

）；

y

=

y

+

a

*

y

*

delta_t

+

b

*

y

。*

r

。。。

+

0。5

*

b

^

2

*

y

。*

（

r

。^

2

-

delta_t

）

。。。

+

0。5

*

a

^

2

*

y

*

delta_t

^

2

+

a

*

b

*

y

。*

r

*

delta_t

；

i

=

i

+

1

；

end

XS

=

mean

（

y

）；

err

=

［

err

abs

（

XS

-

EXT

）］；

end

delta

=

log2

（

delta

）；

err

=

log2

（

err

）；

plot

（

delta

，

err

，

‘o-’

）；

order

=

polyfit

（

delta

，

err

，

1

）；

title

（

‘fitted slope： 1。8368’

）；

xlabel

（

‘log2 -dt’

）

ylabel

（

‘log2 -Weak Errors’

）

end

function

order

=

SimplifiedWeakGBM

（

a， b， X0， T

）

rng

（

2019

）

EXT

=

X0

*

exp

（

a

*

T

）；

N

=

16000000

；

delta

=

2

。^

［

-

1

-

2

-

3

-

4

-

5

-

6

］；

err

=

［］；

for

delta_t

=

delta

y

=

X0

*

ones

（

1

，

N

）；

i

=

0

；

size

=

［

1

，

N

］；

while

i

*

delta_t

<

T

r

=

rv

（

size

，

delta_t

）；

y

=

y

+

a

*

y

*

delta_t

+

b

*

y

。*

r

。。。

+

0。5

*

b

^

2

*

y

。*

（

r

。^

2

-

delta_t

）

。。。

+

0。5

*

a

^

2

*

y

*

delta_t

^

2

+

a

*

b

*

y

。*

r

*

delta_t

；

i

=

i

+

1

；

end

XS

=

mean

（

y

）；

err

=

［

err

abs

（

XS

-

EXT

）］；

end

delta

=

log2

（

delta

）；

err

=

log2

（

err

）；

plot

（

delta

，

err

，

‘o-’

）；

order

=

polyfit

（

delta

，

err

，

1

）；

title

（

‘fitted slope： 1。8179’

）；

xlabel

（

‘log2 -dt’

）

ylabel

（

‘log2 -Weak Errors’

）

end

function

w

=

rv

（

size， delta_t

）

% P（w=\pm \sqrt {3 delta_t）） = 1/6， P（w=0）=2/3

runi

=

rand

（

size

）；

w

=

zeros

（

size

）；

w

（

runi

<

=

1

/

6

）

=

-

sqrt

（

3

*

delta_t

）；

w

（

runi

>

=

5

/

6

）

=

sqrt

（

3

*

delta_t

）；

end

最後執行上面的函式：

o1

=

WeakGBM

（

a

，

b

，

X0

，

T

）；

o2

=

SimplifiedWeakGBM

（

a

，

b

，

X0

，

T

）；

最後得到的結果分別是：

隨機微分方程的模擬（二）

Weak Order 2。0

和：

隨機微分方程的模擬（二）

Simplified Weak Order 2。0

而如果用尤拉方法去做的話，擬合出來的直線斜率差不多應該是1。

總結：本文給出了幾何布朗運動在弱收斂階數2.0下的模擬方法，並設計程式碼驗證了收斂階數。實驗證明，更高階的方法模擬出來的效果更好。

不過，值得一提的是，並不是所有的隨機微分方程用高階方法去模擬，得到的效果都更好，因為高階方法會大量涉及求導，如果方程是高維非線性的，那麼實際還要去估計導數，會導致效果反而變差。