楊燕曦

（湖南省直黨校 經(jīng)濟(jì)學(xué)教研室，湖南 長(zhǎng)沙 410079）

1 引 言

隨著科學(xué)技術(shù)的不斷發(fā)展，在現(xiàn)代金融工程領(lǐng)域愈來(lái)愈重視定量的數(shù)理分析，大量的實(shí)際問(wèn)題，如動(dòng)態(tài)最優(yōu)定價(jià)、金融衍生產(chǎn)品的定價(jià)、投資風(fēng)險(xiǎn)的規(guī)避等，經(jīng)過(guò)數(shù)理建模，最終都?xì)w結(jié)為對(duì)隨機(jī)微分方程（組）的求解［1－3］．這些微分方程（組）中很多都不易求得解析解，發(fā)展相應(yīng)的數(shù)值解法具有重大意義．傳統(tǒng)的數(shù)值求解方法主要包括二叉樹(shù)方法，蒙特卡洛方法、有限差分法［4］，這些方法對(duì)計(jì)算機(jī)的計(jì)算能力要求較低，計(jì)算精度不高．近年來(lái)，國(guó)內(nèi)外學(xué)者又將有限元法應(yīng)用于金融工程計(jì)算領(lǐng)域［5］，提高了計(jì)算的精度和效率，但其收斂性和穩(wěn)定性還有待進(jìn)一步研究．當(dāng)前，金融活動(dòng)的風(fēng)險(xiǎn)及復(fù)雜性進(jìn)一步加劇，數(shù)理建模得到的微分方程規(guī)模更大、復(fù)雜程度更高，有的還具有一定的非線性，迫切需要一種簡(jiǎn)潔、準(zhǔn)確、高效的數(shù)值計(jì)算方法．

求積元方法是一種結(jié)合了高效數(shù)值積分和微分求積法二者優(yōu)勢(shì)的新的求解常（偏）微分方程（組）的高階數(shù)值方法．該方法自2007年由清華大學(xué)鐘宏志教授提出以來(lái)，在工程結(jié)構(gòu)分析領(lǐng)域中已得到較為廣泛地應(yīng)用［6－9］，展現(xiàn)出其相比傳統(tǒng)有限元法的獨(dú)特優(yōu)勢(shì)．

工程結(jié)構(gòu)計(jì)算分析所涉及的微分方程（組）一般均具有線性自伴隨的特性，因而具有相應(yīng)的變分形式．而對(duì)于金融工程計(jì)算分析中所涉及的微分方程（組）一般不具有自伴隨的特性，對(duì)于求積元方法的應(yīng)用還是一個(gè)新的領(lǐng)域．

針對(duì)金融工程計(jì)算領(lǐng)域的靜態(tài)一維問(wèn)題，將求積元方法應(yīng)用于非自伴隨的微分方程的數(shù)值求解，建立相應(yīng)的求積元單元．選取3個(gè)典型問(wèn)題進(jìn)行計(jì)算，與解析解、有限差分解和有限元解分別進(jìn)行比較，驗(yàn)證求積元方法的適應(yīng)性、準(zhǔn)確性和高效性．為該方法在金融工程計(jì)算領(lǐng)域動(dòng)態(tài)問(wèn)題（期權(quán)定價(jià)問(wèn)題）、二維問(wèn)題中的深入應(yīng)用奠定基礎(chǔ)．

2 一維邊值問(wèn)題的求積元離散

一般地，金融工程中的靜態(tài)一維問(wèn)題可用如下微分方程

和相應(yīng)的邊界條件表示，

式（1）中，u（x）為定義在區(qū)域）上的未知（待求）函數(shù)，（x）（x）分別表示對(duì)x求二階、一階導(dǎo)數(shù)．a(chǎn)1（x）、a2（x）、f（x）為已知函數(shù)．式（2）、式（3）為邊界條件．

