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.
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 . When , the value characterizes a superdiffusive process, a subdiffusive one, and 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 , and the second mechanism introduces the trapping events and is represented by the so-called inverse α-stable subordinator The subordinated process combines both mechanisms and gives the subdiffusive dynamics. The inverse α-stable subordinator is defined in the following way:where is a strictly increasing α-stable process. Here, the α-stability means the has independent increment and , in law. And, its Laplace transform is given by . Since is a pure-jump process, then, for every jump of , there is a corresponding flat period of its inverse These heavy-tailed flat periods of represent long waiting times in which the subdiffusive particle gets immobilized in the trap. As approaches 1, the random time reduces to the normal time t. In other words, the probability density function (PDF) of the subordinated process 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 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. is usually called price return. The model is a geometric random walk with drift μ and diffusion σ. 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 Section 2, 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: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 , besides the regular time derivative on the left-hand side. The convolution in equation (3) leads to memory effects. Here, can be any continuous function, and is the total density of tracer, including mobile and immobile phases.
Let be an increasing Lévy motion with Laplace transform . is the increasing Lévy motion with drift, . Let be the inverse time of , defined as denotes the Laplace transform of the function . Then, we have the following theorem.
Theorem 1. The PDF of the subordinated process is the stochastic representation of the modified advection-dispersion equation (equation (3)) for , where the parent process is defined asHere is independent of .
Proof. Since the Laplace transform of is given as , then with for all τ. Firstly, we establish the relation between the PDF of and the PDF of . From the definition of , we have [25, 26]; therefore,So, the Laplace transform of can be expressed aswhere we have used the property that for . Using the total probability formula and the independence between and , we get the PDF of , given bywhere is the PDF of the parent process . So, the Laplace transform of the above equation yieldsThus, the following relation between and holdsSince the process is given by the Itô stochastic differential equation (3), its PDF obeys the classical advection-dispersion equation [27]:The Laplace transform of the above equation with respect to τ yieldsBy changing the variable k to and using the relation (10), the above equation yieldsComparing the above equation with the Laplace transform of equation (3), we get if , 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, if is given, we can calculate by using equation (8), and then inverting the Laplace transform in , the expression of can be obtained. Substituting 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. .
We can get and . So, from equation (7), we have , and then . Applying inverse Laplace transform on yieldsSubstituting it into equation (10), we haveTaking advantage of the relation between and , we can get the following evaluation formula for the mean square displacement:Therefore,So, if , the process is a normal diffusion.
Case 2. .
For this case, equation (3) yieldswhere the operator is the Caputo fractional derivative of order α [28].
Following the same procedure shown in Case 1, we can get , so is a strictly increasing α-stable Lévy motion. As we know, the PDF of inverse time α-stable Lévy motion can be expressed in the form of a Fox function, i.e., (see [29]). By using the relation between and , after some calculations, we can getwhere 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 of equation (7), and then inverting the Laplace transform, we havewhere is the Mittag-Leffler function [31]. Here, we have used the equality:From equation (20), we can know if , the model resembles a normal diffusion for , and since , the model resembles a subdiffusion for . This result can be obtained in another way, to see this, note that in distribution and grows faster than for , so the stable subordinator dominates as and the drift dominates as . From the definition of and , we also should note that and are equivalent to and , respectively. By the way, as we know, the generalized Einstein relation is defined by [32]. From equations (8) and (16), we should notice that if the Einstein relation for the parent process holds, then it also holds for the subordinated process.

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 motion , where the parent process is given asand the subordinator is the inverse time of . Since is -similar, its inverse time is α-similar. The self-similarity is a quite important property in financial market. In order to employ this feature, we choose . Here, is also independent of .
In fact, from Section 2, we can get the PDF of log-price of this underlying should satisfy the following FFPE:with initial condition . Here, the fractional derivative is also the Caputo type. In [33], Ren et al. have proved that the detailed structure of the solution depends generally on the special shape of the underlying geometry. However, the interesting part of has the asymptotic behavior , where ; here, , i.e., the solution 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 follows the fractional diffusion equation, where the Caputo fractional operator is given as
The definition shows the Caputo fractional operator is a convolution integral, which implies 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 equation 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.
Let be the return of the underlying from time 0 to time t. Then, we have
From the discussion, in the Section 2, we know the PDF of the returns is leptokurtosis and fat-tailed, given aswhere is showed in equation (19) with and . 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 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 .
Theorem 2. Assume that the underlying is driven by the subordinated geometric Brownian motion , then the European call option can be given aswhere the function is the cumulative probability distribution function for standard normal distribution, and
Proof. We assume that the riskless bond price dynamics satisfieswhere r is the riskless interest rate. From Taylor’s formula, we haveSince the is independent of , we getwhere we have used the property of , i.e., is a α-similar process and [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 underlyingSo,by substituting equations (30)–(33) into (29), we havewhere . So, the solution of (34) with boundary condition can be given aswhere the function is the cumulative probability distribution function for standard normal distribution, andSubstituting 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 .
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 Figure 3. 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 , where and are the highest and lowest prices of the same day, respectively. We get the parameter . 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.



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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that no conflicts of interest exist regarding this manuscript.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 11801288), Zhejiang Provincial Natural Science Foundation (nos. LY17A010020 and LY20G030025), and Ningbo Natural Science Foundation (no. 2017A610130).