摘要： 首先， 用原始對偶活躍集方法求解期權定價正問題， 將相應的數值解作為監(jiān)督學習的輸出， 然后用訓練好的神經網絡替代期權定價正問題模型. 其次， 結合Bayes推斷與神經網絡進行Metropolis-Hastings采樣， 求解隱含波動率反問題. 該方法減少了采樣過程中正問題計算量龐大的問題， 從而加速了反問題求解過程.

關鍵詞： 隱含波動率； Bayes推斷； 神經網絡； 替代模型

中圖分類號： O241.8""文獻標志碼： A""文章編號： 1671-5489（2024）06-1363-07

Solving "Implied Volatility of American Lookback Options byBayesian Inference and Neural Network

TAO Li1， ZHU Benxi2， QIAN Yiyuan2， XU Jiaqi2

（1. College of International Business， Hainan University， Haikou 570228， China;

2. College of Mathematics， Jilin University， Changchun 130012， China）

Abstract： Firstly， "we used the original dual active set method to solve the forward problem of option pricing， with the corresponding numerical

solutions as the output for supervised learning， and then replaced the "forward problem model of option pricing with a well-trained neural network. Secondly， we combined

Bayesian inference with neural networks for Metropolis-Hastings sampling "to solve the inverse problem of implied volatility. This method reduced the problem of large computational

complexity "of the forward problem during the sampling process， thereby accelerating the solution process for the inverse problem.

Keywords： implied volatility； Bayesian inference； neural network； surrogated model

隱含波動率是期權的一個重要指標， 它反映了市場對于標的價格的變動幅度. 投資者可以利用它來預測標的價格的未來走勢和供求關系. 求解期權的隱含波動率， 可視為解決期權

定價模型的反問題.Achdou等^［1］使用最小二乘公式和下降算法實現了對波動率的求解; Burkovska等^［2］用縮基法代替經典的有限元法， 從而減少了求解的復雜度; Kutner^［3］利用基于二次近似法的數值方法， 計算了美式期權的隱含波動率. 傳統(tǒng)數值方法雖然能保證求解精度， 但當模型為非線性或高維時， 通常需要耗費大量計算時間， 計算效率較低. 近年來， 深度學習的快速發(fā)展提供了解決高維問題的新思路. Han等^［4］采用了一種基于深度學習的方法， 處理一般的高維拋物型偏微分方程; Liu等^［5］建立了一個校準神經網絡框架， 計算了隱含波動率. 上述基于深度學習求解高維偏微分方程的方法， 在保證精度的情況下， 顯著減少了計算所需時間. 本文提出一個利用神經網絡替代正問題的模型， 并結合Bayes推斷的數值方法， 快速求解隱含波動率.

1"正問題描述和數值算法

設V是期權價格， S是當前股票價格， J是股票價格在過去到當前時間的最大值， t是到期日與當前時刻的差值， σ是波動率， r是無風險利率， q

是股息收益率. 美式回望看跌期權價格V=V（S，J，t;σ）滿足如下自由邊界問題^［6］：

LV∶=Vt+12σ2S

22VS2+（r-q）SVS-rV=0，"B（t）lt;S≤Jlt;+∞，"0≤tlt;T，（1）

其中B（t）為未知自由邊界， 收益函數為F（S，J，T）=J-S且J=max0≤ζ≤t Sζ， 波動率函數為

σ（t）=∑di=1σiai（t）.

這里， {ai（t）}di=1是L2（［0，T］）有限維逼近空間的一組已知基底（例如{sin（it）}^d-1i=0， {e^it}^d-1i=0）， {σi}di=1為參數. 則相應的線性互補模型為

LV·（V-G）=0，JX≤S≤Jlt;+∞，"0≤tlt;T，

LV≤0，"V≥G，JX≤S≤Jlt;+∞，"0≤tlt;T，F（S，J，T）=（J-S），JX≤S≤Jlt;+∞，（2）

其中G=J-S是贖回價格， X滿足：

X=-r+q+σ2/2-（-r+q+σ2/2）2+2rσ2-r+q-σ2/2-（-r+q+σ2/2）2+2rσ2.

利用變換x=ln JS， τ=T-t和u（x，τ）=V（S，J，t）S， 線性互補模型（2）可改寫為

（κ（x）uτ-（γκ（x）ux）x+qκ（x）u）·（u-g）=0，0≤ xlt;L，"0lt;τ≤T，

