跳躍—擴散模型資產(chǎn)定價公式的數(shù)值計算方法

摘 要 假定資產(chǎn)價格變化過程服從跳躍擴散過程，那么基于它的歐式期權就滿足一個偏積分微分方程(PIDE)，本文利用差分法來離散這個PIDE方程，用兩種迭代方法得到方程的數(shù)值解:基于雅可比正則分裂法和預條件共軛梯度法.

關鍵詞 跳躍擴散模型;差分法;FFT算法;歐式看漲期權

中圖分類號 O241.82文獻標識碼:A

1 引 言

美國芝加哥大學教授Black和Scholes^[1]在1973年發(fā)表了“The Pricing of Options and Corporate Liabilities”一文，提出了著名的BlackScholes期權定價公式，在BS公式中，假設股票的價格過程是連續(xù)的幾何布朗運動，但是，在現(xiàn)實市場中，一些突發(fā)情況會引起股票價格發(fā)生跳躍.基于上述考慮，Merton在1976年首先提出了跳躍擴散模型，在Merton模型中，資產(chǎn)價格在沒有受到外界重大影響時服從布朗運動，當資產(chǎn)價格受到突發(fā)事件的影響而發(fā)生跳躍時，就用跳躍過程來描述.

本文首先介紹PIDE^[2]的具體形式，在Merton模型下對其差分離散，得到一個Toeplitz矩陣方程，用兩種方法解這個矩陣方程，一是基于雅可比正則分裂的迭代方法^[3-4]，二是預條件共軛梯度方法.考慮到Toeplitz^[5]矩陣的特殊性，在迭代的過程中，將其植入到一個循環(huán)矩陣中，利用循環(huán)矩陣和向量的乘積來計算Toeplitz和向量的乘積，而循環(huán)矩陣向量乘積可以通過快速富里葉變換(FFT)快速計算，這樣就加快了迭代速度.共軛梯度法是解決Toeplitz線性方程組的主要方法之一，在利用共軛梯度法的情況下，快速傅里葉變換的作用是加快共軛梯度法的迭代速度，但不改變其收斂速度，共軛梯度法的收斂速度取決于線性方程組系數(shù)矩陣的條件數(shù)，基于此考慮，本文采用預條件共軛梯度算法，選用R.Chan優(yōu)化循環(huán)預條件器^[6]，預條件器的使用是為了改善系數(shù)矩陣的條件數(shù)，以便提高收斂速度.

2 跳躍擴散模型

假設市場是完備無套利的市場，在跳躍擴散模型下，資產(chǎn)的價格變化過程服從隨機微分方程:dS(t)=S(t)(υ(t)dt+σ(t)dω(t)+(η-1)dq(t)).(1)

其中，υ(t)是漂移率，σ(t)是波動率，ω(t)標準布朗運動，dq(t)是泊松過程，dq(t)=0的概率是1-λdt，dq(t)=1的概率是λdt，λ是泊松到達強度，η-1是由S跳躍到Sη的跳躍幅度函數(shù)，是一個隨機變量，用ζ表示平均跳躍幅度E(η-1)，泊松過程dq(t)與布朗運動ω(t)是相互獨立的.

由文獻[7]可知，在上述假設下，基于資產(chǎn)價格S與時間τ的未定權益V(S，τ)滿足PIDE:

Vt=σ2S22VSS+(r-λζ)SVS-(r+λ)V+λ∫

SymboleB@ 0V(Sη)g(η)dη.(2)

式中，t=T-τ是到期時間為T的時間，r是無風險利率，g(η)是跳躍幅度η的密度函數(shù).

對式(2)的積分部分進行指數(shù)變量變換，令S=ex，η=ey，則式(2)變?yōu)楠И?/p>

Vt=σ2S22VSS+(r-λζ)SVS-(r+λ)V+λ∫

SymboleB@ 0V(Sη，t)g(η)dη.(3)

再對其余部分進行變換，令f(y)=g(ey)ey，υ(x+y，t)=V(ex+y，t)，函數(shù)f(y)是跳躍幅度y=log(η)的密度函數(shù)，則式(3)變?yōu)楠И?/p>

υt=σ22υxx+(r-12σ2-λζ)υx-(r+λ)υ+λ∫

SymboleB@ -

SymboleB@ f(y)υ(x+y，t)dy，

