# Solving Parametric Volterra Integral Equation from Distributed-Order Rough Heston Model and Option Pricing

Authors: Zhengguang Shi
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The rough Heston model has recently attracted the attention of many financial practitioners and researchers because it maintains the basic structure of the classic Heston model and has an advantage in describing the microstructure foundation of the market. In this paper, we study the distributed-order rough Heston model with an exponential tempered factor, and from the characteristic function of log-price in this model, we obtain a nonlinear parametric Volterra integral equation. Finally, the Fourier-cosine methods are combined with the Adams methods to price the option.


---

## Body

## 1. Introduction

It is well known that the fractional Brownian motion of Hurst parameterH∈0,1 is a class of zero-mean Gaussian processes BH=BtH,t≥0, with covariance functions given by(1)Covt,s=12t2H+s2H−t−s2H.We can see that ifH=1/2, BH is the ordinary Brownian motion, H<1/2 yields negatively correlated increments and rather erratic paths suitable for modeling systems with fast changing states, and H>1/2, the increments are positively correlated, and the process exhibits long-range dependence. In nature, the dynamics of a system are often observed over time. For example, during market shocks, one can deduce unusual movements in the volatility of a stock price. It seems like that fractional Brownian motion with H>1/2 seems to be a good choice for modeling the volatility, but some researches show that different H values for different time intervals are necessary, for example, Vehel and Riedi [1]; Chronopoulou and Viens [2]; Ryvkina [3] and the reference therein.Meerschaert and Sabzikar [4] provide that a new stochastic process called tempered fractional Brownian motion Bα,λH with tempering parameter λ≥0 is introduced by representation as(2)Bα,λH=∫−∞te−λt−st−sH−1/2Bds−∫−∞0e−λ−s−sH−1/2Bds,where Bdx is independently scattered Gaussian random measure on ℝ. Tempered fractional Brownian motion can be defined by modifying the moving average representation of a fractional Brownian motion with the inclusion of an exponential tempering factor to the power-law kernel. In contrast to the fractional Brownian motion, tempered fractional Brownian motion has some nice features, such as scaling property, stationary increments, moving average representation. Also, it has an advantage in modeling realistic data when the tempering parameter modifies the power law kernel which can be chosen to be the desired degree of accuracy over a finite interval. The kernel e−λt−sH−1/2 in (2) plays a central role in the rough dynamic of the tempered fractional Brownian. In particular, one can show that the process(3)∫0te−λt−st−sH−1/2dBs,where B is standard Brownian motion. Then, Meerschaert and Sabzikar [5] develop the theory of stochastic integration for tempered fractional Brownian motion. Chen, Wang, and Deng [6] discuss the tempered fractional Brownian motion, its ergodicity, and the derivation of the corresponding Fokker-Planck equation. Sabzikar and Surgailis [7] study a closely related but different tempered process called tempered fractional Brownian/stable motion of the second kind. Lim and Eab [8] study the properties of tempered fractional Brownian motion via fractional Ornstein–Uhlenbeck processes.The classical Heston model is one of the most flexible and widely used parametric asset price models, where the volatility process follows the stochastic differential equations driven by Brownian motion. However, empirical studies show that considering a fractional Brownian motion with a small Hurst parameter also enables us to obtain remarkable fits for the whole volatility surface. In particular, contrary to most stochastic volatility models, rough volatility models generate an exploding term structure for the at-the-money skew when maturity goes to zero, which is very commonly observed in practice (e.g., Bayer, Friz, and Gatheral [9]; Jeng and Kilicman [10] and references therein).El Euch and Rosenbaum [11] used an original link between nearly unstable Hawkes processes and fractional volatility models, and then, the nonlinear fractional Riccati equation for the characteristic function of the asset log-price is derived. Moreover, they use the Adams scheme to solve the fractional Riccati equations. Recently, Gatheral and Radoicic [12] use multipoint Padé approximants on the asymptotical solutions (t⟶0 and t⟶∞) of the fractional Riccati equation and apply it to the rough Heston model. However, the multipoint Padé method performs extremely well with H⟶0 but not H⟶0.5. Jeng and Kilicman [13] used the Adams method and the multipoint Padé approximation method for solving the fractional Riccati equation. In addition, they investigate the performance of the multipoint Padé approximation method for the small u values in Dαhu−i/2,t as it plays a huge role in the computation for the option prices. Jeng and Kilicman [14] show that under the Laplace-Adomian decomposition method, the infinite series expansion of the fractional Riccati equation’s solution corresponds to the existing expansion method from previous work for at least up to the fifth order, and then, they show that the fourth-order Padé approximants can be used to construct an extremely accurate global approximation to the fractional Riccati equation in an unexpected way.Our main goal is to design the distributed-order rough Heston model with exponential tempered factor as in El Euch and Rosenbaum [11] by combining the variable Hurst exponent H as well as the advantages of tempered fractional Brownian motion; namely, we replace kernel t−sα−1 in El Euch and Rosenbaum [11] by e−λt−st−sα−1 where α\≔H+1/2 is a random variable. Then, from the characteristic function of log-price in this model, we obtain a nonlinear parametric Volterra integral equation in the expectation sense; that is,(4)ht=∫0tEt−sα−1ΓαLt,s,hsds,where α follows a uniform distribution on 1/2,1, and Lt,s,hs is given in the next section. Nevertheless, to the extent of our knowledge, published papers on the Adams method for solving Volterra integral (4) on graded meshes are very rare. Hence, we will use the Adams method (refer to Liu, Roberts, and Yan [15]) to compute the Volterra integral equation on graded meshes and use the Fourier-cosine methods (refer to Fang and Oosterlee [16] in details) to price the European option under the rough model.We organize the rest of the paper as follows: in Section2, we introduce the distributed-order rough Heston model with an exponential tempered factor, and then, we derive the Volterra integral equation from the characteristic function of the asset log-price. In Section 3, we investigate the Adams method on the graded mesh to solve the Volterra integral equations. In Section 4, we use the Fourier-cosine methods to compute the option prices. Numerical examples are given in Section 5, and conclusions in Section 6.


