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

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. +en, 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.


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][3][4]. In general, anomalous diffusion may be classified by employing the second moment 〈x 2 〉. When 〈x 2 〉 ∝ t c , the value c > 1 characterizes a superdiffusive process, c < 1 a subdiffusive one, and c � 1 a normal diffusion. In order to describe this phenomenon clearly, one needs to introduce the fractional Fokker-Planck equation (FFPE). e 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).
e subordinated process X(S α (t)) combines both mechanisms and gives the subdiffusive dynamics. e inverse α-stable subordinator is defined in the following way: 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 is a pure-jump process, then, for every jump of U α (τ), there is a corresponding flat period of its inverse S α (t). ese 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 X(S α (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]: where S t 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. e 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. is 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. is has led to the formulation of longmemory time-dependent conditional heteroscedastic processes such as GARCH and also of corresponding longmemory 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][19][20]. e 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. en, we apply the subordinated geometric Brownian motion into the option pricing problem and deduce the option pricing formula.
is paper is organized as follows. In Section 2, we derive a subordinated process and prove that the PDF of this process is rightly the solution of the modified advectiondispersion 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. e 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.

Stochastic Representation
e main aim of this section is to find the stochastic representation for the following differential equation: which was originally proposed by Schumer et al. [24] to describe the mobile/immobile (MIM) solute transport phenomenon. e 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. e model was indeed developed by hydrologists, often facing situations where only the "mobile" phase is accessible to field measurements. e 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. e equation appears as a modified advection-dispersion equation (ADE), involving a convolution with f(t), besides the regular time derivative on the lefthand side. e convolution in equation (3) leads to memory effects. Here, f(t) can be any continuous function, and C(x, t) is the total density of tracer, including mobile and immobile phases.
Let D τ be an increasing Lévy motion with Laplace f(k) denotes the Laplace transform of the function f(t). en, we have the following theorem.

Theorem 1. e PDF of the subordinated process Y t � X(E t ) is the stochastic representation of the modified advection-dispersion equation (equation (3)) for ψ(k) � βkf(k), where the parent process X(τ) is defined as
Here E t is independent of B(τ).
Proof. Since the Laplace transform of D τ is given as [25,26]; therefore, So, the Laplace transform of g(τ, t) can be expressed as 2 Discrete Dynamics in Nature and Society where we have used the property that u(t, τ) � 0 for t < mτ.
Using the total probability formula and the independence between X(τ) and E t , we get the PDF p(x, t) of Y(t), given by where f(x, τ) is the PDF of the parent process X(τ). So, the Laplace transform of the above equation yields us, the following relation between f(x, τ) and p(x, τ) holds Since the process X(τ) is given by the Itô stochastic differential equation (3), its PDF f(x, τ) obeys the classical advection-dispersion equation [27]: e Laplace transform of the above equation with respect to τ yields By changing the variable k to (mk + ψ(k)) and using the relation (10), the above equation yields Comparing the above equation with the Laplace transform of equation (3), e 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, if f(t) 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. □ We can get f(k) � 1 and ψ(k) � βk. So, from equation (7), we have g(τ, k) � (m + β)e − τ(1+β)k , and then E t � t/(1 + β). Applying inverse Laplace transform on g(τ, k) yields Substituting it into equation (10), we have Taking advantage of the relation between p(x, t) and f(x, τ), we can get the following evaluation formula for the mean square displacement: erefore, So, if f(t) � δ(t), the process is a normal diffusion.
is 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 〈X 2 (t)〉 F�0 � const × 〈X(t)〉 F≠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.

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 motion S t � X(E t ) , where the parent process X(τ) is given as and the subordinator E t is the inverse time of A τ . Since U α (τ) is 1/α-similar, its inverse time S α (t) is α-similar. e self-similarity is a quite important property in financial market. In order to employ this feature, we choose Here, E t is also independent of B(τ). In fact, from Section 2, we can get the PDF p(x, t) of logprice of this underlying should satisfy the following FFPE: with initial condition p(x, 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 p(x, t) depends generally on the special shape of the underlying geometry. However, the interesting part of p(x, t) has the asymptotic behavior log p( here, u � 1/(1 − α/2), i.e., the solution p(x, 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 p(x, t) follows the fractional diffusion equation, where the Caputo fractional operator is given as e definition shows the Caputo fractional operator is a convolution integral, which implies B(S α (t)) is a longmemory 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 equation X(E t ) presents many important financial properties, such as selfsimilarity, leptokurtic feature, and long memory. So, we assert this model characterizes the price of underlying well.
Let r(t) be the return of the underlying from time 0 to time t. en, we have From the discussion, in the Section 2, we know the PDF of the returns is leptokurtosis and fat-tailed, given as 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. Let C(t, S t ) 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 C(T, S T ) � (S T − K) + .

Theorem 2.
Assume that the underlying is driven by the subordinated geometric Brownian motion S t � X(E t ), then the European call option C(t, S t ) can be given as where the function N(x) is the cumulative probability distribution function for standard normal distribution, and Proof. We assume that the riskless bond price dynamics satisfies where r is the riskless interest rate. From Taylor's formula, we have Since the B(τ) is independent of E t , we get where we have used the property of S α (t), i.e., S α (t) is a α-similar process and E[S n α (t)] � t nα Γ(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 "selffinancing" investment strategy, involving the underlying and an option on this underlying So, by substituting equations (30)-(33) into (29), we have where σ 2 (t) � (1 + t α− 1 /Γ(α))σ 2 . So, the solution of (34) with boundary condition C(T, S T ) � (S T − K) + can be given as where the function N(x) is the cumulative probability distribution function for standard normal distribution, and Substituting σ 2 (t) 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. e return of the underlying assets is supposed to be heavy tailed, which means the underlying assets is more risky. e risk also increases with the increase of α, that leads to the smaller the parameter α is, the higher the option price will be. e differences between this option pricing formula and the classical Black-Scholes model for different expiration times are presented in Figure 3. One can observe that the price differences become smaller as the maturity time increases. at is caused by the fact that the α-stable term dominates in the short time and the drift term dominates in the long time. is 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. e empirical analysis on SPDR S&P500 ETF call option (SPY) is introduced to verify our results. e 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/4n ln 2) { n i�0 (S ih − S il ) 2 } 1/2 , where S ih and S il 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.

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 lefthand side. e 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. e 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. en, we apply this model to option pricing problem. e 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. e 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.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that no conflicts of interest exist regarding this manuscript.