多層次蒙特卡洛(multilevel monte carlo) · 自適應多層次蒙特卡洛(adaptive MLMC) · 障礙期權

04-29

以賭場著名的摩納哥的一個城區，蒙特卡洛，連接了數學裡互相瞧不起的numeric和 stochastic方向。

蒙特卡洛方法在金融衍生品定價的應用，不言而喻，以致人們淡化了數學裡追求的解析解。

這裡，作為mc的愛好者，苟膽也說說蒙特卡洛的一種優化，複習下多層次蒙特卡洛（multilevel monte carlo），以及其中的一個擴展adaptive的mlmc。前者，由牛津大牛 giles得以發展弘揚，後者發表於之前斯德哥爾摩大學的一個小組（H?kon Hoel， Erik von Schwerin, Anders Szepessy, Raúl Tempone）。自適應多層次可以把歐式障礙期權的計算成本從普通的

$mathcal{O}(TOL^{-4}), 降低為mathcal{O}(TOL^{-2})，$

且不受多障礙、多維度的影響。

在金融裡面，會對標的價格（比如股票價格，spot rate）假設到某些stochastic process，所以有很多模型的名字，經典的black scholes, heston, HJM, Vasicek, Hull-white模型，而蒙特卡洛方法，就是模擬這些標的物（路徑path）價格的，目的是求期末payoff的期望值

$mathbb{E}[g(X_T)]$ ，

其中g是payoff的functional函數， $X_T$ 是標的物在期末T的價格。

比如，假設股票價格服從 $ilde{d}$ 維幾何布朗運動 $egin{equation*} X^{(i)}_t = X^{(i)}_0 + int_0^t mu^{(i)} X^{(i)}_sds+ sum_{k=1}^{{ ilde{d}}}int_0^t sigma_{k}^{(i)}X^{(i)}_sdW^{(k)}_s end{equation*}$ 。我們看到，要模擬一條path，自然就必須考慮到其時間的離散化問題。在小區間里 $[t_n,t_{n+1})in [0,t]$ 里，我們有 $egin{equation*} X^{(i)}_{t_{n+1}} = X^{(i)}_{t_n} + int_{t_n}^{t_{n+1}} mu^{(i)} X^{(i)}_sds+ sum_{k=1}^{{ ilde{d}}}int_{t_n}^{t_{n+1}} sigma_{k}^{(i)}X^{(i)}_sdW^{(k)}_s end{equation*}$ ，於是我們有最簡單的Euler-Maruyama離散方法，因為

$egin{align} X^{(i)}_{t_{m+1}}&approx X^{(i)}_{t_m}+ mu_i(t_m,X^{(i)}_{t_m})int_0^t ds + sum_{k=1}^{{ ilde{d}}}sigma_{i,k}(t_m,X^{(i)}_{t_m})int_0^t dW^{(k)}_{t_m} otag\[0.2cm] &= X^{(i)}_{t_m} + mu_i(t_m, X^{(i)}_{t_m})(t_{m+1}-t_m)+ sum_{k=1}^{{ ilde{d}}} sigma_{i,k}(t_m, X^{(i)}_{t_m}) (W^{(k)}_{t_{m+1}} - W^{(k)}_{t_m}) otag\[0.2cm] &:=X^{(i)}_{t_m} + mu_i(t_m, X^{(i)}_{t_m})Delta t_m+ sum_{k=1}^{{ ilde{d}}} sigma_{i,k}(t_m, X^{(i)}_{t_m}) Delta {W^{(k)}_{t_m}}. label{euler} end{align}$

Euler-Maruyama離散方法對上面的SDE有弱收斂性convergence order 1。

假如我們利用ito公式繼續對上面展開，會有milstein schema，以及general的taylor展開式離散方式，感興趣的同學可以看看大牛P. Kloeden, E. Platen的Numerical Solution of Stochastic Differential Equations。（上過Kloeden的2門課，目前他退休後在華中科大撈錢）。

用Euler-Maruyama離散方式的蒙特卡洛方法稱之為euler monte carlo。我們現在要用euler mc模擬期末的payoff，其實就是模擬N條path，對每個path的期末的payoff值求均值：

$egin{align}label{mceuler} mathbb{E}[g(X_T)]&approx sum_{i=1}^N frac{gleft(hat{X}_T(Delta t, w_i) ight)}{N}\ &:=mathcal{A}_{MC} otag, end{align}$

對任何的數值分析，判斷好壞的標準當然是誤差與收斂級數。這裡，誤差有