## 2. The Rough Heston Model

The distributed-order rough Heston model with exponential tempered factor is as follows:(5)dSt=StVtdWt,Vt=V0+γΓα∫0te−λt−st−sα−1θ−Vsds+σΓα∫0te−λt−st−sα−1VsdBs,where the parameters γ,θ,σ,S0,V0 are positive constants, Wt;t≥0 and Bt;t≥0 are two correlated Brownian motions with dWt,dBt=ρdt, ρ∈−1,1, and α follows a uniform distribution on 1/2,1, that is, α∼U1/2,1. Then, by means of Theorem 3.2 [17], we deduce the characteristic function of (5) as(6)Θu;lnS,V≔EexpulnST=expulnS+γθ∫0Thsds+∫0T12u2−u+uρσ−γhs+σ22h2sdsV,where T<∞, u∈ℂ, h∈L20,T,ℂ, S and V represent the values of stock and volatility at t=0, respectively, and ht satisfies a Volterra integral equation as follows:(7)ht=∫0tEt−sα−1ΓαLt,s,sds,where Lt,s,hs=e−λt−s1/2u2−u+uρσ−γhs+σ2/2h2s and h0=0.Remark 1.
The proofs of (6) and (7) are similar to Theorem 3.2 in [17], here noting that we define a convolution in the expectation sense.


## 3. Adams Method for Volterra Integral Equation

Now, we use the two-step fractional Adams scheme on graded meshes to solve the Volterra integral (7). Graded meshes on 0,T are defined as(8)tk=TkNp,k=0,1,…,N,withp≥1.On the graded meshes, we estimate(9)htk+1=∫0tk+1Etk+1−sα−1ΓαLtk+1,s,hsds,≈∫0tk+1Etk+1−sα−1ΓαL∧tk+1,s,hsds,k=0,1,…,N−1,where(10)L∧tk+1,t,h1=tj+1−ttj+1−tjLtk+1,tj,htj+t−tjtj+1−tjLtk+1,tj+1,htj+1.It then leads to(11)htk+1≈∑j=0kaj,k+1Ltk+1,tj,htj+ak+1,k+1Ltk+1,tk+1,htk+1,where(12)aj,k+1=ETαN−pαΓ2+αEk+1pα1+α+k+1p−11+α−k+1p1+α,j=0,Ek+1p−j−1p1+α−k+1p−jp1+αjp−j−1p+k+1p−j+1p1+α−k+1p−jp1+αj+1p−jp,1≤j≤k,Ek+1p−kpαj=k+1.Moreover, using product rectangle rules, we derive that(13)htk+1=∫0tk+1Etk+1−sα−1ΓαLtk+1,s,hsds≈∑j=0kbj,k+1Ltk+1,tj,htj≔hPtk+1,where(14)bj,k+1=ETαN−pαΓ1+αk+1p−jpα−k+1p−j+1pα.Combining (6) and (8), we give the two-step fractional Adams scheme, which is an explicit scheme, as(15)h∧tk+1=∑j=0kaj,k+1Ltk+1,tj,h∧tj+ak+1,k+1Ltk+1,tk+1,h∧Ptk+1,h∧0=0,where h∧ is the approximation of h, the weights aj,k+1 are defined in (7), and the correction term h∧Ptk+1 denotes the hPtk+1 in (13) with replacement of h by h∧.Remark 2.
Noting thatα∼U1/2,1, for (7) and (9), we will use the Gauss–Legendre numerical integration formula to solve them; that is,(16)Eℱα=12∫1/21ℱxdx=12∑ℓ=1n−1tℓ+1−tℓ2∫−11ℱtℓ+1−tℓ2x+tℓ+1+tℓ2dx≈12∑ℓ=1n−1tℓ+1−tℓ2∑q=0mGqℱtℓ+1−tℓ2xq+tℓ+1+tℓ2,where t1=1/2,tn=1, Gqq=0,1,…,m is the weight and xq is the Gaussian point.


