正向隨機微分方程的經典數值格式模擬

• • 前言
  • 一：隨機微分方程（SDE）的相關性質
  • 二：方程的數值解格式
  • 三：實驗模擬
  • 四：總結

前言

隨機微分方程的發展分支可以說到現在非常廣泛了，從19世紀的布朗運動發現到上世紀40年代的伊藤積分再到正向伊藤隨機微分方程解的存在唯一性定理，后又跨越到上世紀90年代的倒向隨機微分方程解的存在唯一性定理發現，最后到近現在的帶跳、帶馬爾科夫、帶延遲等一系列正倒向隨機微方程解的存在唯一性定理等發現，可以說體系越發廣泛，這也是順應時代發展的必然趨勢。

本篇文章將數值模擬經典下的正向隨機微分方程數值格式的模擬，由于不是存粹數學研究，有些假設和定理和推導會簡單概括，重點是計算機編程模擬，也算是學習隨機微分方程數值法的入門介紹。

一：隨機微分方程（SDE）的相關性質

• 簡介

隨機過程 $X_t$ 被稱作一個擴散過程，如果它滿足下面的隨機微分方程：

$X_t=X_0+{\int }_{t_0}^{t}a(s,X_s)ds+{\int }_{t_0}^{t}b(s,X_s)d{W}_s,t\in [t_0,T],$ 其中 $a_s=a(s,X_s)$ 和 $b_s=b(s,X_s)$ 為 $[0,T]\times R^m$ 上的可測函數， $W_t $為概率空間 $(\Omega ,F,P)$ 上的 $d?$維標準維納過程。但為了此隨機模型在實際中得以應用，保證上述擴散過程得以存在，我們將對此擴散過程給與一定條件，此$SDE$的解才會存在且唯一！

• SDE解的存在唯一性定理

令 $T>0$ ，且對于擴散過程和漂移過程分別有

$a(.,.):[0,T]\times R^m\to R^m,b(.,.):[0,T]\times R^m\to R^{m\times d}$ 為可測函數，

滿足 $\left|a(s,x)+b(s,x)\right|\le C(1+|x|),x\in R^m,s\in [t_0,T]$,

$\left|a(t,x)?a(t,y)\right|+\left|b(t,x)?b(t,y)\right|\le L|x?y|,x,y\in R^m,s\in [t_0,T]$

其中 $C,L$ 為有界常數，令 $\Phi$ 為獨立于由 $W_s,s\ge 0$ 產生的 $\sigma$ 代數 $F$ ，使得

$E\left[|Z^2|\right]<\infty$。

則上述隨機微分方程存在唯一的連續解，其中 $X_0\in \Phi$ ,并且我們還知道

$X_t$ 適應于 $\Phi$ 和 $W_s,s\le t$ 產生的域流 $F_t^{\Phi }且E\left[{\int }_0^T|X_t{|}^{2}dt\right]<\infty$。

• 結合伊藤公式的一般形式的隨機算子表示

對于 m 維擴散過程和 d 維的維納過程，且

$\phi (t,X)\in C^{2,1}\left([0,T]\times R^m\right)$ ，定義算子

$L^0\phi ,L^j\phi 從[0,T]\times R^m\to R^m$ :

$L^0\phi =\frac{?\phi }{?t}+{\sum }_{k=1}^m\frac{?\phi }{?x_k}{a}_k+\frac{1}{2}{\sum }_{q=1}^d{\sum }_{l,k=1}^q\frac{{?}^2\phi }{\phi x_k\phi x_l}{b}_{l,q}{b}_{k,q},L^j\phi$

$={\sum }_{k=1}^m\frac{?\phi }{\phi x_k}{b}_{k,j}$。

二：方程的數值解格式

• 伊藤展開

使用上述隨機算子后，對于含任$d$ 為維納過程的 $\phi (t,X)\in C^{2,1}\left([0,T]\times R^m\right)$ ，并且在時間 $t_n$ 處展開，由伊藤公式我們有如下積分形式為：