$egin{align}label{errorEntw} mathcal{E}&=mathbb{E}[g(X_T)]-mathcal{A}_{MC} otag\ &=mathbb{E}[g(X_T)]-mathbb{E}[g(hat X_T)]+mathbb{E}[g(hat X_T)]-mathcal{A}_{MC} otag\ &= left(mathbb{E}left[g(X_T)-g(hat X_T) ight] ight) +left( mathbb{E}[g(hat X_T)]-mathcal{A}_{MC} ight)\ &:=mathcal{E} _T + mathcal{E}_S otag, end{align}$

他可以拆為因為離散而造成的離散誤差 $mathcal{E} _T$ ，以及N個sample求均值的統計誤差 $mathcal{E}_S$ 。這樣人們就可以通過各自的tolerance $TOL_T,$ 和 $TOL_S$ 來控制它。其中

$egin{align} mathcal{E} _T &=mathbb{E}left[g(X_T)-g(hat X_T) ight] otag\ &=mathcal{O}(maxDelta t)=mathcal{O}(Delta t)label{zeiterrorordnung} end{align}$

和由於中心極限定理得到的

$egin{align*} frac{sum_{i=1}^N g(hat{X}_T(Delta t, w_i))-Nmathbb{E}[g(hat{X}_T)]}{sqrt{N Var(g(hat{X}_T))}}&=frac{sqrt{N}mathcal{E}_S}{sqrt{Var(g(hat{X}_T))}}overset{d}{sim}mathcal{N}(0,1), end{align*}$

$egin{align}label{Nordnung} mathcal{E}_Soverset{d}{sim}mathcal{N}left(0, frac{Var(g(hat{X}_T))}{{N}} ight). end{align}$ 所以

$egin{align*} mathcal{E}_S in left(-Phi^{-1}(frac{1+gamma}{2})sqrt{frac{Var(g(hat{X}_T))}{{N}}}, Phi^{-1}(frac{1+gamma}{2})sqrt{frac{Var(g(hat{X}_T))}{{N}}} ight) end{align*}$ ,所以

$egin{align}label{s_error} |mathcal{E}_S|lesssim C_C sqrt{frac{Var(g(hat{X}_T))}{{N}}} end{align} , C_C是mathrm{ Quantil}$ ，所以

統計誤差的收斂order只有0.5。因為我們定義蒙特卡洛的成本(complexity)是每個path的離散點M與模擬多少條path N的乘積，所以euler MC 對lipschit 連續的payoff的成本是 $egin{align}label{mckomplexit?t} mathcal{C}=mathcal{O}(Mcdot N)=mathcal{O}left(frac{1}{TOL}cdotfrac{1}{TOL^2} ight)=mathcal{O}(frac{1}{TOL^3}). end{align}$

多層次蒙特卡洛（multilevel monte carlo）對lipschitz連續的payoff g有complexity order $mathcal{C}=mathcal{O}left(TOL_{RMSE}^{-2}(log TOL_{RMSE})^2 ight)$ ，相對於普通的mc $mathcal{C}=mathcal{O}(TOL_{RMSE}^{-3})$ 。

$egin{align*} mathbb{E}[g(X_T)]&=mathbb{E}left[g_{L}left(hat{X}_T ight) ight]\[0.2cm] &=mathbb{E}left[g_0left(hat{X}^{}_T ight) ight]+sum_{l=1}^Lmathbb{E}left[g_{l}left(hat{X}^{}_T ight)-g_{l-1}left(hat{X}^{}_T ight) ight]\[0.2cm] &approx sum_{i=1}^{N_0} frac{g_0left(hat{X}_T^{}(W_i) ight)}{N_0}+sum_{l=1}^L left( sum_{i=1}^{N_l} frac{g_lleft(hat{X}_T^{}(W_i) ight)-g_{l-1}left(hat{X}_T^{}(W_i) ight)}{N_l} ight)\[0.2cm] &:= sum_{l=0}^L sum_{i=1}^{N_l} frac{Delta^{}g_lleft(hat{X}(W_i) ight)}{N_l}\[0.2cm] &=sum_{l=0}^L mathcal{A}(Delta^{}g_l,N_l)\[0.2cm] &=mathcal{A}_{ML}, end{align*}$