## 4. Fourier-Cosine Methods for Option Pricing

In this section, we investigate the Fourier-cosine methods for the prices of the European put option for the distributed-order rough Heston model with exponential tempered factor.Firstly, using the Richardson extrapolation technique of trapezoidal formula and scheme (10) for approximating h, we derive the numerical formula for the characteristic function (4) as follows:(17)Θu;lnsk,V=expulnsk+γθ∫0Thsds+∫0TgsdsV≈expul0skn+4I2−I13≕Θ∧u;lnsk,V,where K is strike price, gt=1/2u2−u+uρσ−γht+σ2/σ2h2t,(18)I1=12∑k=0N/2−1γθh∧t2k+h∧t2k+2+gh∧t2k+gh∧t2k+2Vt2k+2−t2k,I2=12∑k=0N−1γθh∧tk+h∧tk+1+gh∧tk+gh∧tk+1Vtk+1−tk,and N is an even number.Under risk-neutral measureℚ, the pricing formula of rough Heston model with exponential tempered factor is(19)vt,S,V=Et,S,VℚK−ST+,where Et,S,V is a conditional expectation operator given that St=S,Vt=V. Let x≔lns/k and y≔lnST/k. According to Fang and Oosterlee [16] (interesting readers should refer to this literature in detail), we have(20)vt,S,V=Et,x,VℚK1−ey+=∫ℝK1−ey+fy|x,Vdy≔v˜t,x,V,where fy/x,V is the probability density of y given x,V.Then, using the cosine expansion forfy/x,V and (11), we can obtain that numerical formula of the European put option pricing (12) for the rough Heston model with exponential tempered factor as follows:(21)v˜t,x,V≈2η2−η1∑′k=0M−1ReΘ∧ikπη2−η1;x,Vexp−ikπη1η2−η1Vk,where Re⋅ denotes taking the real part of the argument and M represents the number of truncation, and we can choose the boundary of the integral(22)η1,η2=c1−Lc2,c1+Lc2,where L is a positive constant which is large enough to guarantee η1<0<η2,(23)cn=1in∂nlnΘ∧iω;x,T∂iωn|iω=0,n=1,2,Vk=∫η1η2v˜T,x,Vcoskπy−η1η2−η1dy.Remark 3.
Since characteristic function (4) has no analytic solution, we calculate the cumulants cn numerically using finite differences. Let Ciω≔lnψiω be the characteristic exponent of the distributed-order rough Heston model (3), and then, the first second finite differences denoted by dnCiω are given by(24)d1Ciω=Ciω+h−Ciω−h2h,d2Ciω=Ciω+h−2Ciω+Ciω−hh2,where h is a small value.


## 5. Numerical Examples