假設(shè)未知函數(shù)u（x）可以用近似函數(shù)u（x）來(lái)表示，基于Galerkin加權(quán)殘值積分近似為零和求積元法求解思想，權(quán)函數(shù)選定為近似函數(shù)u的變分δu，令式（1）殘值在加權(quán)積分意義下為零，即

對(duì)式（4）中的二階導(dǎo)數(shù)進(jìn)行分部積分

式（5）中，b．t．?（）u表示邊界條件．

將式（5）中積分進(jìn)一步離散，根據(jù)求積元求解基本步驟，首先將待求解物理域坐標(biāo)系通過(guò)式（6）轉(zhuǎn)換到標(biāo)準(zhǔn)域，如圖1所示，圖中1，2，3，…，N－1，N為L(zhǎng)obatto數(shù)值積分［10］點(diǎn)．

圖1 物理域到標(biāo)準(zhǔn)域映射

利用Lobatto數(shù)值積分［10］計(jì)算式（5）中的積分，

其中，N表示積分點(diǎn)數(shù)，右側(cè)下標(biāo)i表示該變量在積分點(diǎn)處的值，Hi為相應(yīng)積分點(diǎn)對(duì)應(yīng)的積分權(quán)系數(shù)．需指出，式（7）中導(dǎo)數(shù)均為對(duì)標(biāo)準(zhǔn)域坐標(biāo)ξ求導(dǎo)．結(jié)合微分求積法則［11］，

將式（7）中所含積分點(diǎn)處的函數(shù)值和函數(shù)導(dǎo)數(shù)值表示為積分點(diǎn)處基本自由度（近似函數(shù)值）的線性加權(quán)代數(shù)和．式（8）中，為m階微分求積系數(shù)．

物理域坐標(biāo)系下Lobatto數(shù)值積分點(diǎn)處組成的列向量?u構(gòu)成了待求解問(wèn)題的單元基本自由度，

右上角（e）即表示一個(gè)求積元單元，則

式（10）中，

其中為Kronecker符號(hào)，即

則式（7）可進(jìn)一步表示為

則式（13）中

則式（5）最終離散為

由于變分具有任意性，式（15）可轉(zhuǎn)化為一個(gè)線性代數(shù)方程組，

對(duì)于邊界條件b．t．?（）u，當(dāng)β1≠0且β2≠0時(shí)，邊界條件可表示為

可對(duì)矩陣K、F修正如下：其余元素分別與K、F一致，則式（16）轉(zhuǎn)化為

進(jìn)行求解．

當(dāng)β1＝0且β2＝0時(shí)，邊界條件可表示為

由于β1＝0且β2＝0，由式（2）和（3）可知，u1、uN為常量，

則

只需修正K、F，使其滿足式（21）即可．故修正如下：

、其余元素分別與K、F，則式（15）仍轉(zhuǎn)化為

進(jìn)行求解．其余邊界條件，如β1≠0而β2＝0，亦可類(lèi)似處理．

3 實(shí)證分析

選取金融工程計(jì)算分析中較為典型的3個(gè)實(shí)例，采用求積元方法進(jìn)行計(jì)算，驗(yàn)證求積元方法的準(zhǔn)確性和高效性．計(jì)算程序采用Matlab軟件編制．

3．1 壟斷動(dòng)態(tài)最優(yōu)化問(wèn)題

壟斷企業(yè)的目標(biāo)是尋找產(chǎn)品價(jià)格P的一條最優(yōu)路徑，從而在一個(gè)有限的時(shí)間內(nèi)［0，T］內(nèi)實(shí)現(xiàn)利潤(rùn)最大化．假設(shè)這個(gè)時(shí)期足夠短，以保證固定的需求成本函數(shù)以及忽略折現(xiàn)的設(shè)定是合理的．這個(gè)問(wèn)題可以通過(guò)變分法采用一個(gè)歐拉方程來(lái)描述［12］，

該方程是一個(gè)二階線性微分方程，其解析解為

將邊值條件代入式（26），可得