(t，x)∈[0，T)×R，邊界條件υ(T，x)=g(ex).(4)令u(τ，#8226;)=υ(T-τ，#8226;)，則式(4)變?yōu)楠И?/p>

uτ-12σ2uxx-(r-12σ2-λζ)ux+(r+λ)u-λ∫

SymboleB@ -

SymboleB@ u(τ，x+y)f(y)dy=0，

(τ，x)∈[0，T)×R， u(0，x)=g(ex).(5)

3 Merton模型下的有限差分離散

在PIDE中，由于S=ex，則x=ln(S)，通常限定x的取值范圍是Ω=(-x，x)，x*稱為截斷點，Ω*稱為截斷區(qū)域，式(5)中的積分部分可以分割為兩部分∫R=∫Ω+∫R\\Ω.在計算歐式看漲期權的情況下，在R\\Ω上的積分u(τ，z)可以進行近似代替:

u(τ，x)→ex-Ke-rτ， 當x→+

SymboleB@ 時.(6)

u(τ，x)→0， 當x→-

SymboleB@ 時.(7)經(jīng) 濟 數(shù) 學第 27 卷

第2期張鴻雁等:跳躍擴散模型資產(chǎn)定價公式的數(shù)值計算方法

對于式(5)中的積分部分，進行變量變換z=x+y，則

∫Ru(τ，x+y)f(y)dy=∫Ru(τ，z)f(z-x)dz.(8)

定義函數(shù):

ε(τ，x，x)=∫Ω(ez-Ke-rτ)f(z-x)dz.(9)

在Merton模型中:

f(x)=12πσJe-(x-μJ)2/2σ2J.(10)

在μJ=0的情況下，有表達式:

ε(τ，x，x)=ex+σ2J2(x-x+σ2σJ)-Ke-rτ(x-xσJ).(11)

其中，(y)是標準分布函數(shù):

(y)=12π∫y-

SymboleB@ e-x22dx.(12)

考慮時間和空間節(jié)點，使xi=-x*+(i-1)h，i=1，…，n，

h=(x-(-x))/(n-1)=2x/(n-1)，和τm=(m-1)k，m=1，…，q，

k=T/q.記umi=u(τm，xi)，fij=f(xj-xi).

在[-x，x]上用復合梯形原則，有積分近似:

∫Ru(τm，z)f(z-xi)dz≈h2[um1fi，1+2∑n-1j=2umjfi，j+umnfi，n]+ε(τm，xi，x) ，i=2，…，n-1.(13)

對于時間變量與空間變量，作近似:

uτ(τm，xi)≈(32umi-2um-1i+12um-2i)，m≥2，

(umi-um-1i)/k， m=1.(14)

uxx(τm，xi)≈(umi+1-2umi+umi-1)/h2，(15)

ux(τm，xi)≈(umi+1-umi-1)/2h，(16)

這里，對于時間的偏導，當m≥2時，用向后的二階差分來近似對時間的微分，當m=1時，用向后的一階差分近似;對于空間的偏導，用中心差分來近似.

定義向量um=(um1，…，umn)T.由初始條件，初始向量:

u1=(g(ex1)，…，g(exn))T.

由式(13)~(16)，則式(5)的有限差分離散能寫成矩陣形式:

Aum=bm，(17)

其中，A=γ0I+C+D，γ0=1， m=1，

3/2，m≥2，I是單位矩陣，C和D定義為

cij=-kσ2/2h2+k(r-σ2/2-λζ)/2h，如果i=j-1， 2≤i≤n-1，kσ2/h2+(r+λ)k，如果i=j， 2≤i≤n-1，-kσ2/2h2-k(r-σ2/2-λζ)/2h，如果i=j+1， 2≤i≤n-1，

0，其他情形;

dij=-khλfij/2，如果2≤i≤n-1，且j=1，n，-khλfij如果2≤i≤n-1，且2≤j≤n-1，0，其他情形，且式(17)中，向量bm=(b1，b2，…，bn)T定義為:

bi=kλε(τm，xi，x)+γ1um-1i+γ2um-2i， i=2，…，n-1，(18)

其中

γ1=1，如果m=1，2，如果m≥2，(19)

γ2=0，如果m=1，-1/2，如果m≥2.(20)

由初始邊界條件可知:b1=0，bn=γ0(ex-Ke-rτm).

4 基于雅可比正則分裂的迭代方法

定義1 假設矩陣A可用分裂成形式:

A=Q-R，(21)

其中，Q是單調(diào)矩陣(Q-1≥0)且R≥0，則稱A可以正則分裂.

對于每一個形如式(21)的分裂，都存在相應的迭代方法:

Vl+1=Q-1RVl+Q-1b，l=0，…，V0=0.(22)

若A是單調(diào)矩陣，則迭代式(22)是收斂的(ρ(Q-1R)<1)，證明過程見文獻[8].

給出雅可比正則分裂的形式:

(A)A=Q1-R1，其中Q1是A的對角矩陣.

如果滿足:

(i)-(k/h)σ2/2-k(λζ+σ2/2)/2+h(ω0+λk)/6≤0，

(ii)-(k/h)σ2/2+k(λζ+σ2/2)/2+h(ω0+λk)/6≤0，

則分裂(A)是正則的，且

ρ(Q-12R2)≤ρ(Q-11R1)<1，(23)

證明過程見文獻[9].

在有限差分法中，若:

(iii)-(k/h)σ2/2-k(λζ+σ2/2)/2≤0，

(iv)-(k/h)σ2/2+k(λζ+σ2/2)/2≤0，

則可以得到一個精確穩(wěn)定的解.

若保持k/h固定不變而讓h→0，則存在一個h0>0使得在h≤h0時條件(i)~(iv)同時成立.

5 預條件共軛梯度算法

本文中系數(shù)矩陣A是一個Toeplitz矩陣，現(xiàn)選擇R.Chan優(yōu)化循環(huán)預條件器^[10]加快迭代過程中的收斂速度.預條件器C:

C=c0cn-1…c2c1c1c0…c3c2

cn-2cn-3…c0cn-1cn-1cn-1…c1c0.(24)

其中，cj=(n-j)aj+jaj-nn，aj是矩陣A中的元素，j=0，…，n-1.

在所有的n階循環(huán)矩陣中，C極小化Frobenius范數(shù)‖C-T‖F，在這個意義下，C被視為A的一個近似矩陣.在預條件共軛梯度法下，每一次迭代都要計算矩陣向量積Ax和C-1y(x和y是n維向量)，可以利用快速富里葉變換(FFT)快速計算.C可以被n階離散富里葉矩陣對角化，即

C=F*n∧Fn，∧=diag[λ0，λ1，…，λn-1]，(25)

且λk=∑n-1j=0cje-i2πjk/n(k=0，1，…，n-1)是C的特征值(i是虛數(shù)單位)，于是

C-1y=IFFT(∧-1#8226;FFT(y)).(26)

對于計算Ax，可以先將A嵌入到一個2n階的循環(huán)矩陣，然后借助2n階快速富里葉變換來計算，即

ABBAx0=AxBx.(27)

其中B=0an-1…a1a1-n0…a2

a-2a-3…an-1a-1a-2…0.(28)

Aum=bm等價于C-1Aum=C-1bm

對于Toeplitz矩陣方程Aum=bm，循環(huán)優(yōu)化預條件器是式(24)，共軛梯度法采用析因形式^[11]，不生成系數(shù)矩陣，迭代算法為:

r(0)=b-AC-1u(0)，u(0)是初值;(29)

p(0)=s(0)=C-*A*r(0);(30)

γ0=‖s0‖22，

for k=0，1，…，

q(k)=AC-1p(k)，

αk=γk/‖q(k)‖22，

u(k+1)=u(k)+αkp(k)，

r(k+1)=r(k)-αkq(k)，

s(k+1)=C-*A*r(k+1)，

γk+1=‖s(k+1)‖22，

βk=γk+1/γk，

p(k+1)=s(k+1)+βkp(k).

終止條件是‖s(k)‖/‖s(0)‖<10-8，A*表示A的共軛轉置.

6 數(shù)值實驗

在Merton模型^[12]下做數(shù)值實驗，當μJ=0時，歐式看漲期權有解:

ω(t，s)=∑

SymboleB@ m=0e-λ(1+η)τ(λ(1+η)τ)mm!VBS(τ，s，K，rm，σm)，(31)

其中，τ=T-t，η=eσ2J-1，σ2m=σ2+mσ2J/τ，

rm=r-λη+mlog (1+η)/τ，VBS表示歐式看漲期權的價格.

用Matlab編程進行數(shù)值試驗.在所有的實驗中，式(22)的迭代停止時刻由前后兩個迭代矩陣之間的差的l

SymboleB@ -范數(shù)決定，即當‖Vl+1-Vl‖

SymboleB@ <ε時停止，這里取ε=10-8.

在Merton模型中用FD和BDF2對空間與時間進行差分，并用三對角分裂法處理Toeplitz矩陣，到期時刻T=1，截斷點x*=4，r=0，波動率σ=0.2，跳躍方差σJ=0.5，跳躍強度λ=0.1，協(xié)定價格K=1，xK=log (K).結果為:

由表1和表2，隨著差分節(jié)點數(shù)的增加，計算的誤差越來越小，從空間差分節(jié)點數(shù)129開始，差分節(jié)點數(shù)每增加一倍，雅克比正則分裂迭代算法的計算誤差就下降一個數(shù)量級，對于預條件共軛梯度法，當差分節(jié)點數(shù)從65增加到129時，計算誤差下降了兩個數(shù)量級，

而在節(jié)點數(shù)從129增加到1 025的過程中，誤差都是在同一數(shù)量級內(nèi)減少，這說明預條件共軛梯度法的收斂速度很快.從計算的精確度來說，雅克比正則分裂法和預條件共軛梯度法相差不大，但是從計算的速度來看，后者要比前者快.

本文討論了當市場有跳躍時歐式期權定價的數(shù)值計算方法，期權的價格是一個PIDE方程的解，本文用差分法對這個方程進行離散，得到一個Toeplitz矩陣系統(tǒng)，本文用兩種方法來處理這個系統(tǒng)，由表1和表2可以看出，二者在計算精度上差別不大，預條件共軛梯度法比雅可比正則分裂法的迭代結果誤差更小些，而且迭代過程中的迭代次數(shù)更少，分析這些差別的原因，是由于預條件共軛梯度法對系數(shù)矩陣進行了處理，使系數(shù)矩陣的條件數(shù)減小，因而加快了迭代的收斂速度.

參考文獻

[1] BLACK F，SCHOLES M.The price of options and corporate liabilities[J].Journal of Political Economy，1973，81(3)，637-654.

[2] ANITA Mayo. Methods for the rapid solution of the pricing PIDE in exponential and merton models[J]. Journal of Computational and Applied Mathematics:2008，22(34):128-143.

[3] CONT R，VOLTCHKOVA E.A finite difference scheme for option pricing in jumpdiffusion and exponential levy model[J].SIAM J 2005，43(67):1596-1626.

[4] 楊向群，吳巒東.帶跳的冪型支付歐式期權定價[J].廣西師范大學學報:自然科學版，2007，25(34)，56-58.

[5] STANG G.A proposal for toeplitz calculations[J].Stud Appl Math，1986，74(39):171-176.

[6] CHAN T.An optimal circulant preconditioner for Toeplitz systems[J].SIAM，J，Sci，Stat，Comput，1988，9(13):766-771.

[7] BRIANT M，NATALINI R，RUSSO G.Implicitexplicit numerical schemes for jumpdiffusion process[J].Technical Report，2004，38(37):35-45.

[8] YOUNG D M.Iterative solution of large systems[J].New York:Academic，1971，5(23):25-35.

[9] ARIEL Almendral，CORNELIS W.Oosterlee. Numercial valuation of options with jumps in the underlying[J]. Applied Numercial Mathematics，2005，53(29):1-18.

[10]CHAN R，NAGY J，PLEMMONS R.Circulant preconditioned:toeplitz least squares iterations[J]. SIAM J Matrix Appl，1994，15(8):80-97.

[11]BRIANI M，Numerical methods for option pricing in jumpdiffusion markets[D].Universita Degli Studi Di Roma “La Sapienza”Dottor to Di Ricerca in Miatematica Per Le Applicazioni Economiche e Finanziarie，2003.

[12]ANDERSEN L，ANDREASEN J.Jumpdiffusion process:volatility smile fitting and numerical methods for option pricing[J].Rew.Derivatives Res，2000，4(17):231-262.

Numerical Solution of Assets Pricing Equation under Jumpdiffusion Model



ZHANG Hongyan， LI Qiang， ZHANG Zhi

(School of Mathematical and Calculating Technology，Central South University，ChangSha，Hunan 410083，China)

Abstract The paper assume that the price process of the assets is a jumpdiffusion process， then， the value ofEuropean optaon satisfies a general partial integrodifferential equation(PIDE) under this assumption.The equation was discretized by difference formula.The result was obtainedby two iterative methods:Jacobi regular splitting method and preconditioned conjugate gradient method.

Keywordsjumpdiffusion model; finite differences; FFT algorithm; European call option