κ（x）uτ-（γκ（x）ux）x+qκ（x）u≥0，"u-g≥0，0≤xlt;L，"0lt;τ≤T，u（x，0）=g（x），0≤x≤L，

ux（0，τ）=0，"u（L，τ）=g（L），0lt;τ≤T，（3）

其中γ=σ22， κ（x）=e^{（q-r-γ）x/γ}， g（x）=ex-1， L=ln 1X.

設H1L（Ω）∶={v∈H1（Ω）; v≥g（x）， v（L）=g（L）}，H1E0（Ω）∶={v∈H1（Ω）; v≥g（x）， v（L）=0}，""Ω=［0，L］.

引理1"設u∈L2（0，T;H2（Ω））， uτ∈L2（0，T;L2（Ω））， 且u（0，τ）gt;g（0）， 0lt;τ≤T. 則u是模型（3）的解當且僅當

u是下列變分不等式問題的解：

（κ（x）uτ，v-u）+（γκ（x）ux，vx-ux）+（qκ（x）u，v-u）≥0，"v∈H1E0（Ω），"0lt;τ≤T，

且u（x，0）=g（x）.

令時間和空間的剖分分別為

Tt： 0=τ0lt;t1lt;…lt;tM=T，"Δt=TM，"tj=jΔt，"j=0，1，…，M，

Th： 0=x0lt;x1lt;…lt;xN=L，"h=LN，"xi=ih，"i=0，1，…，N，

其中M和N為正整數. 設Vh（Ω）H1（Ω）為分片線性函數空間， 取有限元試探函數空間和檢驗函數空間分別為

VhL（Ω）∶={vh∈Vh（Ω）; v（xn）≥g（xn）， n=0，1，…，N-1， vh（L）=g（L）}=span{φi（x）}Ni=1，

VhE0（Ω）∶={vh∈Vh（Ω）; v（xn）≥g（xn）， n=0，1，…，N-1， vh（L）=0}.

則半離散格式為： 求uh（τ）∈VhL（Ω）， 使得uh（x，0）=gI（x）， 且

（κ（x）uhτ，vh-uh）+（γκ（x）uhx，vhx-uhx）+（qκ（x）uh，vh-uh）≥0，"vh∈VhE0（Ω），"0lt;τ≤T，

其中gI（x）為g（x）在VhL（Ω）中的插值.

進一步， 在時間方向采取向后Euler格式， 則全離散格式為

（κ（x）（（1+qΔτ）umh-u^m-1h），vh-umh）+Δτ（γκ（x）umhx，vhx-umhx）≥0，"vh∈VhE0（Ω），"m=1，2，…，M，

其相應的矩陣形式為

（DUm，V-Um）≥（Wm，V-Um），

V≥F1，"m=1，2，…，M，（4）

其中D=γΔτA+（1+qΔτ）B， Wm=BU^m-1-F2，

V=（v0，v1，2，…，vN-1）T， Um=（um0，um1，…，umN-1）T， A=（Ai，j）N×N，

B=（Bi，j）N×N， i，j=0，1，…，N-1， F1=（g（x0），…g（xN-1））T， F2=（0，…，0，g（L）Δτ（qBN-1，N+γAN-1，N））T，

Ai，j=（κ（x）φ′j，φ′i）， Bi，j=（κ（x）φj，φi）， i，j=0，1，…，N-1.

算法1^［7-8］"原始對偶活躍集方法.

輸入： U0，F1，Wm;

輸出： Um，Xml，Xmr.

for m=1∶M， do

k=0， ρgt;0， εgt;0， U^（0）=max{U^m-1，F1}， U^（1）=0， λ^（0）≥0.

while ‖U^（k+1）-U^（k）‖∞lt;ε， do

① 令活躍集為IS^（k+1）， 非活躍集為AS^（k+1）：

IS^（k+1）={i∈N， λ^（k）i+ρ（Gmi-U^（k）i）≤0}，

AS^（k+1）={i∈N， λ^（k）i+ρ（Gmi-U^（k）i）gt;0}，

N={0，1，…，N-1}.

② 求解下列方程， 得到U^（k+1）和λ^（k+1）：

DU^（k+1）-λ^（k+1）=Wm，

U^（k+1）=Gmi， i∈AS^（k+1），

λ^（k+1）i=0， i∈IS^（k+1）.

end while

Um=U^（k+1），

