# Option Pricing Based on Modified Advection-Dispersion Equation: Stochastic Representation and Applications

Authors: Longjin Lv; Luna Wang
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, we first investigate the stochastic representation of the modified advection-dispersion equation, which is proved to be a subordinated stochastic process. Taking advantage of this result, we get the analytical solution and mean square displacement for the equation. Then, applying the subordinated Brownian motion into the option pricing problem, we obtain the closed-form pricing formula for the European option, when the underlying of the option contract is supposed to be driven by the subordinated geometric Brownian motion. At last, we compare the obtained option pricing models with the classical Black–Scholes ones.


---

## Body

## 1. Introduction

Recently, the diffusion equations that generalize the usual one have received considerable attention due to the broadness of their physical applications, in particular, to the anomalous diffusion [1]. For instance, surface growth, transport of fluid in porous media, two-dimensional rotating flow, laser cooling, diffusion on fractals, or even in multidisciplinary areas such as in analyzing the behavior of CTAM micelles dissolved in salted water or econophysics [2–4]. In general, anomalous diffusion may be classified by employing the second moment x2. When x2∝tγ, the value γ>1 characterizes a superdiffusive process, γ<1 a subdiffusive one, and γ=1 a normal diffusion. In order to describe this phenomenon clearly, one needs to introduce the fractional Fokker–Planck equation (FFPE). The FFPE for anomalous diffusion are obtained as particular cases of the Kolmogorov’s equation [5]. In [6], Dubkov used the functional analysis approach to derive the FFPE directly from Langevin equation with symmetric α-stable Lévy noise.On one hand, several methods were introduced to get the solution of the FFPE in [7]. However, the limitation of such an approach is that it does not allow one to construct and analyze sample paths of the underlying stochastic process. In [8], the authors introduced a simple and efficient method for computer simulation of sample paths of anomalous diffusion process described by the FFPE. It reveals that subdiffusion is actually a combination of two independent mechanisms: the first mechanism is the standard diffusion represented by some Itô process Xτ, and the second mechanism introduces the trapping events and is represented by the so-called inverse α-stable subordinator Sαt. The subordinated process XSαt combines both mechanisms and gives the subdiffusive dynamics. The inverse α-stable subordinator is defined in the following way:(1)Sαt=infτ>0:Uατ>t,where Uατ is a strictly increasing α-stable Lévy process. Here, the α-stability means the Uατ has independent increment and Uατ+Δτ−Uατ=UαΔτ=Δτ1/α, in law. And, its Laplace transform is given by Ee−kUατ=e−τkα,0<α<1. Since Uατ is a pure-jump process, then, for every jump of Uατ, there is a corresponding flat period of its inverse Sαt. These heavy-tailed flat periods of Sαt represent long waiting times in which the subdiffusive particle gets immobilized in the trap. As $\α approaches 1, the random time Sαt. reduces to the normal time t. In other words, the probability density function (PDF) of the subordinated process XSαt is the solution of the FFPE [9].On the contrary, in financial market, the classical Black–Scholes model is based on the diffusion process called geometric Brownian motion [10, 11]:(2)dSt=μStdt+σStdBt,where St is the stock’s price and μ,σ are the expected average growth for the price and the expected noise intensity (the volatility) in the market dynamics, respectively. dS/S is usually called price return. The model is a geometric random walk with drift μ and diffusion σ. Bτ is the standard Brownian motion.However, the empirical studies show that many characteristic properties of markets cannot be captured by this model [12]. For example, (1) the leptokurtic feature: in other words, the historical data show that the return distribution has a higher peak and heavier tails than that of the normal one [13]. (2) Implied volatility smiled: if the Black–Scholes model is completely correct, then the implied volatility should be constant. In reality, it is widely recognized that the implied volatility curve resembles a smile, meaning it is a convex curve of the strike price. This empirical phenomenon is called volatility smile in option markets [14]. (3) Long memory: more recently, a number of authors have noted the apparent long memory property of powers of absolute returns and also of the volatility process of high frequency asset returns data. This has led to the formulation of long-memory time-dependent conditional heteroscedastic processes such as GARCH and also of corresponding long-memory stochastic volatility processes [15].In order to capture these financial characteristics, Bonanno [16, 17] applied a modified nonlinear Heston model to analyze the dynamics of stock prices with stochastic volatility. Magdziarz applied the time-changed Brownian motion into the option pricing problem [18–20]. The correlated continuous time random walk is also employed to describe the stock price, and the option pricing formula is obtained in [21]. In [22, 23], Aguilar showed that the price of an European call option, whose underlying asset price is driven by the space-time fractional diffusion, can be expressed in terms of rapidly convergent double-series.Based on the previous scholars’ research studies, in this paper, we first use the stochastic representation method to study the modified advection-dispersion equation and find the subordinated Brownian motion model which can capture the financial characteristics well. Then, we apply the subordinated geometric Brownian motion into the option pricing problem and deduce the option pricing formula.This paper is organized as follows. In Section2, we derive a subordinated process and prove that the PDF of this process is rightly the solution of the modified advection-dispersion equation, where the parent process is a classical diffusion process and the subordinator is the inverse of a Lévy motion with drift. Two special cases are also introduced to help understand the process clearly. The mean-squared moment of the equation is also discussed in this section. Taking advantage of the stochastic representation method, in Section 3, we apply this model to option pricing problem. In Section 4, we present our conclusions.