應(yīng)用求積元方法對(duì)該問(wèn)題在t＝［0，2］定義域內(nèi)進(jìn)行求解，各時(shí)刻t價(jià)格P的計(jì)算結(jié)果與解析解的對(duì)比如表1所示．計(jì)算相關(guān)參數(shù)：產(chǎn)出函數(shù)中的系數(shù)，a＝160，b＝8，h＝100；總成本函數(shù)中的系數(shù)，α＝0．1，β＝100；P0＝11，PT＝15．由表1可見(jiàn)求積元方法僅需劃分1個(gè)求積元單元4個(gè)積分點(diǎn)（N＝4）共計(jì)4個(gè)自由度即可達(dá)到良好的求解精度，小數(shù)點(diǎn)后4位有效數(shù)字與解析解完全一致，體現(xiàn)出求積元方法的準(zhǔn)確性．

表1 價(jià)格P計(jì)算結(jié)果對(duì)比（1個(gè)單元4自由度）

3．2 幾何布朗運(yùn)動(dòng)的首出時(shí)

考察幾何布朗運(yùn)動(dòng)

在給定標(biāo)的物價(jià)格范圍內(nèi)的首出時(shí)是有實(shí)踐意義的．可以得到給定標(biāo)的物價(jià)格偏離某一確定界限的平均時(shí)間，進(jìn)而評(píng)估相關(guān)雙障礙期權(quán)的風(fēng)險(xiǎn)．該問(wèn)題可以描述為

該方程的解析解為

應(yīng)用求積元方法對(duì)該問(wèn)題進(jìn)行求解，計(jì)算結(jié)果與解析解及有限元解［5］的對(duì)比如表2所示．計(jì)算相關(guān)參數(shù)為收益率a＝0．1，波動(dòng)率σ＝0．2，xmin＝20，xmax＝60．由表2可見(jiàn)求積元僅需劃分1個(gè)單元23個(gè)積分點(diǎn)共計(jì)23個(gè)自由度即可達(dá)到良好的求解精度，小數(shù)點(diǎn)后8位有效數(shù)字與解析解完全一致，而有限元法則需要?jiǎng)澐?9個(gè)單元共計(jì)200個(gè)自由度才能達(dá)到以上精度，求積元法的計(jì)算自由度僅約為有限元法的十分之一，而計(jì)算大規(guī)模問(wèn)題時(shí)，計(jì)算自由度是影響計(jì)算機(jī)計(jì)算效率的重要因素．因此．求積元法相比有限元法具有更為高效的特點(diǎn)．

表2 首出時(shí)u計(jì)算結(jié)果對(duì)比

3．3 對(duì)流占優(yōu)問(wèn)題

對(duì)流占優(yōu)問(wèn)題在金融工程中具有很強(qiáng)的實(shí)際意義［13］，比如當(dāng)標(biāo)的物價(jià)格較低且（／或）波動(dòng)率較低時(shí)，股票期權(quán)、外匯期權(quán)的定價(jià)將成為對(duì)流占優(yōu)問(wèn)題．以如下的邊值問(wèn)題

為例進(jìn)行說(shuō)明，當(dāng)k減小時(shí)，該微分方程橢圓型方程特征逐漸減弱，雙曲型方程特征逐漸增強(qiáng)．此時(shí)，由于“對(duì)流項(xiàng)”主要影響方程的特性，該問(wèn)題稱(chēng)為對(duì)流占優(yōu)問(wèn)題．

該方程的解析解為

應(yīng)用求積元方法對(duì)該問(wèn)題進(jìn)行求解，計(jì)算結(jié)果與解析解及有限差分解的對(duì)比如圖2所示．計(jì)算相關(guān)參數(shù)為k＝0．002．由圖2可見(jiàn)，該問(wèn)題的解析解曲線具有很強(qiáng)的非線性，表現(xiàn)為在［0，0．99］范圍內(nèi)非常平緩，而在［0．99，1］范圍內(nèi)急劇上升．