lm=min（IS^（k+1））， rm=max（IS^（k+1）），

xml=xlm， xmr=xrm.

end for

本文采用算法1對式（4）進行數值求解. 為方便， 上述正問題整體求解過程， 可以抽象為如下算子：

V=T（σ），（5）

其中σ=（σ1，…，σd）T為波動率參數向量， V=（U0，…，UM）為期權離散矩陣， T為數值求解算子.

2"反問題描述和Bayes算法

在式（5）中， 如果已知期權離散矩陣V和G， 求波動率參數向量σ， 則稱為反問題. 該類問題是有確定性的， 通常也是病態(tài)的，

一般用Tikhonov正則化方法^［9］求解該類問題. 但從統(tǒng)計學的角度， 更關注的是真解的分布， 而不再是一個點. 式（5）的統(tǒng)計模型版本可寫為

V=T（σ）+δ，（6）

其中δ={δij}^{j=1，2，…，M+1}i=1，2，…，N， δij～δ·N（0，1）， N（0，1）是均值為0、 方差為1的標準正態(tài)分布.

2.1"Bayes公式

Bayes公式為

πV（σ）∶=π（σV）=

π（Vσ）π0（σ）∫Zπ（Vσ）π0（σ）dσ，（7）

其中πV（σ）和π0（σ）分別是波動率后驗和先驗的概率密度函數. 于是有

πV（σ）=πη（V-T

（σ））π0（σ）∫Zπη（V-T（σ）π0（σ）

dσ∝πδ（V-T（σ））π0（σ），（8）

其中∝表示兩者成比例的關系， πδ（·）為最大似然估計， 并通?？紤]先驗為正態(tài)分布.

在線性正演模型和高斯先驗的假設下， 最大似然估計等價于Tikhonov正則化^［10］. 從確定性的角度可見， 通過優(yōu)化框架只能得到反問題相關的一個點估計.

而從統(tǒng)計學的角度， 可根據Bayes公式獲得更多的統(tǒng)計信息. 而且這個后驗分布， 使得反演層參數的刻畫更全面并包含不確定性^［11］

. 本文采用Metropolis-Hastings（MH）算法， 從后驗分布πV（σ）中進行采樣.

2.2"Metropolis-Hastings算法

Markov Chain Monte Carlo （MCMC）算法是一類用于從復雜分布中抽樣的統(tǒng)計方法， 其中MH算法是常用的MCMC方法之一^［12］.

算法2^＼V（σ）， 一個選定分布q， 最大樣本數 M1和一個初始點σ0;

輸出： 后驗樣本集合{σm}^M1m=1.

for m=1∶M1

提議一個候選點： σc～q（·σm-1）;

計算接受概率： α（σc，σm-1）=min1，πV

（σc）q（σm-1σc）πV（σm-1）q

（σcσm-1）;

生成一個樣本u～U（0，1）;

if ult;α（σc，σm-1）

σm=σc;

else

σm=σm-1;

end if

end for

雖然MH采樣算法可以計算真解的后驗分布， 但無論是Tikhonov正則化方法， 還是MH采樣算法， 都需要求解大量的正問題（5）. 雖然每個正問題求解只需幾秒鐘， 但計算大量正問

題， 會增加整個算法的計算時間. 本文用神經網絡的方法構造替代模型， 以減少整體算法的計算量.

2.3"神經網絡替代模型

全連接神經網絡（FNN）是一種特殊結構的神經網絡， 它由一系列全連接層組成， 對非線性函數可以有效地逼近. 全連接神經網絡由（L+1）層網絡組成， 其中前一層的節(jié)點x^l-1與

下一層的所有節(jié)點xl都有連接. 在每個節(jié)點上， 有兩個函數依次復合構成. 第一個是一次函數

L^l（x^l-1）=Wlx^l-1+bl：^nl-1→^nl， 第二個是激活函數σl：^nl→^nl（例如Sigmoid，tanh，ReLU）， 其中Wl表示權重， bl表示偏置. 圖1為FNN的示意圖. 為簡便，

令Θ={Wl，bl}Ll=1表示FNN的所有網絡參數， 則可以定義一個由Θ參數化的映射FNN：m→n：

y=FNN（z;Θ）=L^Lσ^L-1L^L-1…σ1L¹（z）.（9）

用式（9）求解正問題（5）， 即z=σ， y=V， FNN≈T， m=d， n=N×（M+1）. 設訓練集為{（σi，Vi）}ⁿ¹i=1， 測試集為{（σi，Vi）}ⁿ²