In this section, we carry out numerical examples to verify the convergence of Adams scheme (10) for solving the Volterra integral equation (5) and convergence rates of Fourier-cosine methods for (13) under the distributed-order Heston model with exponential tempered factor. The codes for the examples are run in MATLAB R2017a on a PC with the configuration: Intel(R)Core(TM) i7-8550UCPU@1.80 GHz 1.99 GHz and 8.0 GB RAM.Example 1.
The aim of this example is to verify the two-step fractional Adams scheme (10) for solving the Volterra integral equation (5). The values of parameters are given as follows:(25)γ=0.3,ρ=−0.7,σ=0.3,θ=0.2,V=0.2,and the time horizon is taken as T=1. The complex number u is taken as u=i, and the graded mesh parameter is set by p=3.
From Table1, we can see that the numerical results are convergent and effective.Table 1
Convergence of the two-step fractional Adams scheme (10) for Example 1.
Nλ=0.04λ=0.06λ=0.08ErrorErrorError80———1003.35e − 43.36e − 43.37e − 41202.31e − 42.32e − 42.32e − 41401.69e − 41.69e − 41.70e − 41601.29e − 41.29e − 41.30e − 41801.02e − 41.02e−41.02e−4Example 2.
In this example, we verify the convergence of the Fourier-cosine methods (13) to calculate the prices of the European put options for the distributed-order rough Heston model (3). The strike price K=80,100,120, the initial stock price S=100, and the other parameters are the same as Example 1.
Table2 gives the values of options, the errors, the computational time, and we can see that the numerical results are convergent.
As shown in Figure1, under the same parameter λ, the Fourier-cosine methods present an exponential convergence for the different strike price K and vice versa.Table 2
Computational results for European put option pricing for Example2 with Adams schemes. The number of Fourier nodes is taken as M=40.
KNλ=0.3λ=0.5λ=0.7ValueErrorTime (sec.)ValueErrorTime (sec.)ValueErrorTime (sec.)80807.672 790—177.683 448—177.693 013—171007.672 1766.13e − 4267.682 8416.06e − 4267.924 1256.01e − 4281207.671 7833.94e − 4387.682 4513.90e − 4397.692 0243.87e − 4411407.671 5092.74e − 4527.682 1782.72e − 4527.691 7542.70e − 4541607.671 3082.01e − 4697.681 9782.00e − 4777.691 5551.99e − 4721807.671 1531.54e − 4987.681 8251.54e − 4947.691 4021.53e − 4921008017.333 798—1717.358 220—1817.380 226—1710017.332 7251.07e − 32717.357 1321.09e − 32617.379 1241.10e − 33012017.331 9957.31e − 44017.356 3927.40e − 43817.378 3767.49e − 44414017.331 4655.30e − 45417.355 8555.36e − 45417.377 8335.42e − 45316017.331 0634.02e − 47217.355 4494.07e − 47117.377 4234.11e − 47018017.330 7473.15e − 49017.355 1303.20e − 49117.377 1003.22e − 4901208030.583 570—1730.613 037—1730.639 661—1710030.582 4431.13e − 33030.611 8751.16e − 32630.638 4681.19e − 32712030.581 6487.95e − 44130.611 0588.17e − 44130.637 6338.36e − 44414030.581 0575.90e − 45530.610 4536.05e − 46030.637 0156.17e − 45816030.580 6024.56e − 47230.609 9884.66e − 47730.636 5414.74e − 47318030.580 2403.62e − 49930.609 6193.69e − 49430.636 1653.76e − 494Figure 1
Errors of Fourier-cosine methods for pricing European put options under the rough Heston model withλ=0.2 and λ=0.8, respectively, and the other parameters are the same as Example 1. We choose N=50 and M=10,15,…,45, and the reference values are obtained by the Fourier-cosine method with N=400 and M=115.


## 6. Conclusions

In this paper, we introduce the distributed-order rough Heston model with exponential tempered factor, and from the characteristic function of log-price in this model, we obtain a parametric Volterra integral equation in the expectation sense. Then, we use a two-step fractional Adams scheme to solve the Volterra integral equation on graded meshes. Moreover, the European options are calculated by the Fourier-cosine methods. Numerical examples are provided to verify the efficiency and accuracy of the computations. From Table2, we see that the program is not running very fast because it contains the expectation in the integral equation, and the expectation itself is a numerical integral, and the whole part is a double integral which requires an approximated double integral with a suitable numerical Volterra integral equation and gives an approximated solution. Here, we only chose the Adams method because we get high convergence with the Adams method on graded meshes for Volterra integral equation without expectation (refer to Shi, Ma, and Liang [17]), and compared with other methods, the Adams method is more simple and easy to implement. Next, we will investigate convergence rates of the two-step fractional Adams scheme for the authors of (10). In the future, we will try to find more efficient numerical methods to solve the Volterra integral equation (10) as well as using meshless methods (refer to Assari and Dehghan [18, 19], Assari and Cuomo [20], and references therein) to study the properties of the Volterra integral equation (10).



---

*Source: 1979003-2022-03-30.xml*