其構造的核心是構造一個隨層數步數幾何遞增的grid（ $Delta t_{l}equivfrac{T}{ ilde{M}^l}， ilde{M}:=2$ ），通過拉格朗日方法再優化每層模擬次數 $N_l$ 和最高層數L，從而減少模擬的成本，而達到所設定精度TOL。我們對第i個標的物，第j層的模擬 $hat{X}_l^{}(T,W_i)$ 通關對 $j=0,1,..., ilde{M}^l$ 的遞歸iterative：

$egin{align*} hat{X}_{l}^{}(t_{j+1},W_{i,j+1})=hat{X}_{l}^{}(t_{j},W_{ij})+aleft(t_j,hat{X}_{l}^{}(t_{j},W_{ij}) ight)cdot frac{T}{ ilde{M}^l}+bleft(t_j,hat{X}_{l}^{}(t_{j},W_{ij}) ight)sqrt{frac{T}{ ilde{M}^l}} z_{i,j} end{align*}$ ，

以及對 $hat{X}_{l-1}^{}(T,W_i)$ 的遞歸有

$egin{align*} hat{X}_{l-1}^{}left(t_{(k+1) ilde{M}},W_{i,(k+1) ilde{M}} ight)&=hat{X}_{l-1}^{}left(t_{k ilde{M}},W_{i,k ilde{M}} ight)+aleft(t_{k ilde{M}},hat{X}_{l-1}^{}left(t_{k ilde{M}},W_{i,k ilde{M}} ight) ight)cdot frac{T}{ ilde{M}^{l-1}}\ &+bleft(t_{k ilde{M}},hat{X}_{l-1}^{}left(t_{k ilde{M}},W_{i,k ilde{M}} ight) ight)sqrt{frac{T}{ ilde{M}^l}} sum_{j=(k-1) ilde{M}+1}^{k ilde{M}} z_{i,j}, end{align*}$

因為l-1層的增量有：

$egin{align*} Delta W^{(l-1)}_{i,k ilde{M}}&=W^{(l-1)}_{i,k ilde{M}}-W^{(l-1)}_{i,(k-1) ilde{M}}\&=W^{(l)}_{i,k ilde{M}}-W^{(l)}_{i,(k-1) ilde{M}}\ &=W^{(l)}_{i,k ilde{M}}-W^{(l)}_{i,k ilde{M}-1}+W^{(l)}_{i,k ilde{M}-1}-...+W^{(l)}_{i,(k-1) ilde{M}+1}-W^{(l)}_{i,(k-1) ilde{M}}\ &=sqrt{frac{T}{ ilde{M}^l}}z_{i,k ilde{M}}+sqrt{frac{T}{ ilde{M}^l}}z_{i,k ilde{M}-1}+...+sqrt{frac{T}{ ilde{M}^l}}z_{i,(k-1) ilde{M}+1}\ &=sqrt{frac{T}{ ilde{M}^l}} sum_{j=(k-1) ilde{M}+1}^{k ilde{M}} z_{i,j}. end{align*}$

其模擬的圖如下

記L為最高層， $ilde{M}^l$ 為第l層的步數， $N_l$ 為第N層需要模擬的次數。根據上面的grid的構造，有

$egin{align*} M_{l}:= ilde{M}^l=mathcal{O}(TOL^{-1})， end{align*}$

所以有 $egin{align}label{MM_ordnung} ilde{M}^l+ ilde{M}^{l-1}=mathcal{O}(TOL^{-1})+frac{1}{ ilde{M}}mathcal{O}(TOL^{-1})=mathcal{O}(TOL^{-1})， end{align}$ 因此在最高處L，存在一個常數c，使得

$egin{align}label{L_ordnung} &M_{L}:= ilde{M}^L le c^* TOL^{-1} otag\ &Rightarrow Lle c frac{log TOL^{-1}}{log ilde{M}}. end{align}$

通過拉格朗日求最優解

$egin{align*} egin{cases} min_{N_l}Var(mathcal{A}_{ML})=sum_{l=0}^L frac{mathcal{V}_l}{N_l} quadquadquadmathrm{(目標方程)}\ mathcal{C}=N_0+sum_{l=1}^L N_l( ilde{M}^l+ ilde{M}^{l-1}).quadquadmathrm{(條件約束)} end{cases} end{align*}$

其中第l層模擬方差

$egin{align*} Varleft(mathcal{A}(Delta ^{}g_l) ight)&=N_l^{-2} sum_{i=1}^{N_l}Varleft(Delta ^{}g_l(X_T(W_i)) ight)\ &=frac{mathcal{V}_l}{N_l} ,mathrm{規定}mathcal{V}_l:=Varleft(Delta^{}g_l(X_T(W_i)) ight)， end{align*}$