i=1， n*=n1+n2表示樣本集的總數. 對于監(jiān)督學習， 考慮平方損失函數， 即

Loss（Θ）=12n1∑n1i=1‖Vi-FNN（σi;Θ）‖2.（10）

本文主要目標是通過適當的優(yōu)化算法， 找到最優(yōu)參數Θ*， 即Θ*=argminΘ Loss（Θ）.（11）

解決問題（10）的常見優(yōu)化算法目前已有很多， 如梯度下降法、 牛頓法和隨機梯度下降法等. 對本文所考慮的問題， 通常n1很大， 導致計算梯度的代價很高. 為減少相應的計算量， 這里使用最小批量梯度下降方法. 該方法將訓練集劃分為小批次， 用于計算損失函數誤差和梯度， 并更新Θ， 綜合了隨機梯度下降法和梯度下降法的優(yōu)勢.

關于替代模型， 本文有如下收斂性結果. 假設：

（H1） 散度函數f滿足

f（t）≤C1log t，0lt;tlt;1，C2tlog t，t≥1，

其中C1和C2是依賴于f的兩個常數.

引理2^［11］"對εgt;0， 存在C=C（ε）gt;0， 使得對任意σ∈d， 有

‖T（σ）-FNN（σ;Θ）‖≤Cexp{ε‖σ‖2}ΦT（Θ），

其中ΦT（·）表示與T相關的任何函數， 只需滿足lim#{Θ}→∞ ΦT（Θ）=0， #{Θ}表示神經網絡參數的個數.

設μ0為先驗測度， μ為離散算子T（·）對應的真實模型后驗測度， μF為神經網絡FNN（·;Θ）替代的后驗測度， 下列定理給出了在f散度下兩者誤差的上界.

定理1"假設（H1）成立， 設測度μ，μFμ0， T（·）和FNN（·;Θ）都是

Lipschitz函數. 則測度μ和μF的誤差在f散度下是收斂的， 即存在常數Cgt;0， 使得

Df（μ‖μF）∶=∫dfdμ1/dμ0dμ2/dμ0

dμ2dμ0dμ0（σ）≤CΦG（Θ）.（12）

定理1的證明類似文獻［11］， 首先利用正演模型和替代模型是Lipschitz函數的性質， 以及相應逼近的假設， 并基于Fernique定理， 得到歸一化常數的估計. 然后， 根據散度函數f具有分段上界的假設， 以及正演模型和替代模型的假設， 將f散度進行分段估計， 從而可得結論. 定理1表明， 當本文替代神經網絡模型與正演模型足夠逼近時， 在f散度意義下相應的替代后驗測度收斂到真實模型的后驗測度， 故替代模型有一定的理論保障.

3"數值實例

下面進行數值實驗（d=1，3）驗證本文提出算法的有效性. 對于正問題（5）， 選擇T=1， r=0.03， q=0.01， 并將轉換后的截斷長度設為L=1.6， 以確保截斷區(qū)域足夠大. 對于離散化問題， 時間和空間分割分別為M=150和N=100.用MH和MH+FNN算法針對式（6）進行采樣， 并設先驗分布π0（σ）～Nd（-1.8×1d，0.1×Id）， 取密度q（·;σ）～Nd（σ，10^-6×Id）. 此外， 采樣噪聲振幅為δ=0.012. 在產生σ的若干樣本后， 可利用算法1對現性互補問題（4）進行求解， 得到正問題的若干組解（即期權離散矩陣V）， 并作為監(jiān)督學習目標泛函中的真解. 然后， 以σ的若干樣本為輸入參數， 正問題的若干組解為輸出， 經過優(yōu)化訓練，得到神經網絡FNN的參數， 并用該網絡代替正問題的求解.

取神經網絡結構為{d，50，50，50，50，151×100}， d=1，3. 使用Adam優(yōu)化器對網絡進行訓練， 每次迭代使用128個樣本的小批量， 并進行5 000次循環(huán)的訓練， 而權重初始化通過Glorot均勻初始化.這里， 每次使用原始對偶活躍集方法， 對式（6）進行一次求解， 可得一個期權離散矩陣V， 即為一個樣本. 取總樣本數n*=20 000， 其中訓練集的樣本數n1=16 000， 測試集的樣本數n2=4 000. 所有的數值實驗都在一臺Intel Core i7的計算機工具包MATLAB R2020b上實現.