$\phi (t_{n+1},X_{t_{n+1}})=\phi (t_n,X_{t_n})+{\int }_{t_n}^{t_{n+1}}{L}^0\phi (s,X_s)ds+$

${\sum }_{j=1}^d{\int }_{t_n}^{t_{n+1}}{L}^j\phi (s,X_s)d{W}_s^j$,

為了簡化模型，我們對上述的只含一維維納過程的 \varphi(.,.) 繼續使用伊藤公式可得：

$\phi (t_{n+1},X_{t_{n+1}})=\phi (t_n,X_{t_n})+$

${\int }_{t_n}^{t_{n+1}}\left[L^0\phi (s,X_s)+{\int }_{t_n}^s{L}^0{L}^0\phi (\tau ,X_{\tau })+{\int }_{t_n}^s{L}^1{L}^0\phi (\tau ,X_{\tau })d{W}_{\tau }\right]ds+$

${\int }_{t_n}^{t_{n+1}}\left[L^1\phi (s,X_s)+{\int }_{t_n}^s{L}^0{L}^1\phi (\tau ,X_{\tau })d\tau +{\int }_{t_n}^s{L}^1{L}^1\phi (\tau ,X_{\tau })d{W}_{\tau }\right]d{W}_s$.

PS:只要上述 $\phi (.,.)$ 函數足夠光滑，我們可以繼續展開，可以得到更多的復雜高階數值格式！

• 收斂階的定義

給定時間分割 $0=t_0<t_1<...<t_{N?1}<t_N=T$ ，假設 $\Delta t=t_{n+1}?t_n$ ,

$\Delta W_n=W_{t_{n+1}}?W_{t_n}(n=0,1,2,...,N)$,

且存在光滑函數 $G\in C_p^{2(\beta +1)}$ 和正常數 $C$ 滿足不等式

$\left|E[G(X_t)?G(X^N)]\right|\le C\left(1+E[|X_0{|}^b]\right)(\Delta t{)}^{\beta }$ ,

其中 $\beta ,b\in N$ ,則 $X^N$ 以階數 $\beta$ 弱收斂于 $X_T$. 當然有弱收斂就有強收斂，當我們不以整個軌跡誤差取期望再絕對值，取而代之的是以每個點的絕對誤差取期望后，就是我們的強收斂誤差要求，高階強收斂格式更加復雜，條件更加苛刻，階數可以是0.5的倍數關系。

• 經典數值格式介紹

設初始條件 $X_0$ 已知，對于 $0\le n\le N?1$ ，對于 $m$ 維的伊藤過程的第 $k$ 個分量，我們有如下格式：
1：歐拉格式（弱1.0階，強0.5階）：

$X_k^{n+1}=X_k^n+a_k\Delta t+b_k\Delta W_n$

2：米爾斯坦格式（強弱都是1.0階）：

$X_k^{n+1}=X_k^n+a_k\Delta t+b_k\Delta W_n+\frac{1}{2}{L}^1{b}_k\left((\Delta W_n{)}^2?\Delta t\right)$

3：弱二階格式：

$X_k^{n+1}=X_k^n+a_k\Delta t+b_k\Delta W_n+\frac{1}{2}{L}^1{b}_k\left((\Delta W_n{)}^2?\Delta t\right)+\left(b_k+\frac{1}{2}\Delta t(L^0{b}_k+L^1{a}_k)\Delta W_n\right)+\frac{1}{2}{L}^0{a}_k(\Delta t{)}^2$

三：實驗模擬

我們以幾何布朗運動 $X_t=a{X}_{t}dt+b{X}_{t}d{W}_t$ ，只要簡單的對 $ln(X_t)$ 使用伊藤公式，就可得此模型真解為$X_t=X_{0}exp[(a?\frac{1}{2}{b}^2)t+b{W}_t]$ ，所以在實驗中能清楚地比較真解和數值解的誤差結果，正向隨機微分方程經典模型還有$Ornstein?Uhlenbeck$過程 $X_t=a(\mu ?X_t)dt+bd{W}_t$ 等，也是存在真解的，但實驗內容是跟幾何布朗運動完全一致的。