所以所有的MLMC方差是

$egin{align*} Var(mathcal{A}_{ML})&=Varleft(sum_{l=0}^L mathcal{A}(Delta^{}g_l) ight)\ &=sum_{l=0}^L frac{mathcal{V}_l}{N_l}。 end{align*}$

實際上，給定 $TOL_{RMSE}$ ， $MSE(mathcal{A}_{ML})=Var(mathcal{A}_{ML})+left(E[g-mathcal{A}_{ML}] ight)^2le TOL_{RMSE}^2$ ,各自平攤一半誤差，即 $egin{align} Var(mathcal{A}_{ML})le frac{TOL_{RMSE}^2}{2}label{Var_2}\ mathbb{E}[g-mathcal{A}_{ML}]^2le frac{TOL_{RMSE}^2}{2}, end{align}$

求得當滿足 $egin{align}label{Stopkriterium} mathcal{A}(Delta g_L, N_L)approxmathbb{E}(g_L-g_{L-1})le frac{1}{sqrt{2}}( ilde{M}-1)TOL_{RMSE} end{align}$ 時，就可以滿足設定的精度。

最後那個拉格朗日有最優解：

$egin{align}label{N_l_end} N_l=leftlceil 2 TOL_{RMSE}^{-2}sqrt{mathcal{V}_lDelta t_l}left(sum_{i=0}^Lsqrt{mathcal{V}_i/Delta t_i} ight) ight ceil. end{align}$

實際演算法就是從L=0層開始，模擬10000條path，求得empirically $N_l$ ,假若求得的 $N_l$ 大於10000，則在該層繼續模擬 $N_l - 10000$ 條path，L=L+1，直到滿足 $egin{align}label{Stopkriterium} mathcal{A}(Delta g_L, N_L)approxmathbb{E}(g_L-g_{L-1})le frac{1}{sqrt{2}}( ilde{M}-1)TOL_{RMSE}。 end{align}$

自適應多層次蒙特卡洛

自適應演算法是種後驗式的演算法，最重要的是獲得後驗誤差密度函數 $ho$ ，一般有如下形式

$egin{align*} mathbb{E}[g(X(T))-g(hat{X}(T))]lesssimmathbb{E}left[sum_{m=1}^M ho_m(Delta t_m)^2 ight] end{align*}$ ，其中 $ho_m$ 是個比較複雜的backwards function。根據後驗誤差判斷是否對步長細分。

$egin{align*} { ho}left(t_m,hat{X}(t_m) ight):&=frac{1}{2}left(frac{partial}{partial_t} mu_k+a_jpartial_j mu_k + d_{ij}partial_{ij}^2 mu_k ight)varphi_k(t_{m+1})\[0.2cm] &+frac{1}{2}left(frac{partial}{partial_t} d_{kq}+ a_jpartial_j d_{kq} + d_{ij}partial_{ij}^2 d_{kq} + 2 d_{jq}partial_{j} a_k ight)varphi_{kq}(t_{m+1})\[0.2cm] &+ (d_{jr}partial_{j}d_{kq})(hat{X}(t_m))varphi_{kqr}(t_{m+1})), end{align*}$

其中 $d_{ij}:=frac{1}{2}sigma_i^{(l)} sigma_j^{(l)}quad i,j=1,2,...,d,$ 小標的表達方式是愛因斯坦格式，意思是下表重複時對其去和，比如

$mu_jpartial_j mu_k+d_{ij}partial_{ij}^2 mu_k :=sum_{k=1}^d left( sum_{j=1}^d mu_j partial_j mu_k+sum_{l=1}^{ ilde{d}} sum_{i=1}^dsum_{j=1}^dfrac{1}{2}sigma_i^{(l)}sigma_j^{(l)}partial_{x_i x_j}^2 mu_k ight)$ 。

密度函數裡面其他的式子就不說明，因為比較長（看文獻1）.

對於 $(X)_t$ 服從幾何布朗過程的歐式call option的payoff的話，很簡單：

$egin{align*} { ho}(t_m,hat{X}(t_m))=frac{1}{2}mu^2X(t_m)varphi(t_{m+1}) end{align*}$ ，以及backward 方程

$egin{align*} &varphi(t_m)=(1+mu Delta t_m+ sigmaDelta W_m)varphi(t_{m+1})quadquadmathrm{for } t_m<T,\ &varphi(T)= egin{cases} 1quadquadmathrm{for }hat{X}_T>K\ 0quadquadmathrm{others}. end{cases} end{align*}$