例1"考慮如下一維的波動率：

σ（t）=σa（t），""a（t）=0.3et.（13）

表1列出了當σ=0.15，0.3，0.45時， Bayes方法（MH）和Bayes+替代模型（MH+FNN）的數值結果， 其中時間包括線下時間和線上時間. 這里， FNN的實驗輸出的是期權價

格V， 但MH和MH+FNN實驗輸出的是波動率σ. 由表1可見， 兩種方法的精度近似， 但后者花費的時間只有前者的1/3左右， 故神經網絡作為替代模型， 加速效果顯著.

表2列出了MH+FNN算法固定每層神經元為50時， 不同隱藏層的數值結果. 表3列出了MH+FNN算法固定隱藏層為3時， 不同神經元個數的數值結果. 表2和表3的

數值結果表明， 隨著隱藏層和神經元個數的增加， 精度會有些提高， 但會增加一定的時間. 對于增加樣本， 有類似的結果， 這里不再列出.

例2"考慮含3個參數的波動率， 即σ（t）=∑3i=1σiai（t），其中a1（t）=1， a2（t）=t， a3（t）=3^-t， 波動率的真解為σ1=0.1， σ2=0.15， σ3=0.2. 表4列出了（σ1，σ2，σ3）=（0.1，0.15，0.2）時， Bayes方法（MH）和Bayes+替代模型（MH+FNN）的數值結果. 與一維的結果類似， MH+FNN的計算結果和MH類似， 但時間上要快2/3左右.

圖2給出了MH和MH+FNN算法采樣的相關圖. 由圖2可見： （A），（C），（F）是自相關圖， 其兩條線基本重合， 說明自相關度非常高; （B），（D），（E）是協(xié)相關圖， 其各自的兩個環(huán)

線也基本重合， 說明協(xié)相關度很高. 實驗結果表明， 兩種采樣效果一致.

參考文獻

［1］"ACHDOU Y， INDRAGOBY G， PIRONNEAU O.

Volatility Calibration with American Options ［J］. Methods Appl Anal， 2004， 11（4）： 533-556.

［2］"BURKOVSKA O， GLAU K， MAHLSTEDT M， et al. Complexity Reduction for Calibration

to American Options ［J］. J Comput Finance， 2019， 23（1）： 25-60.

［3］"KUTNER G W. Determining the Implied Volatility for A

merican Options Using the QAM ［J］. "Financ Rev， 1998， 33（1）： 119-130.

［4］"HAN J Q， JENTZEN A， E W N. Solving High-Dimensional Partial Differential Equations Using Deep L

earning ［J］. Proc Natl Acad Sci USA， 2018， 115（34）： 8505-8510.

［5］"LIU S Q， LEITAO ， BOROVYKH A， et al. On a Calibration Neural Network to Extract Implied Information fr

om American Options ［J］. Appl Math Finance， 2020， 28（5）： 449-475.

［6］"BLACK F， SCHOLES M. The Pricing of Options and Corporate Liabilities ［J］. J Polit Econ， 1973， 81（3）： 637-654.

［7］"GAO Y， SONG H M， WANG X S， et al. Primal-Dual Activ

e Set Method for Pricing American Better-of Option on Two Assets ［J］. Commun Nonlinear Sci Numer Simulat， 2020， 80： 104976-1-104976-15.

［8］"SONG H M， WANG X S， ZHANG K， et al. Primal-Dual Active Set Method for American Lookback Put Option Pricing ［J］. "East Asian J Appl Math， 2017， 7（3）： 603-614.

［9］"BOUCHOUEV I， ISAKOV V. The Inverse Problem of Option Pricing ［J］. Inverse Probl， 1997， 13（5）： L11-L17.

［10］"LEHTINEN M S， PIVRINTA L， SOMERSALO E. Linear Inverse Problems for Generalised Random Variables ［J］.

Inverse Probl， 1989， 5（4）： 599-612.

［11］"GAO Y， LIU H Y， WANG X C， et al. A Bayesian Scheme

for Reconstructing Obstacles in Acoustic Waveguides ［J］. J Sci Comput， 2023， 97（3）： 53-1-53-45.

［12］"HASTINGS W K. Monte Carlo Sampling Methods Using Markov Chains and Their Applications ［J］. Biometrika， 1970， 57（1）： 97-109.

（責任編輯： 李"琦）