## 2. Stochastic Representation

The main aim of this section is to find the stochastic representation for the following differential equation:(3)m∂Cx,t∂t+β∂Cx,t∂t∗ft=−∂∂xvCx,t+D∂2∂x2Cx,t,which was originally proposed by Schumer et al. [24] to describe the mobile/immobile (MIM) solute transport phenomenon. The original MIM model builds upon the fact that, in many natural porous media, tracers may temporarily be stopped at immobile sites, which may be disseminated in the solid matrix. The model was indeed developed by hydrologists, often facing situations where only the “mobile” phase is accessible to field measurements. The simplest form of the model assumes first-order kinetics for exchanges between mobile and immobile phases. Solving the kinetic equation ruling the immobile concentration allows the determination of the total density of tracer in the form of a single equation. The equation appears as a modified advection-dispersion equation (ADE), involving a convolution with ft, besides the regular time derivative on the left-hand side. The convolution in equation (3) leads to memory effects. Here, ft can be any continuous function, and Cx,t is the total density of tracer, including mobile and immobile phases.LetDτ be an increasing Lévy motion with Laplace transform Ee−kDτ=e−τψk. Aτ is the increasing Lévy motion with drift, Aτ=mτ+Dτ m>0. Let Et be the inverse time of Aτ, defined as(4)Et=infτ>0:Aτ>t,f^k denotes the Laplace transform of the function ft. Then, we have the following theorem.Theorem 1.
The PDF of the subordinated processYt=XEt is the stochastic representation of the modified advection-dispersion equation (equation (3)) for ψk=βkf^k, where the parent process Xτ is defined as(5)dXτ=vdτ+2DdBτ.
HereEt is independent of Bτ.Proof.
Since the Laplace transform ofDτ is given as Ee−kDτ=e−τψk D0=0, then Ee−kAτ=e−τk+ψk with Aτ≥mτ for all τ. Firstly, we establish the relation between the PDF gτ,t of Et and the PDF ut,τ of Aτ. From the definition of Et, we have PEt<τ=PAτ>t [25, 26]; therefore,(6)gτ,t=∂∂τ∫t∞uy,τdy=−∂∂τ∫mτtuy,τdy.
So, the Laplace transform ofgτ,t can be expressed as(7)g^τ,k=∫0∞e−ktgτ,tdt=−∂∂τ∫mτ∞uy,τ∫y∞e−ktdtdy=−1k∂∂τ∫mτ∞e−kyuy,τdy=−1k∂∂τe−τmk+ψk=m+ψkke−τmk+ψk,where we have used the property that ut,τ=0 for t<mτ. Using the total probability formula and the independence between Xτ and Et, we get the PDF px,t of Yt, given by(8)px,t=∫0∞fx,τgτ,tdτ.where fx,τ is the PDF of the parent process Xτ. So, the Laplace transform of the above equation yields(9)p^x,k=∫0∞fx,τg^τ,kdτ=m+ψkk∫0∞fx,τe−τmk+ψkdτ.
Thus, the following relation betweenfx,τ and px,τ holds(10)p^x,k=m+ψkkf^x,mk+ψk.
Since the processXτ is given by the Itô stochastic differential equation (3), its PDF fx,τ obeys the classical advection-dispersion equation [27]:(11)∂fx,τ∂τ=−∂∂xv+D∂2∂x2fx,τ.
The Laplace transform of the above equation with respect toτ yields(12)kf^x,k−fx,0=−∂∂xv+D∂2∂x2f^x,k.
By changing the variablek to mk+ψk and using the relation (10), the above equation yields(13)m+ψkkkp^x,k−px,0=−∂∂xv+D∂2∂x2p^x,k.
Comparing the above equation with the Laplace transform of equation (3), we get Cx,t=px,t if ψk=βkf^k, that ends our proof.
The superiority of the stochastic representation approach to the fractional differential equation is that it not only provides a way to get the solution of the corresponding equations, but also helps us to understand the physical process by providing a description of the dynamical system governed by the fractional differential equation. So, we can use Monte Carlo method to simulate the stochastic process, and then we can see how the particle moves from the figures. In other words, ifft is given, we can calculate ψk by using equation (8), and then inverting the Laplace transform in g^τ,k, the expression of gτ,t can be obtained. Substituting gτ,t into equation (8) leads to the analytical expression of the solution of equation (3). In order to describe these clearly, we present the following two cases.Case 1.
f
t
=
δ
t.
We can getf^k=1 and ψk=βk. So, from equation (7), we have g^τ,k=m+βe−τ1+βk, and then Et=t/1+β. Applying inverse Laplace transform on g^τ,k yields(14)gτ,t=δtm+β−τ.
Substituting it into equation (10), we have(15)Cx,t=fx,t/m+β.
Taking advantage of the relation betweenpx,t and fx,τ, we can get the following evaluation formula for the mean square displacement:(16)X2t=∫0∞x2px,tdx=∫​0∞gτ,t∫0∞x2fx,τdxdτ=2D∫0∞τgτ,tdτ.
Therefore,(17)X2t∝2D1+βt.
So, ifft=δt, the process is a normal diffusion.Case 2.
f
t
=
t
−
α
/
Γ
1
−
α.
For this case, equation (3) yields(18)m∂Cx,t∂t+β∂αCx,t∂tα=−∂∂xvCx,t+D∂2∂x2Cx,t,where the operator ∂αCx,t/∂tα is the Caputo fractional derivative of order α [28].
Following the same procedure shown in Case1, we can get Ee−kDτ=e−τkα, so Dτ=β1/αUατ is a strictly increasing α-stable Lévy motion. As we know, the PDF of inverse time α-stable Lévy motion Sαt can be expressed in the form of a Fox function, i.e., 1/τH1110τ/tα1,α1,1 (see [29]). By using the relation between Et and Dτ, after some calculations, we can get(19)gτ,t=1τ+αmt−mταH1110β1/ατt−mτα1,α1,1,where Hpqmnxa1,A1,…,ap,Apb1,B1,…,bq,Bq is the Fox function [30]. The solution can be simulated by Monte Carlo method (see Figure 1).
In order to calculate the mean square displacement, we can first compute its Laplace transform. By using the Laplace transform ofgτ,t equation (7), and then inverting the Laplace transform, we have(20)X2t∝2DtmE1−α,2−βt1−αm,where Eα,βx=∑k=0∞xk/Γαk+β is the Mittag-Leffler function [31]. Here, we have used the equality:(21)∫0∞e−pttαk+β−1Eα,βk±atαdt=k!pα−βpα∓ak+1.
From equation (20), we can know if m≠0, the model resembles a normal diffusion for t⟶0+, and since 0<α<1, the model resembles a subdiffusion for t⟶∞. This result can be obtained in another way, to see this, note that Uατ=τ1/αUα1 in distribution and τ1/α grows faster than mτ for 0<α<1, so the stable subordinator dominates as τ⟶∞ and the drift dominates as τ⟶0+. From the definition of Uατ and Sαt, we also should note that τ⟶0+ and τ⟶∞ are equivalent to t⟶0+ and t⟶∞, respectively. By the way, as we know, the generalized Einstein relation is defined by X2tF=0=const×XtF≠0 [32]. From equations (8) and (16), we should notice that if the Einstein relation for the parent process Xτ holds, then it also holds for the subordinated process.Figure 1
The probability density functionCx,t with different m and α. Here, v=0, D=0.5, β=1, and t=1.