在對path的步數M運算

$egin{align} M&=int_0^T frac{1}{ ilde{Delta t}( au)}d au otag\ &=int_0^T frac{sqrt{|{ ho}( au)|} int _0^T sqrt{| ho(s)|}ds}{TOL_T}d au otag\ &=frac{1}{TOL_T}left(int_0^Tsqrt{|{ ho}( au)|}d au ight)^2 otag\ &=frac{1}{TOL_T}Vert hoVert_{L^{frac{1}{2}}(0,T)}label{Adaptiv_M_DichtNorm}， end{align}$

其中 $ilde{Delta t}( au)$ 是細化了的步長，得到誤差密度與離散誤差的關係 $egin{align}label{Adaptiv_Indikator_Konstante} |{ ho}_m|( ilde{Delta t}_m)^2=frac{TOL_T}{M}， end{align}$

從而就有了刻畫步長誤差的 $egin{align}label{verfeinerung_delta_t} |{ ho}_m|( ilde{Delta t}_m)^2=r(m)lefrac{TOL_T}{M} quadforall m=0,1,2,...,M-1 end{align}$

如果， $egin{align}label{C_R_Kriterium} r_m> C_Rfrac{TOL_T(C_S)}{mathbb{E}[M]}quadmathrm{(criterion for refinement)} end{align}$ ，則對該步長細分的，如果所有的每個小區間（步長）的誤差都比平均每個小區間應有的誤差小

$egin{align}label{C_S_Stopp} max_{1le mle M}r_mle C_Sfrac{TOL_T(C_S)}{mathbb{E}[M]}quad mathrm{(criterion for stop refinement)}, end{align}$ 就可以保證整個的離散誤差在 $TOL_T$ 內。其中 $C_R:=2, C_S:=3$ 兩個控制細分速度的參數（不是2跟3也行）。

用 $TOL_T$ 控制步長，接下來就用 $TOL_S$ 控制要模擬多少條path，一致達到所要的精度，通過方差最小化，有

$egin{align}label{M_Update_1} N=leftlceil frac{C_C^2Var(g(hat{X}(T)))}{TOL_S^2} ight ceil, end{align}$ 其中 $C_C:=1.65$ 可以是95%的分位數。

接下來就要把這adaptive演算法打包到多層次蒙特卡洛上去。

（演算法此處省略，如果感興趣可以看文獻1，如果對文獻1的演算法還是不滿意，私信阿廣，呵呵）

障礙期權定價

假設是在black scholes假設下，標的物價格服從幾何布朗過程 $egin{align*} dX_t=rX_t dt+sigma X_t dW_t end{align*}$ ，單一敲出Call障礙期權的payoff是

$egin{align*}egin{cases} max(X_T-K, 0)&mathrm{if } X_t<B, forall tin [0,T]\ 0&mathrm{others,} end{cases} end{align*}$

在風險中性測度的話，有 $egin{align} label{Auszahlung_Barrier_1} g(X_T)&=exp(-rT)max(X_T-K, 0){1}_{{ augeq T}} end{align}$ ，其中 $au$ 是首次突破障礙的時間 $egin{align*} au:=inf{t:(X_t,t) otin D imes [0,T]}，D=(-infty,B], Binmathbb{R}. end{align*}$

在障礙期權里，標準euler mc的弱收斂性order只有0.5，（正常option的為1） $egin{align}label{Barriere_MCE_Ordnung} mathbb{E}[g(X_{ au}, au)-g(hat{X}_{hat{ au}},hat{ au})]=mathcal{O}(frac{1}{sqrt{M}}) end{align}$ ，所以標準Euler mc的complexity在這裡是 $mathcal{O}left(frac{1}{TOL^4} ight),$ 但是有一個藉助已知兩點，在這兩點的布朗橋突破已知障礙的解析概率 $egin{align}label{Austrittswahrscheinlichkeit} P_m&=mathbb{P}left(max_{tin[t_m,t_{m+1})} X_t>B| X({t_m})=hat{X}_{m}, X(t_{m+1})=hat{X}_{m+1} ight) otag\[0.2cm] &=expleft(-2frac{(B-hat{X}_m)^+(B-hat{X}_{m+1})^+}{sigma^2 hat{X}_m^2 Delta t_m} ight), 其中Delta t_m=t_{m+1}-t_m end{align}$