[隨機分析] Euler-Maruyama 法求 隨機微分方程 數值解 (利用 MATLAB)

在此我們介紹 Euler-Maruyama Method 求解 隨機微分方程   (Stochastic Differential Equation, SDE)，現在 考慮  SDE 的積分型式 可寫為：
$$X(t) = X_0 + \int_0^t f(X(s))ds + \int_0^t g(X(s))dW(s), \;\; 0 \le t \le T$$ 其中 $f, g$ 為 純量函數 且 初始值 $X_0$ 為隨機變數；另外 第二項積分為 Ito integral。上式若有解則其解 $X(t)$ 對任意 $t$ 皆為隨機變數。

一般而言，上式可改寫為較為簡潔的 隨機微分方程的微分形式：
$$dX(t) = f(X(t))dt + g(X(t)) dW(t),\;\; X(0)=X_0, \; 0 \le t \le T \ \ \ \ \ (*)$$
Comments: 
1. 讀者需注意我們不可寫 $dW(t)/dt$ 因為 Browian motion 為 處處連續但處處不可微 (with probability 1)
2. 若 $g = 0$ 且 $X_0$ 為常數，則 SDE 退化成一般的 ODE；亦即
$$dX(t) = f(X(t))dt, \;\; X(0)=X_0, \; 0 \le t \le T$$ (此時即可用 Euler method 求數值解。)
3. 關於 SDE 何時有解，讀者可參考BLOG 相關系列文章 ：
[隨機分析] Uniqueness and Existence theorem for S.D.E. (1)- Uniqueness
[隨機分析] Uniqueness and Existence theorem for S.D.E. (2) - Picard Iteration for SDE
[隨機分析] Uniqueness and Existence theorem for S.D.E. (3) - An intermediate result (Upper bound) of Picard Iteration
[隨機分析] Uniqueness and Existence theorem for S.D.E. (4) - The Existence of Solution



那麼現在回歸主題，在此我們的目的是：
對 SDE 在關心的區間 $t \in [0,T]$ 之間進行數值求解 (透過 MATLAB )。

第一步首先將此區間 $[0,T]$ 離散化，亦即 對某些 $N$ 而言，定義 $\Delta t := T/N$ 且 對 $j=1,...,N$ ，定義 $\tau_j := j \Delta t$。 並且稱 數值近似解 $X_j :=X(\tau_j)$。

利用 Euler-Maruyama (EM) 法，我們可將前述 SDE 在離散化區間內 改寫如下 $${X_j} = {X_{j - 1}} + f({X_{j - 1}})\Delta t + g({X_{j - 1}})(W({\tau _j}) - W({\tau _{j - 1}})),\;\;j = 1,2,...,N$$
Comment: 注意到上式若寫成積分形式即為
$$X({\tau _j}) = X({\tau _{j - 1}}) + \int_{{\tau _{j - 1}}}^{{\tau _j}} f (X(s))ds + \int_{{\tau _{j - 1}}}^{{\tau _j}} g (X(s))dW(s)$$

Example
利用 Euler-Maruyama 法求解下列線性 SDE
$$dX(t) = \mu X(t) dt + \sigma X(t) dW(t),\;\; X(0) = X_0$$ 其中 $\mu, \sigma$ 為 常數 (亦即在此例我們選 $f(X(t)) = \mu X(t)$ 且 $g(X(t)) = \sigma X(t)$)。此例為 Geometric Brownian Motion (GBM) 模型，且此 SDE 有解析解如下：
$$X\left( t \right) = {X_0}{e^{\left( {\mu  - \frac{{{\sigma ^2}}}{2}} \right)t + \sigma W\left( t \right)}}$$ 以下我們將利用 MATLAB 實現 上述 解析解 與透過 EM 法 求數值解並比較其差異。


現在我們利用 MATLAB 實現 Euler-Maruyama 法：

首先計算 discretized Brownian path over $[0,1]$
定義 stepsize $\Delta t = R \delta t$ (一般泛取 $R=1$)且 由於 EM 法 $${X_j} = {X_{j - 1}} + f({X_{j - 1}})\Delta t + g({X_{j - 1}})(W({\tau _j}) - W({\tau _{j - 1}})),\;\;j = 1,2,...,L$$ 上式還需要計算 $W(\tau_j) - W(\tau_{j-1})$，故我們計算：
$$\begin{array}{l} W({\tau _j}) - W({\tau _{j - 1}}) = W\left( {j\Delta t} \right) - W\left( {\left( {j - 1} \right)\Delta t} \right)\\ \begin{array}{*{20}{c}} {}&{}&{}&{}&{}&{}&{}&{} \end{array} = W\left( {jR\delta t} \right) - W\left( {\left( {j - 1} \right)R\delta t} \right)\\ \begin{array}{*{20}{c}} {}&{}&{}&{}&{}&{}&{}&{} \end{array} = \sum\limits_{k = jR - R + 1}^{jR} {d{W_k}} \end{array}$$ 上式計算出現在下列程式碼 (第14行)


我們將 MATLAB 程式碼附上如下： (考慮 $\mu=2$, $\sigma =1$ 且 $X_0=1$。 $R=1$)
(line 10 : line 14: 產生 standard Brownian Motion)
(line 16: 解析解 for GBM SDE)
(line 19: line 25: EM-method)

執行結果如下：

上圖中 紅線為 EM 法結果，藍線為 解析解結果。讀者可調整 $R$ 值來檢驗兩者差距。


ref: Desmond J. Higham, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM REVIEW Vol 43, No. 3, pp. 525-546, 2001.