import matplotlib.pyplot as plt import numpy as np plt.rcParams['font.sans-serif']=['SimHei']##中文亂碼問題！ plt.rcParams['axes.unicode_minus']=False#橫坐標負號顯示問題！def schemes(a,b,names):logtimels = []logerrorls = []CRls = []msg = '''1：exp函數2：sin函數3：原函數的3次方4：原函數'''print(msg)##定義初始化內容choose = int(input('請選擇g(x)函數：'))echoes = 100time_ls = [8,16,32,64,128]error_ls = []value_matx = np.zeros((echoes,len(time_ls)))Euler = 'y + a * dt * y + b * y * dwt_temp'Milstein = Euler + ' + 0.5 * b * b * y * (dwt_temp * dwt_temp - dt)'weak_2_order = Milstein + ' + 0.5 * a * a * y * dt * dt + 0.5 * dt * (b * a * y + a * b * y) * dwt_temp'schemes_ls = [Euler,Milstein,weak_2_order]for s,n in zip(schemes_ls,names):for e in range(echoes):for i in range(len(time_ls)):x0 = 1##SDE真解初始值y = 1##數值格式初始值T = 1N = time_ls[i]dt = 1 / Nyt = []dwt = np.random.normal(0,1,N) * np.sqrt(dt)Wt = np.cumsum(dwt)###真解WT = Wt[N - 1]vtemp = x0 * np.exp((a - 0.5 * b * b) * T + b * WT)xT = 0if choose == 1:xT = np.exp(vtemp)if choose == 2:xT = np.sin(vtemp)if choose == 3:xT = np.power(vtemp,3)if choose == 4:xT = vtemp##數值格式解for j in range(N):dwt_temp = dwt[j]y = eval(s)if choose == 1:yt.append(np.exp(y))if choose == 2:yt.append(np.sin(y))if choose == 3:yt.append(np.power(y,3))if choose == 4:yt.append(y)error_temp = yt[N - 1] - xT ##真解與數值解末端值作差error_ls.append(error_temp)value_matx[e,:] = error_lserror_ls = []error_fin = np.abs(sum(value_matx) / echoes)##弱收斂最終求期望log_time = np.log(1 / np.array(time_ls))log_error_fin = np.log(error_fin)CR = np.polyfit(log_time,log_error_fin,1)##一次多項式擬合print('\033[1;36m{0:*^80}\033[0m'.format('計算結果'))print('格式%s的相對誤差為：%s'%(n,np.round(error_fin,6)))print('格式%s的弱收斂階為：%s'%(n,round(CR[0],6)))logerrorls.append(np.round(log_error_fin,6))logtimels.append(log_time)CRls.append(round(CR[0],6))return logtimels,logerrorls,CRlsmu = 1##漂移系數 sigm = 0.05##擴散系數 scheme_names = ['Euler', 'Milstein', 'weak_2_order'] R = schemes(mu,sigm,scheme_names)

圖1

圖2

圖3 marker_ls = ['-*','-o','-^'] def make_figure():plt.figure(figsize=(8, 6))for i,j,k,l,m in zip(R[0],R[1],R[2],scheme_names,marker_ls):plt.plot(i,j,m,label=l)plt.plot(np.array([-4,-3]),np.array([-4,-3]) + 1,'--',label='slope=1')plt.plot(np.array([-4, -3]), 2 * np.array([-4, -3]) - 1.5,'--',label='slope=2')plt.xlabel('stepsize logN')plt.ylabel('log(|E(Y_T-X^N)|)')plt.legend()plt.show() make_figure()

圖4

四：總結

經典的正向隨機微分方程數值格式和收斂階的模擬大致如上，但最開始也提到過，由正向衍生出來的各種方程在如今是非常多的。與正向相對性的倒向隨機微分方程也是一個大類研究方向，今后如果有機會，作者還會發布些關于倒向隨機微分方程的數值模擬，感興趣的朋友希望點贊關注加收藏啦！

總結

以上是生活随笔為你收集整理的正向随机微分方程的经典数值格式模拟的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。