本例中求積元方法（QEM）共劃分8個(gè)單元，每個(gè)單元采用15個(gè)積分點(diǎn)，共計(jì)113個(gè)自由度，達(dá)到了較好的計(jì)算結(jié)果．而有限差分法（FDM）在劃分單元數(shù)較少時(shí)，計(jì)算結(jié)果出現(xiàn)了明顯的震蕩［5］，即使劃分200個(gè)單元（201個(gè)自由度），也存在震蕩現(xiàn)象（如圖2所示）．若采用有限元方法，得到滿意的計(jì)算結(jié)果也需要200個(gè)自由度以上［5］．相比有限差分法和有限元法，求積元法的計(jì)算自由度縮減了近一半，再次體現(xiàn)出準(zhǔn)確高效的特點(diǎn)．

圖2 求積元法計(jì)算結(jié)果與解析解和有限差分解對(duì)比

4 結(jié) 論

針對(duì)金融工程計(jì)算領(lǐng)域的靜態(tài)一維問(wèn)題，將求積元方法的應(yīng)用領(lǐng)域從線性自伴隨微分方程的求解拓展到非自伴隨微分方程的求解．首先，基于Galerkin加權(quán)殘值法思想建立了相應(yīng)的求積元單元；之后，選取了三個(gè)典型問(wèn)題進(jìn)行編程求解計(jì)算，并與解析解、有限差分解和有限元解分別進(jìn)行了比較．

計(jì)算結(jié)果表明，相比有限元方法和有限差分法，求積元方法在得到相同精度計(jì)算結(jié)果的同時(shí)，大幅減少了自由度數(shù)，提高了計(jì)算效率．對(duì)于一般性問(wèn)題，僅需劃分一個(gè)單元，也可視問(wèn)題的復(fù)雜性進(jìn)行多單元拼接求解．是一種準(zhǔn)確、高效和靈活的數(shù)值方法．用于金融工程領(lǐng)域的靜態(tài)一維問(wèn)題計(jì)算分析有較大的優(yōu)勢(shì)，可進(jìn)一步用于該領(lǐng)域動(dòng)態(tài)問(wèn)題（期權(quán)定價(jià)問(wèn)題）、二維問(wèn)題的計(jì)算分析．

［1］M ROSS．An Elementary Introduction to Mathematical Finance．［M］．London：Cambridge University Press，2011．

［2］郭宇權(quán)．金融衍生產(chǎn)品數(shù)學(xué)模型［M］．北京：世界圖書(shū)出版公司北京公司，2010．

［3］科森多爾．隨機(jī)微分方程［M］．北京：世界圖書(shū)出版公司北京公司，2006．

［4］蔣致遠(yuǎn)，張?zhí)?，龔閃閃．基于拉普拉斯變換有限差分方法的B－S期權(quán)定價(jià)［J］．經(jīng)濟(jì)數(shù)學(xué)，2014，31（3）：18－22．

［5］T JURGEN．Financial engineering with finite elements［M］．Chichester：John Wiley ＆Sons Ltd，2005．

［6］H ZHONG，T YU．A weak form quadrature element method for plane elasticity problems［J］．Applied Mathematical Modelling，2009，33（10）：3801－3814．

［7］H ZHONG，M GAO．Quadrature element analysis of planar frameworks［J］．Archive of Applied Mechanics，2010，80（12）：1391－1405．

［8］Z SHEN，H ZHONG．Static and vibrational analysis of partially composite beams using the weak－form quadrature element method［J］．Mathematical Problems in Engineering，Vol．2012，Article ID 974023，23pages．

［9］Z SHEN，H ZHONG．Nonlinear quadrature element analysis of composite beams with partial interaction．Australian Journal of Mechanical Engineering，2013，11（1）：45－52．

［10］P DAVIS，P RABINOWITZ．Methods of numerical integration．［M］．Orlando：Academic Press，1984．

［11］C SHU．Differential quadrature and its application in engineering［M］．London：Springer－Verlag，2000．

［12］A CHIANG．Elements of Dynamic Optimization［M］．New York：McGraw－Hill，Inc，1992．

［13］R SEYDEL．Tools for Computational Finance［M］．Berlin：Springer，2002．