## 3. Application to Option Pricing Problem

In this section, we aim to evaluate the price of an European option when the underlying of the option contract is supposed to be driven by a subordinated geometric Brownian motionSt=XEt , where the parent process Xτ is given as(22)dXτ=μXτdτ+σXτdBτ,and the subordinator Et is the inverse time of Aτ. Since Uατ is 1/α-similar, its inverse time Sαt is α-similar. The self-similarity is a quite important property in financial market. In order to employ this feature, we choose Aτ=τ+Uατ. Here, Et is also independent of Bτ.In fact, from Section2, we can get the PDF px,t of log-price of this underlying should satisfy the following FFPE:(23)∂∂tpx,t+∂α∂tαpx,t=−μ−σ22∂∂x+σ22∂2∂x2px,t,with initial condition px,0=δx. Here, the fractional derivative is also the Caputo type. In [33], Ren et al. have proved that the detailed structure of the solution px,t depends generally on the special shape of the underlying geometry. However, the interesting part of px,t has the asymptotic behavior logpx,t∼Cξu, where ξ≡x/tα/2≫1; here, u=1/1−α/2, i.e., the solution px,t for (22) follows a stretched Gaussian distribution, which has heavy tail and high peak, in contrast to a normal distribution. Heavy tail and high peak phenomenon is common in the financial historical data [34].From equation (23), we know the PDF px,t follows the fractional diffusion equation, where the Caputo fractional operator is given as(24)Dtα0ft=1Γ1−α∫0tt−s−αdfsdsds,0<α<1.The definition shows the Caputo fractional operator is a convolution integral, which impliesBSαt is a long-memory process. However, there is substantial evidence that long-memory processes describe rather well financial data such as forward premiums, interest rate differentials, and inflation rates. Perhaps, the most dramatic empirical success of long-memory processes has been in recent work on modeling the volatility of asset prices and power transformations of returns [15].In summary, the stochastic differential equationXEt presents many important financial properties, such as self-similarity, leptokurtic feature, and long memory. So, we assert this model characterizes the price of underlying well.Letrt be the return of the underlying from time 0 to time t. Then, we have(25)rt=μEt+σBEt,in law.From the discussion, in the Section 2, we know the PDF of the returns is leptokurtosis and fat-tailed, given as(26)qx,t=∫0∞12πσ2τex−μτ2/σ2τgτ,tdτ,where gτ,t is showed in equation (19) with m=1 and β=1. By using the stochastic representation, we can also get the statistical signatures of the stock price returns, such as the mean value, variance, and approximated solution (see equation (20) and Figure 1).Now, let us pass on our main problem. LetCt,St be the value at time t of an European option on the above underlying with expiration data T and exercise price K and the boundary condition CT,ST=ST−K+.Theorem 2.
Assume that the underlying is driven by the subordinated geometric Brownian motionSt=XEt, then the European call option Ct,St can be given as(27)Ct,St=StNd1−Ke−rTNd2,where the function Nx is the cumulative probability distribution function for standard normal distribution, and(28)d1=lnSt/K+rT−t+σ2/2T−t+Tα−tα/Γα+1σT−t+Tα−tα/Γα+1,d2=d1−σT−t+Tα−tαΓα+1.Proof.
We assume that the riskless bond price dynamics satisfies(29)ΔΠt=rΠtΔt,where r is the riskless interest rate. From Taylor’s formula, we have(30)ΔCt,St=∂C∂tΔt+∂C∂SStΔSS+12∂2C∂S2ΔSS2+oΔSS2.
Since theBτ is independent of Et, we get(31)EΔSS2=EμΔEt+σΔBEt2=σ2EΔEt+oΔt=1+tα−1Γασ2Δt+oΔt,where we have used the property of Sαt, i.e., Sαt is a α-similar process and ESαnt=tnαΓn+1/Γnα+1 [35]. Noting that it is possible to hedge the risk of a portfolio, like the Black–Scholes case [10, 11], we assume that the value of the riskless bond can be exactly replicated by following “self-financing” investment strategy, involving the underlying and an option on this underlying(32)Πt=Ct,St−∂C∂SSt,
So,(33)ΔΠt=ΔCt,St−∂C∂SΔSt.by substituting equations (30)–(33) into (29), we have(34)∂C∂t+rS∂C∂S+12σ˜2S2∂2C∂S2=rC,where σ˜2t=1+tα−1/Γασ2. So, the solution of (34) with boundary condition CT,ST=ST−K+ can be given as(35)Ct,St=StNd1−Ke−rTNd2,where the function Nx is the cumulative probability distribution function for standard normal distribution, and(36)d1=lnStK+rT−t+12∫tTσ˜2sds,d2=d1−∫tTσ˜2sds.
Substitutingσ˜2t into the above equation ends the proof.
Noting that this option pricing formula is quite similar to the one obtained in [36], where the price of the underlying is supposed to be driven by a geometric fractional Brownian motion. We should also notice that all the results obtained here are consistent with that obtained by the classical Black–Scholes formula when α=0.
By comparing this formula with the classical Black–Scholes pricing formula for differentα in Figure 2, we find the call option price increases with the increase of α, which is due to the heavy-tailed assets return. Under this situation, the risk is redistributed, which results in an adjustment of the option price. In the real world, the option is an insurance product, which is used to hedge the risk of the underlying assets. The return of the underlying assets is supposed to be heavy tailed, which means the underlying assets is more risky. The risk also increases with the increase of α, that leads to the smaller the parameter α is, the higher the option price will be.
The differences between this option pricing formula and the classical Black–Scholes model for different expiration times are presented in Figure3. One can observe that the price differences become smaller as the maturity time increases. That is caused by the fact that the α-stable term dominates in the short time and the drift term dominates in the long time. This behavior can be observed in situations when we face some kind of unexpected or sudden change of regime, such as a black day on the market, the bankruptcy of a company trading on the market, and a natural disaster [37]. In short time, people are willing to spend more money to protect their assets.
The empirical analysis on SPDR S&P500 ETF call option (SPY) is introduced to verify our results. The option written on this index is actively traded European-style contracts. We choose the data April 5, 2019, as the quotation date and April 12, 2019, as the expiration date. At that moment, the index of S&P500 ETF is 288.57, and the risk-free rate is 2.5%, (source:yahoo.com). We use the Parkinson’s formula [38] to estimate the parameter σ, given as σ=1/4nln2∑i=0nSih−Sil21/2, where Sih and Sil are the highest and lowest prices of the same day, respectively. We get the parameter σ=0.153985. In Table 1, we provide observable market bid (offered) prices for SPDR S&P500 ETF call options with different exercise (strike) prices. In Figure 4, we plot the prices obtained in Table 1. We find the option price obtained by our model is higher than the classical Black–Scholes price and closer to the real price, which verifies our results.Figure 2
The call option price with differentα. Here, r=0.05, σ=0.8, St=10, T=2, and t=1.Figure 3
The call option price calibrated by theα-BS model and BS model for different expiration times T. Here, r=0.05, σ=0.8, and St=10.Table 1
Price for SPDR S&P500 ETF call option.
K
Real price
α-BS (α=0.1)
BS price
286
3.89
3.4576
3.2458
287
3.14
2.7863
2.5489
288
2.44
2.1955
1.9432
289
1.9
1.6888
1.4346
290
1.45
1.2664
1.0232
291
1.02
0.9245
0.7035
292
0.69
0.6562
0.4655
293
0.48
0.4524
0.2959
294
0.3
0.3027
0.1804
295
0.2
0.1963
0.1054
296
0.11
0.1234
0.059Figure 4
Comparison between the real call option prices with the one obtained by different models under different strike pricesK. Here, BS price stands for the price obtained by the classical Black–Scholes model, and α-BS stands for the price obtained by formula (26) with α=0.1.


## 4. Conclusions

In this paper, we employ a modified advection-dispersion equation, including the first-order derivative with respect to time and its convolution integral with a function on the left-hand side. The stochastic representation of this equation is proved to be a subordinated process. Its subordinator is the inverse of a Lévy process, whose characteristic function is dependent on the function presented in the convolution. The solution, mean square displacement and Einstein relation of this equation are all discussed for two special cases. We find they are greatly dependent on those of the corresponding standard advection-dispersion equation. Then, we apply this model to option pricing problem. The evaluation formula of a European option is obtained when the underlying of the option contract is supposed to be driven by a subordinated geometric Brownian motion. The differences between the obtained results and the classical ones are also given. Finally, we expect that the results obtained here may be useful to the discussion of the anomalous diffusion systems and complex financial market.



---

*Source: 7168571-2020-03-12.xml*