Generalised geometric Brownian motion: Theory and applications to option pricing
Viktor Stojkoski, Trifce Sandev, Lasko Basnarkov, Ljupco Kocarev, Ralf Metzler
aa r X i v : . [ q -f i n . P R ] O c t Generalised geometric Brownian motion: Theory andapplications to option pricing
Viktor Stojkoski , , Trifce Sandev , , , Lasko Basnarkov , , Ljupco Kocarev , and Ralf Metzler , *1 Faculty of Economics, Ss. Cyril and Methodius University, 1000 Skopje,Macedonia Research Center for Computer Science and Information Technologies,Macedonian Academy of Sciences and Arts, Bul. Krste Misirkov 2, 1000 Skopje,Macedonia Institute of Physics & Astronomy, University of Potsdam, D-14776Potsdam-Golm, Germany Institute of Physics, Faculty of Natural Sciences and Mathematics, Ss. Cyril andMethodius University, Arhimedova 3, 1000 Skopje, Macedonia Faculty of Computer Science and Engineering, Ss. Cyril and MethodiusUniversity, P.O. Box 393, 1000 Skopje, MacedoniaNovember 3, 2020
Abstract
Classical option pricing schemes assume that the value of a financial asset follows a ge-ometric Brownian motion (GBM). However, a growing body of studies suggest that a simpleGBM trajectory is not an adequate representation for asset dynamics due to irregularities foundwhen comparing its properties with empirical distributions. As a solution, we develop a gener-alisation of GBM where the introduction of a memory kernel critically determines the behaviorof the stochastic process. We find the general expressions for the moments, log-moments, andthe expectation of the periodic log returns, and obtain the corresponding probability densityfunctions by using the subordination approach. Particularly, we consider subdiffusive GBM(sGBM), tempered sGBM, a mix of GBM and sGBM, and a mix of sGBMs. We utilise theresulting generalised GBM (gGBM) to examine the empirical performance of a selected groupof kernels in the pricing of European call options. Our results indicate that the performance ofa kernel ultimately depends on the maturity of the option and its moneyness. * Corresponding author: [email protected] Introduction
Geometric Brownian motion (GBM) frequently features in mathematical modeling. The advantageof modelling through this process lies in its universality, as it represents an attractor of more com-plex models that exhibit non-ergodic dynamics [1–3]. As such, GBM has been used to underliethe dynamics of a diverse set of natural phenomena including the distribution of incomes, bodyweights, rainfall, fragment sizes in rock crushing processes, etc [4, 5]. Nevertheless, perhaps thebest-known application of GBM is in finance, and in particular in terms of the Black-Scholes (BS)model (or Black-Scholes-Merton model) [6–8] for the pricing of European options.By construction, GBM is a simple continuous-time stochastic process in which the logarithm ofthe randomly varying quantity of interest follows a Brownian motion with drift. Its non-ergodicityis manifested in the difference between the growth rate observed in an individual trajectory and theensemble average growth [9]. The time-averaged growth rate is dependent on both the drift andthe randomness in the system, whereas the ensemble growth rate is solely dependent on the drift.If only a single system is to be modeled, on the long run only the time-averaged growth rate, isobserved. This is naturally the case in financial market dynamics, for which only single time seriesexist, and where individual realisations would be expected to be distinctly disparate [10].Moreover, GBM is closely related to the problem of heterogeneous diffusion and turbulent dif-fusion, which are represented by the inhomogeneous advection-diffusion equation with position-dependent diffusion coefficient D ( x ) and velocity field v ( x ) . It is well known that at turbulent diffu-sion the contaminant spreads very fast. For the case of Richardson diffusion the position-dependentdiffusion coefficient behaves as D ( x ) ∼ x / and the relative mean squared displacement (MSD)scales as h x ( t ) i ∼ t [11]. However, the fast spread of contaminants can be essentially increaseddue to multiplicative noise, such that the MSD grows exponentially with time [12, 13].Notably, in a variety of cases GBM has failed to reproduce the properties of real asset prices.For instance, by definition, GBM is not able to adequately reproduce fat tailed distributions ofvarious characteristics widespread in nature [14]. As a solution, three alternating theories havebeen proposed: i) stochastic volatility [15–17]; ii) utilising stochastic processes in which thenoise follows a fat-tailed distribution [18–23]; and iii) generalisations of GBM based on subdiffu-sion [24–26]. In the first approach the volatility is a stochastic process itself. The second approachintuitively leads to the observation of log returns which follow a fat-tailed distribution. The lastapproach, differently from the first two views, assumes anomalous price dynamics. Concretely, theobservation that the distribution of log returns is fat tailed, can be attributed to prolonged periodsin which the price of the asset exhibits approximately constant extreme values. These constantperiods can be considered as trapping of particles, as is done in physical systems which mani-fest anomalous diffusion (subdiffusion) [27, 28]. While the resulting subdiffusive GBM (sGBM)is able to easily reproduce real-life properties, the literature lacks a extensive study in which theexact empirical characteristics of the subdiffusive model are presented.The purpose of this paper is to propose a unifying framework for the application of subdiffusiveGBM models in option pricing. We do this by developing the so-called generalised GBM (gGBM).gGBM is a stochastic process whose behavior is critically determined by a memory kernel. Bychoosing the appropriate kernel, we recover the standard GBM and the typically used subdiffusiveGBM models [24–26]. To understand the behavior of gGBM under various kernels, we performa detailed analysis and show that the dynamics of the model can be easily adjusted to mimicperiods of constant prices and/or fat-tailed observations of returns, thus corresponding to realistic2cenarios. More importantly, we utilise the properties of the model to investigate its capability topredict empirical option values. We find that the performance of a kernel ultimately depends onthe parameters of the option, such as its maturity and its moneyness. The first property describessimply the time left for the option to be exercised, wheres the second characteristic depicts therelative position of the current price with respect to the strike price of the option. On the firstsight, this conclusion appears intuitive – obviously the known information for the properties of theasset greatly impacts its price, the observation that a slight change in the known information maydrastically change the dynamics suggests that there is a need in the option pricing literature formodels that easily allow for such structural changes. We believe that the resolution to this issuelies in applying the concepts of time-averaging and ergodicity breaking to modeling financial time-series, and our gGBM framework offers a computationally inexpensive and efficiently tractablesolution.The paper is organised as follows. In Section 2 we provide an overview of GBM in the BSmodel and its use in option pricing. We also give detailed results for the so-called sGBM in termsof fractional Fokker-Planck equation and its corresponding continuous time random walk (CTRW)model. In Section 3 we present gGBM and describe its properties by using the subordination ap-proach. In particular, we derive the corresponding Fokker-Planck equation with a memory kerneland obtain the respective moments and log-moments. The general function used in the L´evy ex-ponent occurs as a memory kernel in the Fokker-Planck equation, which allows us to recover thepreviously known results for GBM and sGBM. We consider generalisations of GBM and sGBMby introducing tempered sGBM, a mix of GBM and sGBM, as well as a mix of sGBMs. An em-pirical example of application of the gGBM in option pricing is presented in Section 4. Section 5summarises our findings. In the Appendices we give detailed calculations as well as derivation ofthe Fokker-Planck equation for the gGBM within the CTRW theory.
GBM has been applied in a variety scientific fields [9, 29–33]. Mathematically, it is represented bythe Langevin equation dx ( t ) = x ( t ) [ µ dt + σ dB ( t )] , x = x ( ) , (1)where x ( t ) is the particle position, µ is the drift, σ > B ( t ) represents standardBrownian motion. The solution to Eq. (1) in the Itˆo sense is x ( t ) = x e ( µ − σ ) t + σ B ( t ) , x = x ( ) > . (2)When the dynamics of the asset price follows a GBM, then a risk-neutral distribution (proba-bility distribution which takes into account the risk of future price fluctuations) can be easily found3y solving the corresponding Fokker-Planck equation to Eq. (1), ∂∂ t f ( x , t ) = − µ ∂∂ x x f ( x , t ) + σ ∂ ∂ x x f ( x , t ) , (3)with initial condition f ( x , t = ) = δ ( x − x ) . The solution of Eq. (3) is the famed log-normaldistribution f ( x , t ) = x √ πσ t × exp − [ log x − log x − ¯ µ t ] σ t ! . (4)where ¯ µ = µ − σ / h x ( t ) i = x e µ t , (5)and h x ( t ) i = x e ( σ + µ ) t , (6)respectively, and thus, the variance becomes h x ( t ) i − h x ( t ) i = x e µ t (cid:16) e σ t − (cid:17) . (7)The exact derivation of the GBM distribution and its moments is given in Appendix A.Evidently, in GBM the diffusion coefficient scales proportionally with the square of the positionof the particle, i.e., D ( x ) = σ x /
2, and thus the MSD has an exponential dependence on time. Amore convenient measure instead of the MSD for geometric processes is the behaviour of theexpectation of the logarithm of x ( t ) . In the case of GBM the expectation of the logarithm of theparticle position has a linear dependence on time. This can be shown by calculation of the log-moments h log n x ( t ) i = R ∞ log n x P ( x , t ) dx , see Appendix A. The mean value of the logarithm of x ( t ) becomes h log x ( t ) i = h log x i + ¯ µ t , (8)from where for the expectation of the periodic log return with period ∆ t , one finds1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ¯ µ = ddt h log x ( t ) i . (9)The second log-moment is given by h log x ( t ) i = h log x i + ¯ µ t + µ h log x i t + σ t , (10) This Fokker-Planck equation corresponds to the Itˆo interpretation of the multiplicative noise. There are alsoStratonovich and Klimontovich-H¨anggi interpretations, for which the corresponding Fokker-Planck equations areslightly different, see Refs. [12, 34]. In finance math literature the Itˆo convention is the standard interpretation. h log x ( t ) i − h log x ( t ) i = σ t . (11) As previously said, perhaps the best-known application of GBM is in finance, and in particular theBS model for pricing of European options. Formally, a European option is a contract which givesthe buyer (the owner or holder of the option) the right, but not the obligation, to buy or sell anunderlying asset or instrument x ( T ) at a specified strike price K on a specified date T . The sellerhas the corresponding obligation to fulfill the transaction – to sell or buy – if the buyer (owner)“exercises” the option. An option that conveys to the owner the right to buy at a specific priceis referred to as a call ; an option that conveys the right of the owner to sell at a specific priceis referred to as a put . Here we are going to consider the valuation of call options, denoted as C BS ( x , t ) , with the note that the derived results easily extend to put options.In the modeling of financial assets, a standard assumption is that there is a risk-neutral distribu-tion f ( x , t ) for the price of the asset. This measure is simply a probability distribution which takesinto account the risk of future price fluctuations. Once a risk-neutral distribution is assigned, thevalue of the option is obtained by discounting the expectation of its value at the maturity T withrespect to that distribution [6, 35], i.e., C BS ( x , t ) = e − r ( T − t ) Z ∞ K ( x ( T ) − K ) f ( x , T , | x , ) dx , (12)where r is the risk-free rate of return and x is the asset price at the beginning ( t = K the option would not be exercised (i.e. its value is 0).Eqs. (12) and (4) can be combined to derive an analytical formula for the value of the calloption in the BS model for the GBM as C BS ( x , T , K , t ) = N ( d ) x ( t ) − N ( d ) Ke − ( µ − σ )( T − t ) (13) d = σ √ T − t (cid:20) log x ( t ) K + µ ( T − t ) (cid:21) (14) d = d − σ √ T − t , (15)where N ( x ) = √ π R x − ∞ e − u / du is the cumulative distribution function of the Gaussian distributionwith zero mean and unit variance. Put simply, the two terms in the BS formula describe the currentprice of the asset weighted by the probability that the investor will exercise its option at time t andthe discounted price of the strike price weighted by its exercise probability. The terms d , can beseen as measures of the moneyness of the option and N ( d , ) as probabilities that the option willexpire while its value is in the money. The neat BS formulation has allowed the model to be widelyapplied in both theoretical investigations and empirical implementations. However, the BS modelhas failed to adequately reproduce a plethora of real world properties.The European option C BS ( x , t ) (12) is a solution of the Black-Scholes equation, see for example536], (cid:18) ∂∂ t + σ x ∂ ∂ x − r + r x ∂∂ x (cid:19) C BS ( x , t ) = , (16)with initial condition C BS ( x , T ) = max { x − K , } , x ≥
0, and boundary conditions C BS ( x = , t ) = t ≥ T , and C BS ( x → ∞ , t ) → x . By using t = T → t , one finds the equation ∂∂ t C BS ( x , t ) = (cid:18) σ x ∂ ∂ x − r + r x ∂∂ x (cid:19) C BS ( x , t ) . (17)with initial condition C BS ( x , t = ) = max { x − K , } , x ≥
0, and boundary conditions C BS ( x = , t ) = t ≥
0, and C ( x → ∞ , t ) → x .In particular, theoretically predicted option prices with fixed values for drift µ and volatility σ via the BS model are known to significantly deviate from their respective market values in aplethora of cases. To deal with this problem, extensions of the BS model have emerged, whichinclude combination of the GBM with jumps [8, 37], or with stochastic volatility [38, 39]. One of the reasons why the standard GBM is not able to explain empirical data is because it failsto explain periods of constant prices which appear on markets with low number of transactions.The price in these constant periods can be described as a trapped particle in physical systems thatmanifest anomalous diffusion (subdiffusion) [27, 28]. To deal with this problem, the so-called sub-diffusive
GBM (sGBM) has been developed [24], by using the subordination approach. The cor-responding equation for the sGBM becomes the following fractional Fokker-Planck equation [24](see also [25]) ∂∂ t f α ( x , t ) = RL D − α t (cid:20) − µ ∂∂ x x f α ( x , t ) + σ ∂ ∂ x x f α ( x , t ) (cid:21) , (18)where RL D ν t f ( t ) = Γ ( − ν ) ddt Z t ( t − t ′ ) − ν f ( t ′ ) dt ′ (19)is the Riemann-Liouville fractional derivative of order 0 < ν < . To avoid the strange initialcondition, alternatively, we could use the integral version of the equation f α ( x , t ) − g α ( x , ) = RL D − α t (cid:20) − µ ∂∂ x x f α ( x , t ) + σ ∂ ∂ x x f α ( x , t ) (cid:21) , (20) The Laplace transform of the Riemann-Liouville fractional derivative of a given function reads L { RL D ν t f ( t ) } ( s ) = s ν L { f ( t ) } ( s ) − RL I − ν t f ( +) , where RL I ν t f ( t ) = Γ ( ν ) R t ( t − t ′ ) ν − f ( t ′ ) dt ′ is the Riemann-Liouville fractional integral. L (cid:8) RL D − α t f ( t ) (cid:9) ( s ) = s − α L { f ( t ) } ( s ) . In Ref. [25] the time fractional Fokker-Planckequation (18) for sGBM is derived within the CTRW theory for a particle on a geometric lattice inpresence of a logarithmic potential.Here we note that the fractional Fokker-Planck equation (18) can be obtained by using theLangevin equation approach [41], i.e., by considering a CTRW model described by a coupledLangevin equations [42], ddu x ( u ) = µ x ( u ) + σ x ( u ) ξ ( u ) , (21) ddu T ( u ) = ζ ( u ) . (22)Therefore, x ( t ) is parametrised in terms of the number of steps u , and the connection to the phys-ical time t is given by T ( u ) = R u τ ( u ′ ) du , where τ ( u ) is a total of individual waiting times τ for each step. In mathematical terms this is called subordination [43–45]. The noise ξ ( u ) is awhite noise with zero mean and correlation h ξ ( u ) ξ ( u ′ ) i = δ ( u − u ′ ) , while ζ ( u ) is one-sided α -stable L´evy noise with the stable index 0 < α <
1. The inverse process S ( t ) of the one-sided α -stable Levy process T ( u ) with characteristic function h e − s T ( u ) i = e − s α u is given by S ( t ) = inf { u > T ( u ) > t } , i.e., it represents a collection of first passage times [41]. TheCTRW is defined by the subordinated process X ( t ) = x ( S ( t )) . The PDF h ( u , t ) of the inverseprocess S ( t ) can be found from the relation [41] h ( u , t ) = − ∂∂ u Θ ( t − T ( u )) , (23)where Θ ( z ) is the Heaviside theta function. The Laplace transform then yieldsˆ h ( u , s ) = − ∂∂ u s (cid:28) Z ∞ δ ( t − T ( u )) e − st dt (cid:29) = − ∂∂ u s h e − s T ( u ) i = − ∂∂ u s e − s α u = s α − e − s α u . (24)Therefore, f α ( x , t ) = h δ ( x − X ( t )) i = h δ ( x − X ( S ( t )) i = R ∞ f ( x , u ) h ( u , t ) dt , from where onecan easily arrive to the fractional Fokker-Planck equation (18).The mean value for sGBM is given by [25, 42] h x ( t ) i = x E α ( µ t α ) , (25)where E α ( z ) is the one parameter Mittag-Leffler (ML) function [27, 40] E α ( z ) = ∞ ∑ k = z k Γ ( α k + ) , (26)with ( z ∈ C ; ℜ ( α ) > ) , and Γ ( · ) is the Gamma function. The ML function is a generalization of The Laplace transform of the one parameter ML function reads L { E α ( at α ) } ( s ) = s α − s α − a . E ( z ) = e z . The asymptotic behavior of the mean is given by h x ( t ) i ≃ x ( + µ t α / Γ ( + α ) ∼ e µ t α / Γ ( + α ) , t ≪ , α − e µ / α t , t ≫ . (27)The MSD also is given through the one parameter ML function [25, 42] h x ( t ) i = x E α (cid:0) ( σ + µ ) t α (cid:1) ≃ h x ( ) i ( + ( σ + µ ) t α / Γ ( + α ) ∼ e ( σ + µ ) t α / Γ ( + α ) , t ≪ , α − e ( σ + µ ) / α t , t ≫ . (28)From here one concludes that the sGBM is an exponentially fast process.The first log-moment has the form [25] h log x ( t ) i = h log x i + ¯ µ Z t t ′ α − Γ ( α ) dt ′ = h log x i + ¯ µ t α Γ ( + α ) , (29)which gives a power-law dependence with respect to time of the expectation of the log return withperiod ∆ t , i.e., [25] 1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ¯ µ t α − Γ ( α ) . (30)Such models have been used, for example, to explain the dynamics of an asset before a marketcrash [47]. The second log-moment becomes [25] h log x ( t ) i = h log x i + (cid:2) µ h log x i + σ (cid:3) t α Γ ( + α ) + µ t α Γ ( + α ) . (31)from where for the log-variance one finds [25] h log x ( t ) i − h log x ( t ) i = σ t α Γ ( + α ) + ¯ µ (cid:20) Γ ( + α ) − Γ ( + α ) (cid:21) t α , (32)which in the long time limit scales as t α (0 < α < t for regularGBM ( α = In this section we consider a generalization of GBM, under which the standard and subdiffusiveGBM arise as special cases, by using the subordination approach. The continuous time random For the short time limit we use the first two terms from the series expansion of the ML function (26), while forthe long time limit we apply its asymptotic expansion formula [40, 46], E α ( z ) ≃ α e z / α , z ≫
1. Here we note that theasymptotic behavior of the ML function with negative argument has a power-law form, i.e., E α ( − z α ) ≃ z − α Γ ( − α ) for z ≪ < α < h e − s T ( u ) i = e − ˆ Ψ ( s ) u , with ˆ Ψ ( s ) = / ˆ η ( s ) . The generalisation of GBM which we consider is in the form of the stochastic process X ( t ) = x ( S ( t )) , (33)where X ( t ) is the generalised GBM (gGBM) , S ( t ) = inf { u > T ( u ) > t } is the operationaltime, and T ( u ) is an infinite divisible process, i.e., a strictly increasing L´evy motion with h e − s T ( u ) i = e − u ˆ Ψ ( s ) , and ˆ Ψ ( s ) is the L´evy exponent [24, 32, 51]. Here we consider ˆ Ψ ( s ) = / ˆ η ( s ) .Next we find the PDF of gGBM which subordinates the processes from the time scale t (phys-ical time) to the GBM on a time scale u (operational time). Therefore, the PDF P ( x , t ) of a givenrandom process X ( t ) can be represented as [24, 52–55] P ( x , t ) = Z ∞ f ( x , u ) h ( u , t ) du , (34)where f ( x , u ) satisfies the Fokker-Planck equation (3) for the standard GBM. The function h ( u , t ) is the PDF subordinating the random process X ( t ) to the standard GBM. In the Laplace space,Eq. (34) reads ˆ P ( x , s ) = L { P ( x , t ) } = Z ∞ e − st P ( x , t ) dt = Z ∞ f ( x , u ) ˆ h ( u , s ) du , (35)where ˆ h ( u , s ) = L { h ( u , t ) } . By consideringˆ h ( u , s ) = ˆ Ψ ( s ) s e − u ˆ Ψ ( s ) = s ˆ η ( s ) e − u ˆ η ( s ) , (36)we then have ˆ P ( x , s ) = s ˆ η ( s ) Z ∞ f ( x , u ) e − u ˆ η ( s ) du = s ˆ η ( s ) ˆ f (cid:18) x , η ( s ) (cid:19) . (37)By Laplace transform of the Fokker-Planck equation (3) for the GBM, and using relation (37), one The current process should not be confused with the Pagnini-Mainardi generalised grey Brownian motion, seeRef. [48–50]. P ( x , s ) satisfies s ˆ P ( x , s ) − P ( x , ) = s ˆ η ( s ) (cid:20) − µ ∂∂ x x ˆ P ( x , s ) + σ ∂ ∂ x x ˆ P ( x , s ) (cid:21) . (38)After inverse Laplace transform we arrive at the generalised Fokker-Planck equation (see Refs. [42,51] where one-sided α -stable waiting times are considered in detail) ∂∂ t P ( x , t ) = ∂∂ t Z t η ( t − t ′ ) (cid:20) − µ ∂∂ x xP ( x , t ′ ) + σ ∂ ∂ x x P ( x , t ′ ) (cid:21) dt ′ , (39)where η ( t ) is a so-called memory kernel. One observes that for η ( t ) = η ( t ) = t α − Γ ( α ) at the time fractional Fokker-Planck equation(18) for the sGBM. From Eqs. (35) and (36), we find for the PDF in the Laplace domain, see alsoRef. [42],ˆ P ( x , s ) = Z ∞ x √ πσ u × exp (cid:18) − (cid:2) log x − log x − ¯ µ u (cid:3) σ u (cid:19) s ˆ η ( s ) e − u ˆ η ( s ) du = / [ s ˆ η ( s )] x p ¯ µ + σ / ˆ η ( s ) exp (cid:16) − log x − log x σ hp ¯ µ + σ / ˆ η ( s ) − ¯ µ i(cid:17) , x > x , , x = x , exp (cid:16) log x − log x σ hp ¯ µ + σ / ˆ η ( s ) + ¯ µ i(cid:17) , x < x , (40) Remark 1:
Here we note that there are restrictions on the choice of the memory kernel η ( t ) since the PDF (34) should be non-negative. From the subordination integral it follows that thesubordination function h ( u , t ) should be non-negative, which, according to the Bernstein theorem,means that its Laplace transform (36) should be a completely monotone function [56]. Therefore,the PDF (34) will be non-negative if 1 / [ s ˆ η ( s )] is a completely monotone function, and 1 / ˆ η ( s ) is aBernstein function, see Refs. [57, 58]. Remark 2:
We note that Eq. (38) can be written in an equivalent form as Z t γ ( t − t ′ ) ∂∂ t ′ P ( x , t ′ ) dt ′ = − µ ∂∂ x xP ( x , t ) + σ ∂ ∂ x x P ( x , t ) , (41)where the memory kernel γ ( t ) is connected to η ( t ) in Laplace space as γ ( s ) = / [ s η ( s )] [57]. Fromthis relation we find that for GBM ( η ( t ) =
1, i.e., ˆ η ( s ) = / s ) the memory kernel γ ( t ) is given by γ ( t ) = L − (cid:8) s − ˆ η − ( s ) (cid:9) = L − { } = δ ( t ) . For sGBM ( η ( t ) = t α − / Γ ( α ) , i.e., ˆ η ( s ) = s − α )the memory kernel becomes γ ( t ) = L − (cid:8) s α − (cid:9) = t − α / Γ ( − α ) , and thus Eq. (41) reads C D α t P ( x , t ) = − µ ∂∂ x xP ( x , t ) + σ ∂ ∂ x x P ( x , t ) , (42)10here C D ν t f ( t ) = Γ ( − ν ) Z t ( t − t ′ ) − ν ddt ′ f ( t ′ ) dt ′ (43)is the Caputo fractional derivative of order 0 < ν < . We note that with the appropriaterestrictions for η ( t ) and γ ( t ) both formulations are equivalent. Remark 3:
For η ( t ) = t α − Γ ( α ) , 0 < α <
1, gGBM corresponds to sGBM. From the subordinationapproach one finds [24]ˆ h ( u , s ) = s α − e − us α = s α − H , , (cid:20) u s α (cid:12)(cid:12)(cid:12)(cid:12) − ( , ) (cid:21) , (44)where H m , np , q ( z ) is the Fox H -function, see Appendix D. By inverse Laplace transform (D3) weobtained [24] h ( u , t ) = L − (cid:8) ˆ h ( u , s ) (cid:9) = t − α H , , (cid:20) ut α (cid:12)(cid:12)(cid:12)(cid:12) ( − α , α )( , ) (cid:21) = u H , , (cid:20) ut α (cid:12)(cid:12)(cid:12)(cid:12) ( , α )( , ) (cid:21) , (45)where we applied property (D4). The solution in Laplace space then becomesˆ P ( x , s ) = Z ∞ x √ πσ u × exp (cid:18) − (cid:2) log x − log x − ¯ µ u (cid:3) σ u (cid:19) s α − e − us α du = s α − x p ¯ µ + σ s α × exp (cid:16) − log x − log x σ hp ¯ µ + σ s α − ¯ µ i(cid:17) , x > x , , x = x , exp (cid:16) log x − log x σ hp ¯ µ + σ s α + ¯ µ i(cid:17) , x < x , (46)which is obtained in Ref. [42] in a similar way. From here we can plot the PDF by using numericalinverse Laplace transform techniques. If we consider that the asset price follows a gGBM, the generalised BS (gBS) formula for theoption price is [51] C gBS ( x , t ) = h e − r ( S ( T ) − t ) ( x ( S ( T )) − K ) i x = Z ∞ C BS ( x , u ) h ( u , T ) du , (47) The Laplace transform of the Caputo derivative of a given function reads L { C D ν t f ( t ) } ( s ) = s ν L { f ( t ) } ( s ) − s ν − f ( +) . C BS ( x , t ) is taken from the BS formula (15), and h ( x , T ) is the subordination function definedby Eq. (36) in the Laplace domain. By Laplace transform one findsˆ C gBS ( x , s ) = s ˆ η ( s ) ˆ C BS ( x , / ˆ η ( s )) . (48)Therefore, from Eq. (17), the corresponding equation for the option price becomes [42] ∂∂ t C gBS ( x , t ) = ∂∂ t Z t η ( t − t ′ ) (cid:18) σ x ∂ ∂ x − r + r x ∂∂ x (cid:19) C gBS ( x , t ′ ) dt . (49) The n th moment h X n ( t ) i = R ∞ x n P ( x , t ) dx can be calculated by multiplying both sides of Eq. (39)by x n and integration over x , see Appendix C. In the Laplace domain, this results in h ˆ X n ( s ) i = s − − ˆ η ( s ) h σ n ( n − ) + µ n i h x n i . (50)From this result we reproduce the normalization condition h x ( t ) i = h x i =
1. The general resultsfor the mean value ( n =
1) and the MSD ( n =
2) in terms of the memory kernel become [42], h ˆ X ( s ) i = x s − − µ ˆ η ( s ) , (51)and h ˆ X ( s ) i = x s − − ( σ + µ ) ˆ η ( s ) . (52)The log-moments h log n x ( t ) i = R ∞ log n x P ( x , t ) dx , can also be calculated exactly through thememory kernel, see Appendix C. The normalization condition is satisfied, i.e., h log x ( t ) i = h log x ( t ) i = h log x i + ¯ µ Z t η ( t ′ ) dt ′ . (53)From here, we find for the expectation of the periodic log return with period ∆ t ∆ t h log ( x ( t + ∆ t ) / x ( t )) i = ¯ µ ∆ t Z t + ∆ tt η ( t ′ ) dt ′ = ¯ µ I ( t + ∆ t ) − I ( t ) ∆ t ∼ ∆ t → ¯ µη ( t ) , (54)where I ( t ) = R η ( t ) dt , i.e., I ′ ( t ) = η ( t ) . Therefore, the expectation of the periodic log returns12ehaves as the rate of the first log-moment,1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ddt h log x ( t ) i , (55)which is proportional to the memory kernel η ( t ) . Moreover, for the second log-moment we find h log x ( t ) i = h log x i + Z t η ( t − t ′ ) (cid:26) µ (cid:20) h log x i + ¯ µ Z t ′ η ( t ′′ ) dt ′′ (cid:21) + σ (cid:27) dt ′ , (56)from where the log-variance becomes h log x ( t ) i − h log x ( t ) i = σ Z t η ( t ′ ) dt ′ + ¯ µ " Z t η ( t − t ′ ) (cid:18) Z t ′ η ( t ′′ ) dt ′′ (cid:19) dt ′ − (cid:18) Z t η ( t ′ ) dt ′ (cid:19) . (57)From all these general formulas one can easily recover the previous results for the standardGBM ( η ( t ) =
1, i.e., ˆ η ( s ) = / s ) and sGBM ( η ( t ) = t α − / Γ ( α ) , i.e., ˆ η ( s ) = s − α , 0 < α < As an example for another memory kernel in gGBM we consider a power-law memory kernel withexponential truncation, η ( t ) = t α − Γ ( α ) e − t τ , (58)where τ is a characteristic crossover time scale, 0 < α <
1. Such forms are important in manyreal-world applications, in which the scale-free nature of the waiting time dynamics is broken atmacroscopic times t ≫ τ [57]. Therefore,ˆ η ( s ) = ( s + τ − ) − α , (59)where we use the shift rule of the Laplace transform, L { e − at f ( t ) } = ˆ F ( s + a ) , for ˆ F ( s ) = L { f ( t ) } .The mean value reads, h x ( t ) i = x L − (cid:26) s − − µ ( s + τ − ) − α (cid:27) ( t ) = x Z t e − t ′ / τ t ′− E α , (cid:0) µ t ′ α (cid:1) , (60)and the MSD h x ( t ) i = x Z t e − t ′ / τ t ′− E α , (cid:0) ( σ + µ ) t ′ α (cid:1) , (61)13here E α , β ( z ) = ∞ ∑ k = z k Γ ( α k + β ) (62) ( z , β ∈ C ; ℜ ( α ) > ) is the two parameter ML function [40] . From here, for the short time limitwe obtain the results for the sGBM h x ( t ) i ≃ x E α ( µ t α ) , (63) h x ( t ) i ≃ x E α (cid:0) ( σ + µ ) t α (cid:1) , (64)since the exponential truncation has no effect for short times, e − t / τ ≃ − t / τ , t / τ ≪
1. The longtime limit ( s τ ≪
1) yields h x ( t ) i ≃ x µ / α α ( µ / α − τ − ) h e ( µ / α − τ − ) t − i , (65)and h x ( t ) i ≃ x ( σ + µ ) α α (cid:0) [ σ + µ ] / α − τ − (cid:1) n e ( [ σ − µ ] / α − τ − ) t − o . (66)From the general result for the log-mean, we find h log x ( t ) i = h log x i + ¯ µ e − t / τ t α E , α + ( t / τ ) = h log x i + ¯ µτ α γ ( α , t / τ ) Γ ( α ) , (67)where γ ( a , z ) = R z t a − e − t dt = Γ ( a ) e − z z a E , a + ( z ) is the incomplete gamma function. For theexpectation of the periodic log return with period ∆ t we find1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ¯ µ t α − Γ ( α ) e − t / τ = ddt h log x ( t ) i . (68)This leads to a long run log return of 0, whereas on the short time scale the same observablebehaves in the same way as sGBM. As such, the model can be used to model early herd behaviorwhere the price of an asset grows simply as a consequence of investors following trends (short runbehavior), that last until the trade of the asset becomes congested (long run behavior). The second The Laplace transform of the two parameter ML function reads L (cid:8) t β − E α , β ( at α ) (cid:9) ( s ) = s α − β s α − a . Here we use the asymptotic expansion formula for the two parameter ML function E α , β ( z ) ≃ α e z / α z ( − β ) / α , z ≫ E α , β ( − z α ) ≃ z − α Γ ( β − α ) , z ≫ h log x ( t ) i = h log x i + (cid:2) µ h log x i + σ (cid:3) e − t / τ t α E , α + ( t / τ ) + µ e − t / τ t α E , α + ( t / τ ) , (69)from where the log-variance becomes h log x ( t ) i − h log x ( t ) i = µ e − t / τ t α E , α + ( t / τ )+ e − t / τ t α E , α + ( t / τ ) h σ − ¯ µ e − t / τ t α E , α + ( t / τ ) i . (70)Here we note that for t / τ ≪
1, the obtained results correspond to those obtained for sGBM, as itshould be since the exponential truncation has no influence on the process. We observe that on thelong run the log variance becomes constant, i.e., it is equal to σ τ α + ¯ µτ α .The subordination function in this case is given byˆ h ( u , s ) = ( s + τ − ) α s e − u ( s + τ − ) α = (cid:2) + ( s τ ) − (cid:3) ( s + τ − ) α − e − u ( s + τ − ) α , h ( u , t ) = e − t / τ H , , (cid:20) ut α (cid:12)(cid:12)(cid:12)(cid:12) ( , α )( , ) (cid:21) + τ Z t e − t ′ / τ H , , (cid:20) ut ′ α (cid:12)(cid:12)(cid:12)(cid:12) ( , α )( , ) (cid:21) dt ′ , (71)from where one can analyze the PDF P ( x , t ) . As another application, let us consider the combination of GBM and sGBM, represented by thememory kernel η ( t ) = w t α − Γ ( α ) + w , (72)where 0 < α < w + w =
1, and ˆ η ( s ) = w s − α + w s − . (73)This case combines both motions governed by Eq. (3) and (18). In this case, in a jump picturenormal GBM steps occur with weight w while power-law waiting time steps are realised withweight w .The mean value for this case is given by h x ( t ) i = x L − (cid:26) s − − µ ( w s − α + w s − ) (cid:27) ( t ) = x ∞ ∑ n = w n µ n t α n E n + , α n + ( w µ t )= x ∞ ∑ n = w n µ n t α n Γ ( α n + ) F ( n + α n + w µ t ) , (74)15here F ( a ; b ; z ) = ∑ ∞ k = ( a ) k ( b ) k z k k ! is the Kummer confluent hypergeometric function, and E γα , β ( z ) = ∞ ∑ n = ( γ ) n Γ ( α n + β ) z n n ! , (75)is the three parameter ML function [59], ( γ ) n = Γ ( γ + n ) / Γ ( γ ) is the Pochhammer symbol . Fromhere we see that for w = w = n = w = w =
0, yields h x ( t ) i = x ∞ ∑ n = µ n t α n Γ ( α n + ) = E α ( µ t α ) , (76)as it should be for the sGBM. For the second moment we find h x ( t ) i = x ∞ ∑ n = w n ( σ + µ ) n t α n E n + , α n + (cid:0) w ( σ + µ ) t (cid:1) . (77)Following the same procedure as previous, for the log-mean we find h log x ( t ) i = h log x i + ¯ µ (cid:20) w t α Γ ( α + ) + w t (cid:21) . (78)and for the expectation of the periodic log return with period ∆ t ,1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ¯ µ (cid:20) w t α − Γ ( α ) + w (cid:21) = ddt h log x ( t ) i . (79)This model introduces subdiffusive and trapping asset dynamics on short time scales (i.e., then thepart multiplied with w is much bigger), whereas on the long run we recover the standard GBMdynamics. The second log-moment yields h log x ( t ) i = h log x i + µ (cid:20) w t α Γ ( α + ) + w w t α + Γ ( α + ) + w t (cid:21) + (cid:8) µ h log x i + σ (cid:9) (cid:18) w t α Γ ( α + ) + w t (cid:19) , (80)from where the log-variance becomes h log x ( t ) i − h log x ( t ) i = σ (cid:18) w t α Γ ( α + ) + w t (cid:19) + ¯ µ w t α (cid:18) Γ ( α + ) − Γ ( α + ) (cid:19) + µ w w t α + (cid:18) Γ ( α + ) − Γ ( α + ) (cid:19) . (81)Similarly to the behavior of the first log moment, in the log variance, for short time scales the The Laplace transform of the three parameter ML function reads L n t β − E γα , β ( at α ) o ( s ) = s αγ − β ( s α − a ) γ . w w t α + .The subordination function for this case is given byˆ h ( u , s ) = w + w s − α e − uw s − + w s − α , (82)where the L´evy exponent is ˆ Ψ ( s ) = (cid:2) w s − + w s − α (cid:3) − . We may further analyze the case of a mix of two sGBM with different power-law memory func-tions, η ( t ) = w t α − Γ ( α ) + w t α − Γ ( α ) , (83)where 0 < α < α < w + w =
1, andˆ η ( s ) = w s − α + w s − α . (84)This situation corresponds to the case of two different groups of periods of constant prices. Forphysical systems, this situation means that the particles are trapped in traps with different waitingtimes [60], represented by the memory kernel (83).Therefore, for the mean we find h x ( t ) i = x L − (cid:26) s − − µ ( w s − α + w s − α ) (cid:27) ( t ) = x ∞ ∑ n = w n µ n t α n E n + α , α n + ( w µ t α ) , (85)while for the MSD we obtain h x ( t ) i = x ∞ ∑ n = w n ( σ + µ ) n t α n E n + α , α n + (cid:0) w ( σ + µ ) t α (cid:1) . (86)Similarly, the log-mean yields h log x ( t ) i = h log x i + ¯ µ (cid:20) w t α Γ ( α + ) + w t α Γ ( α + ) (cid:21) . (87)The expectation of the log return with period ∆ t , then becomes1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ¯ µ (cid:20) w t α − Γ ( α ) + w t α − Γ ( α ) (cid:21) = ddt h log x ( t ) i . (88)Since 0 < α < α <
1, on short times, the part of first sGBM dominates, whereas on long timesit is the characteristic of the second sGBM that determines the dynamics. The second log-moment17ecomes h log x ( t ) i = h log x i + µ (cid:20) w t α Γ ( α + ) + w w t α + α Γ ( α + α + ) + w t α Γ ( α + ) (cid:21) + (cid:8) µ h log x i + σ (cid:9) (cid:18) w t α Γ ( α + ) + w t α Γ ( α + ) (cid:19) , (89)and for the log-variance we find h log x ( t ) i − h log x ( t ) i = σ (cid:18) w t α Γ ( α + ) + w t α Γ ( α + ) (cid:19) + µ w w t α + α (cid:18) Γ ( α + α + ) − Γ ( α + ) Γ ( α + ) (cid:19) + ¯ µ w t α (cid:18) Γ ( α + ) − Γ ( α + ) (cid:19) + ¯ µ w t α (cid:18) Γ ( α + ) − Γ ( α + ) (cid:19) . (90)In this case, for short times, the kernel with the smaller exponent dominates the variance. Interest-ingly, for long times, this observable is determined by the magnitude of the larger exponent, whichis opposite from the previous kernel examples.For the mix of subdiffusive GBMs the subordination function becomesˆ h ( u , s ) = w s − α + w s − α e − uw s − α + w s − α , (91)where the L´evy exponent is ˆ Ψ ( s ) = [ w s − α + w s − α ] − .Figure 1(a), gives an intuitive illustration of the gGBM dynamics under various choices for thekernel. As argued, for standard GBM we observe smooth dynamics without periods of constantprices, whereas there is more turbulence in the asset price dynamics in the gGBM case. Theperiods of constant prices reproduced by gGBM depend in general on the time scale and, hence,the measuring units of the drift and volatility, with longer time scales also corresponding to longerperiods of constant prices. In Figure 1(b) and Figure 1(c) we plot, respectively, the numericalapproximations for the first moment and the MSD for GBM, sGBM, a mix of GBM and sGBMand a mix of sGBMs. One can easily notice the nonlinear behavior in the generalizations of GBM.For long times all gGBMs give exponential dependence of the first moment and the MSD ontime but with smaller slope than the one of GBM. Finally, Figure 1(d) gives the empirical PDFfor the logarithmic return at t =
1. For each of the studied generalizations of GBM, the PDF ischaracterised with fatter tails (which should increase as the α parameters increase), meaning thatit is more prone to producing values that fall far from the average. This can be easily observedas from the excess kurtosis present in each GBM generalization. This is exactly what makes thegGBM framework useful for understanding the statistical behavior of the asset price dynamics. To illustrate the power of the gGBM framework in the description of option pricing we utiliseempirical data of American options for two companies, Tesla (TSLA) and Apple (AAPL). By18 igure 1: gGBM properties. (a)
An example for simulated individual trajectories of gGBM for differentmemory kernels: standard GBM (blue solid line), sGBM (red dashed line), mix of standard GBM andsGBM (yellow dotted line), mix of sGBM (violet dot-dashed line). (b)
Numerical estimation for the firstmoment in GBM, sGBM, mix of standard GBM and sGBM and the mix of sGBM as a function of time. (b)
Same as (b) , only for the second moment. (c)
Empirical PDF for the logarithmic return at t = (a)-(c) In the simulations, µ = .
03 and σ = .
02. Moreover,for the sGBM case we set α = .
8, for the mix GBM- sGBM case we set α = . w = w = .
5, andfor the mix of sGBM case α = . α = . w = w = . definition, the dynamics of American options differ from European as they allow exercising ofthe option at any time before the option expires. Nevertheless, as given in Hull (2017), one canrely on the fact that American options on non-dividend-paying stocks have the same value as theirEuropean counterpart. This relation has allowed for the empirical examination of a pricing schemeof European options to be widely done via data for American ones.For our analysis we use the freely available data from the Nasdaq’s Options Trading Center.This dataset offers daily data free of charge for options of all companies quoted there. However,the options for most companies have small sample size. Therefore we have restricted the empiricalanalysis to Tesla and Apple, whose options are more frequently traded. In our estimations, the driftparameter µ is simply taken as the 3-Month Treasury Bill Secondary Market Rate at the date ofobservation. The noise parameter, on the other hand, was inferred from the values of the options19n the market as the value which produces the minimum squared error in their fit. In finance, thisis known as use the famous “implied volatility” approach.Let us now turn our attention to Fig. 2 where we use TSLA data gathered on 1st March 2018 onoptions which expire on 16th March 2018 to examine the dependence of the sGBM model on themoneyness of the option in predicting it. Moneyness describes the relative position of the currentprice of TSLA ( x ) with respect to the strike price of the option. An option whose strike price isequal to the current price of the asset is said to be at the money; if the strike price is larger than thecurrent price, the option is “out of the money”; and if the strike price is smaller than the currentprice, the option is described to be “in the money”. In Fig. 2 we vary the subdiffusion parameter α , and plot the absolute difference in the estimated option price C g and the observed option priceas a function of the strike price. We find that for in-the-money-options the best prediction is with α =
1, which corresponds to the BS model. However, as the strike price of the option nears theTSLA price, a transition occurs and α = α . Overall, as shown in the inset plot where we plot the mean squared error of the prediction asa function of α , this analysis suggests that the best prediction for the TSLA data is done with α which is around 0 .
25, thus highlighting the subdiffusive nature of the dynamics of the TSLA stock.Next, we use the AAPL data gathered on 28th February 2018 and examine how the maturity T affects the performance of the same sGBM model in predicting the option price. For this purpose,Fig. 3 depicts the mean squared error of the option price prediction as a function of the parameter α . We observe that, in general, the best prediction occurs when α =
1. This may suggest thatthe dynamics of the AAPL stock price is quite nicely explained with the BS model. However, wealso see that the mean squared error is highly dependent on the maturity, and even that for somematurity very low subdiffusive parameter values exhibit similar performance as the BS model.Hence, one might even argue that different gGBM kernels can lead to similar outcomes in thepricing of options, an interesting finding as such.Evidently, the performance of a kernel ultimately depends on the physical properties of theoption. On the first sight, this conclusion appears intuitive – obviously the known informationfor the properties of the asset greatly impacts its price, the observation that a slight change in theknown information may drastically change the dynamics suggests that there is a need in the optionpricing literature for models that easily allow for such structural changes. In this aspect, we believethat the generalised GBM approach offers a computationally inexpensive and efficiently tractablesolution to this issue. Consequently, we stress that a significant improvement of the description ofthe data in the gGBM framework can be achieved with comparatively few additional parameters.
We investigated the potential of GBM extensions based on subdiffusion to model and predict theprice of options. By assuming that the price of the asset underlying the option undergoes a subdif-fusive process, we introduced the gGBM framework as a potential model for its value.Similar to previous works on subdiffusive GBM models, the dynamics of a particular gGBMinstance is critically determined by a memory kernel. The advantage of gGBM comes in the flavorof allowing various forms for the functional form of the kernel. Depending on its choice, we20 igure 2:
Moneyness in sGBM.
The absolute difference between the predicted TSLA option price C g andits real value C as a function of the strike price of the option for various choices of α . The data is taken on1st March 2020 and describe the value of TSLA options which expire on 16th March 2020. The inset plotgives the mean squared error of the predictions as a function of α . may end up with asset price dynamics whose behavior significantly varies on the short time incomparison to its long run characteristics. This, in turn, may induce observations of the propertiesof the asset price that more closely mimic realistic behavior than standard GBM.We explored the ability of gGBM to fit and predict real option values. Our empirical analysisconfirmed the characteristics of gGBM, as we discovered that the performance of a certain choiceof memory kernel is uniquely determined by the parameters of the option, such as its maturity andits moneyness. Since each kernel produces, in general, different long run and short run dynamics,this suggests that time-averages play an important role in efficient pricing of options. Formally,time-averaging is essential for the analysis of a single time-series (or a set of few) which is charac-terised with non-ergodic dynamics. The non-ergodicity creates non-equilibrium dynamics which,consequently, makes studies of the ensemble behavior irrelevant. This leads to the introduction ofnovel strategies for analysing financial data [9, 61].In line with our conclusions, we believe that the next step in uncovering the properties of gGBMis demonstrating the ergodicity breaking of the process. Since multiplicative processes are fre-21 igure 3: Maturity in sGBM.
Mean squared error of the prediction of the AAPL option price with datataken on 28th February 2018 as a function of α for various maturity periods T (measured in years). quently present in nature, this will not only extend the framework of gGBM in analysing financialdata, but will also provide an avenue for applying the model in other scientific domains. Anotherfruitful research direction would be to incorporate the properties of gGBM in a wider frameworkfor financial modeling which includes the concept of “rough volatility”, where the instantaneousvolatility is driven by a (rough) fractional Brownian motion [62]. Building an explanatory modelfor the volatility in terms of gGBM would bring novel insights about the theoretical and empir-ical characteristics of the asset prices. We also leave for future analysis the problem of gGBMwith stochastic volatility, which can be treated in the framework of the Fokker-Planck equationfor gGBM with time varying volatility σ ( t ) , in analogy of the diffusing-diffusivity models forheterogeneous media [63–67]. Ackdnowledgments
The Authors acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG). TS wassupported by the Alexander von Humboldt Foundation. TS acknowledges Dr. Andrey Cherstvy22or the fruitful discussions and suggestions.
Abbreviations
The following abbreviations are used in this manuscript:GBM Geometric Brownian motionsGBM Subdiffusive geometric Brownian motiongGBM Generalised geometric Brownian motionBS Black-ScholesCTRW Continuous time random walkMSD Mean squared displacementML Mittag-LefflerTSLA TeslaAAPL Apple
A Solution of the Fokker-Planck equation for standard GBM
The solution of Eq. (3) can be found by using the Laplace-Mellin transform method [68]. TheLaplace transform is defined byˆ F ( s ) = L { f ( t ) } ( s ) = Z ∞ f ( t ) e − st dt , while the Mellin transform as [69]˜ F ( q ) = M { f ( x ) } ( q ) = Z ∞ x q − f ( x ) dx . The inverse Mellin transform then reads f ( x ) = M − (cid:8) ˜ F ( q ) (cid:9) ( x ) = π ı Z c + ı ∞ c − ı ∞ x − q ˜ F ( q ) dq . Therefore, by performing Laplace transform in respect to t and Mellin transform in respect to x inEq. (3), we have ˜ˆ F ( q , s ) = x q − × s − h σ ( q − )( q − ) + µ ( q − ) i , (A1)where we use M { δ ( x − x ) } ( q ) = x q − . Then the inverse Laplace transform yields˜ F ( q , t ) = x q − × exp σ (cid:20) q + (cid:18) µσ − (cid:19)(cid:21) t − (cid:16) µ − σ (cid:17) σ t , (A2)23here we use L − (cid:8) s − a (cid:9) ( t ) = e at . Applying the inverse Mellin transform and looking for thesolution in the form of the convolution integral of two functions [69], M { h ( x ) } ( q ) = ˜ H ( q ) and M { g ( x ) } ( q ) = ˜ G ( q ) , M − (cid:8) ˜ H ( q ) ˜ G ( q ) (cid:9) ( x ) = Z ∞ h ( r ) g ( x / r ) drr , we obtain the solution of the Fokker-Planck equation for GBM f ( x , t ) = Z ∞ δ ( r − x ) × exp − h log xr − (cid:16) µ − σ (cid:17) t i σ t ! ( x / r ) √ πσ t drr = x √ πσ t × exp − h log xx − (cid:16) µ − σ (cid:17) t i σ t . (A3)Here we use that h ( x ) = M − n x q − o ( x ) = δ ( x − x ) and g ( x ) = M − exp σ (cid:20) q + (cid:18) µσ − (cid:19)(cid:21) t − (cid:16) µ − σ (cid:17) σ t ( x )= x √ πσ t × exp − h log x − (cid:16) µ − σ (cid:17) t i σ t . We also used the properties of the inverse Mellin transform [69], M − { f ( q + a ) } ( x ) = x a M − { f ( q ) } and M − (cid:8) exp (cid:0) α q (cid:1)(cid:9) ( x ) = √ πα e − x α . Therefore, from the solution (A3) we conclude that the solution of the Fokker-Planck equation isa log-normal distribution.The n th moment h x n ( t ) i = R ∞ x n P ( x , t ) dx of the solution of Eq. (3) can be obtained by multi-plying the both sides of the equation with x n and integration over x . Thus, one has ∂∂ t h x n ( t ) i = (cid:20) σ n ( n − ) + µ n (cid:21) h x n ( t ) i , (A4)from where the n th moment becomes h x n ( t ) i = h x n ( ) i e ( σ n ( n − ) / + µ n ) t . (A5)24or n = h x ( t ) i =
1. The mean value ( n =
1) and the MSD have exponential dependence on time, h x ( t ) i = h x ( ) i e µ t and h x ( t ) i = h x ( ) i e ( σ + µ ) t , respectively, and thus, the variance becomes h x ( t ) i − h x ( t ) i = h x ( ) i e µ t (cid:16) e σ t − (cid:17) . (A6)The log-moments h log n x i = R ∞ log n xP ( x , t ) dx can be obtained by multiplying the both sidesof Eq. (3) with log n x and integration over x . Therefore, one finds the following equation (seeRef. [40] for details) ∂∂ t h log n x ( t ) i = (cid:18) µ − σ (cid:19) n h log n − x ( t ) i + σ n ( n − ) h log n − x ( t ) i . (A7)From here it follows that ∂∂ t h log x ( t ) i =
0, i.e., h log x ( t ) i = h log x ( ) i =
1. The case n = x ( t ) , ∂∂ t h log x ( t ) i = (cid:18) µ − σ (cid:19) h log x ( t ) i | {z } = (A8)i.e., h log x ( t ) i = h log x ( ) i + (cid:18) µ − σ (cid:19) t . (A9)For n = ∂∂ t h log x ( t ) i = (cid:18) µ − σ (cid:19) h log x ( t ) i + σ h log x ( t ) i (A10)which is given by h log x ( t ) i = h log x ( ) i + (cid:18) µ − σ (cid:19) t + (cid:18) µ − σ (cid:19) h log x ( ) i t + σ t . (A11)Therefore, for the log-variance one finds linear dependence on time h log x ( t ) i − h log x ( t ) i = σ t . (A12) B Derivation of the Fokker-Planck equation for gGBM fromCTRW theory
We use the approach given in Refs. [54, 70, 71]. Let us consider a CTRW for a particle at position x i which can move right to the position x i + = h x i or to left at position x i − = h x i , h >
0. Forthe CTRW on a geometric lattice we use h = + u , and at the end we will find the diffusion limit25 →
0. The probability density function (PDF) for the particle to jump to right is p r ( x , t ) , and forjump to left p l ( x , t ) . The total probability is p r ( x , t ) + p l ( x , t ) = λ ( x i , t , x j ) = p r ( x j , t ) δ ( x i − [ + u ] x j ) + p l ( x j , t ) δ ( x i − x j / [ + u ]) , and a waiting time PDF ψ ( t ) , related to the survival probability by φ ( t ) = − Z t ψ ( t ′ ) dt ′ , i.e., ˆ φ ( s ) = − ˆ ψ ( s ) s . By substitution in the master equation [54] ∂∂ t ρ ( x i , t ) = ∑ j λ ( x i , t , x j ) Z t K ( t − t ′ ) ρ ( x j , t ′ ) dt ′ − Z t K ( t − t ′ ) ρ ( x i , t ′ ) dt ′ , (B1)where K ( t ) = L − (cid:2) ˆ ψ ( s ) / ˆ φ ( s ) (cid:3) , one finds ∂∂ t ρ ( x i , t ) = p r (cid:18) x i + u , t (cid:19) Z t K ( t − t ′ ) ρ (cid:18) x i + u , t ′ (cid:19) dt ′ + p l ([ + u ] x i , t ) Z t K ( t − t ′ ) ρ (cid:0) [ + u ] x i , t ′ (cid:1) dt ′ − Z t K ( t − t ′ ) ρ ( x i , t ′ ) dt ′ . (B2)We consider generalised waiting time PDF, which in the Laplace space has the form [57, 72]ˆ ψ ( s ) = + τ η / ˆ η ( s ) , where τ η is a time parameter, which depends on η ( t ) . Therefore,ˆ φ ( s ) = τ η / ˆ η ( s ) s ( + τ η / ˆ η ( s )) , and ˆ K ( s ) = τ η s × ˆ η ( s ) , from where we find that R t K ( t − t ′ ) f ( t ′ ) dt ′ → τ η ddt R t η ( t − t ′ ) f ( t ′ ) dt ′ . From Eq. (B2) then weobtain ∂∂ t ρ ( x i , t ) = τ η p r (cid:18) x i + u , t ′ (cid:19) ∂∂ t Z t η ( t − t ′ ) ρ (cid:18) x i + u , t ′ (cid:19) dt ′ + τ η p l (cid:0) [ + u ] x i , t ′ (cid:1) ∂∂ t Z t η ( t − t ′ ) ρ (cid:0) [ + u ] x i , t ′ (cid:1) dt ′ − τ η ∂∂ t Z t η ( t − t ′ ) ρ ( x i , t ′ ) dt ′ . (B3)Let us now consider the diffusion limit ( u → τ η →
0) of Eq. (B3). From the normalisationcondition of the PDF ρ ( x , t ) given by ∑ i ρ ( x i , t ) = ∆ x i = u x i , one finds ∑ i ρ u ( x i , t ) ux i ∆ x i =
1, such that lim ∆ x i → ∑ i (cid:16) ρ u ( x i , t ) ux i (cid:17) ∆ x i = P u ( x , t ) = ρ u ( x , t ) / [ u x ] , one concludes that P ( x , t ) = lim u → P u ( x , t ) is normalised, i.e., R ∞ P ( x , t ) dx =
1. By introducing B u ( x i , t ) = p r ( x i , t ) − p l ( x i , t ) and b ( x , t ) = lim u → ∂∂ u B u ( x , t ) , inthe diffusion limit u → τ η →
0, where we assume that B ( x , t ) = lim u → B u ( x , t ) = ∂∂ t P ( x , t ) = D ∂∂ t Z t η ( t − t ′ ) ∂∂ x (cid:18) x ∂∂ x − k B T x F ( x ) (cid:19) P ( x , t ′ ) dt ′ , (B4)where D = lim u , τ η → u / [ τ η ] , F ( x ) = − V ′ ( x ) = k B T [ b ( x ) − ] / x , and P ( x , t ) ∝ exp (cid:16) − V ( x ) k B T (cid:17) is obtained from the long time steady state Boltzmann distribution [54]. For a logarithmic potential V ( x ) = v k B T log x , the force becomes F ( x ) = − k B T v / x . By using D = σ / v = − µ / D the Fokker-Planck equation becomes ∂∂ t P ( x , t ) = ∂∂ t Z t η ( t − t ′ ) ∂∂ x (cid:18) σ x ∂∂ x + [ σ − µ ] x (cid:19) P ( x , t ′ ) dt ′ , (B5)which can be rewritten in the form of Eq. (39). C General results for n th moment If we multiply both sides of Eq. (39) by x n , and integrate over x we find the n th moment h x n ( t ) i = R ∞ x n P ( x , t ) dx , ∂∂ t h x n ( t ) i = (cid:20) σ n ( n − ) + µ n (cid:21) µ ddt Z t η ( t − t ′ ) h x n ( t ′ ) i dt ′ , (C1)from where in the Laplace space it reads h ˆ x n ( s ) i = s − − ˆ η ( s ) h σ n ( n − ) + µ n i h x n ( ) i . (C2)From this result we obtain the normalization condition, ∂∂ t h x ( t ) i =
0, i.e., h x ( t ) i = h x ( ) i = n =
1, we find the equation for the mean value ∂∂ t h x ( t ) i = µ ddt Z t η ( t − t ′ ) (cid:10) x ( t ′ ) (cid:11) dt ′ , (C3)and its Laplace pair h ˆ x ( s ) i = s − − µ ˆ η ( s ) h x ( ) i . (C4)27n terms of the memory kernel γ ( t ) , Eq. (C4) reads h ˆ x ( s ) i = ˆ γ ( s ) s ˆ γ ( s ) − µ h x ( ) i . (C5)We note that for the standard case with η ( t ) = η ( s ) = / s ) we recover the previously obtainedresults for the GBM. For n = ∂∂ t h x ( t ) i = ( σ + µ ) ddt Z t η ( t − t ′ ) h x ( t ′ ) i dt ′ , (C6)and its Laplace pair h ˆ x ( s ) i = s − − ( σ + µ ) ˆ η ( s ) h x ( ) i , (C7)or h ˆ x ( s ) i = ˆ γ ( s ) s ˆ γ ( s ) − ( σ + µ ) h x ( ) i . (C8)We also calculate the log-moments h log n x ( t ) i = R ∞ log n x P ( x , t ) dx , which satisfy the followingintegral equation ∂∂ t h log n x ( t ) i = ∂∂ t Z t η ( t − t ′ ) (cid:20)(cid:18) µ − σ (cid:19) n h log n − x ( t ′ ) i + σ n ( n − ) h log n − x ( t ′ ) i (cid:21) dt ′ . (C9)Thus, we find that ∂∂ t h log x ( t ) i =
0, i.e., h log x ( t ) i = h log x ( ) i =
1. For the mean value ( n = ∂∂ t h log x ( t ) i = (cid:18) µ − σ (cid:19) ∂∂ t Z t η ( t − t ′ ) h log x ( t ′ ) i | {z } = dt ′ (C10)from where it follows h log x ( t ) i = h log x ( ) i + (cid:18) µ − σ (cid:19) Z t η ( t ′ ) dt ′ . (C11)For the expectation of the periodic log return with period ∆ t , we find1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i = (cid:18) µ − σ (cid:19) ∆ t Z t + ∆ tt η ( t ′ ) dt ′ = (cid:18) µ − σ (cid:19) I ( t + ∆ t ) − I ( t ) ∆ t ∼ ∆ t → (cid:18) µ − σ (cid:19) η ( t ) , (C12)where I ( t ) = R η ( t ) dt , i.e., I ′ ( t ) = η ( t ) . Therefore, the expectation of the periodic log returns28ehaves as the rate of the first log-moment,1 ∆ t h log ( x ( t + ∆ t ) / x ( t )) i ∼ ∆ t → ddt h log x ( t ) i . (C13)For n = ∂∂ t h log x ( t ) i = ∂∂ t Z t η ( t − t ′ ) (cid:20) (cid:18) µ − σ (cid:19) h log x ( t ′ ) i + σ h log x ( t ′ ) i (cid:21) dt ′ (C14)i.e., h log x ( t ) i = h log x ( ) i + Z t η ( t − t ′ ) (cid:26) (cid:18) µ − σ (cid:19) (cid:20) h log x ( ) i + (cid:18) µ − σ (cid:19) Z t ′ η ( t ′′ ) dt ′′ (cid:21) + σ (cid:27) dt ′ , (C15)and the log-variance becomes h log x ( t ) i − h log x ( t ) i = σ Z t η ( t ′ ) dt ′ + (cid:18) µ − σ (cid:19) " Z t η ( t − t ′ ) (cid:18) Z t ′ η ( t ′′ ) dt ′′ (cid:19) dt ′ − (cid:18) Z t η ( t ′ ) dt ′ (cid:19) . (C16) D Fox H -function The Fox H -function is defined by [73] H m , np , q (cid:20) z (cid:12)(cid:12)(cid:12)(cid:12) ( a , A ) , . . . , ( a p , A p )( b , B ) , . . . , ( b q , B q ) (cid:21) = H m , np , q (cid:20) z (cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) (cid:21) = π ı Z Ω θ ( s ) z s ds , (D1)where θ ( s ) is given by θ ( s ) = ∏ mj = Γ ( b j − B j s ) ∏ nj = Γ ( − a j + A j s ) ∏ qj = m + Γ ( − b j + B j s ) ∏ pj = n + Γ ( a j − A j s ) , 0 ≤ n ≤ p , 1 ≤ m ≤ q , a i , b j ∈ C , A i , B j ∈ R + , i = , ..., p , j = , ..., q . The contour Ω starting at c − ı ∞ and ending at c + ı ∞ separatesthe poles of the function Γ ( b j + B j s ) , j = , ..., m from those of the function Γ ( − a i − A i s ) , i = , ..., n . A special case of the Fox H -function is the exponential function [73], e − z = H , , (cid:20) z (cid:12)(cid:12)(cid:12)(cid:12) − ( , ) (cid:21) . (D2)The inverse Laplace transform of the Fox H -function reads [73] L − (cid:20) s − ρ H m , np , q (cid:20) a s σ (cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) (cid:21)(cid:21) ( t ) = t ρ − H m , np + , q (cid:20) at σ (cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p ) , ( ρ , σ )( b q , B q ) (cid:21) . (D3)29he Fox H -functions have the following property [73] z k H m , np , q (cid:20) z (cid:12)(cid:12)(cid:12)(cid:12) ( a p , A p )( b q , B q ) (cid:21) = H m , np , q (cid:20) z (cid:12)(cid:12)(cid:12)(cid:12) ( a p + kA p , A p )( b q + kB q , B q ) (cid:21) . (D4) References [1] Viktor Stojkoski, Zoran Utkovski, Lasko Basnarkov, and Ljupco Kocarev. Cooperation dy-namics in networked geometric brownian motion.
Physical Review E , 99(6):062312, 2019.[2] Viktor Stojkoski, Marko Karbevski, Zoran Utkovski, Lasko Basnarkov, and Ljupco Kocarev.Evolution of cooperation in populations with heterogeneous multiplicative resource dynam-ics. arXiv preprint arXiv:1912.09205 , 2019.[3] Ole Peters and William Klein. Ergodicity breaking in geometric brownian motion.
Physicalreview letters , 110(10):100603, 2013.[4] John Aitchison and James AC Brown.
The lognormal distribution with special reference toits uses in economics . Cambridge Univ. Press, 1957.[5] Sidney Redner. Random multiplicative processes: An elementary tutorial.
American Journalof Physics , 58(3):267–273, 1990.[6] Fischer Black and Myron Scholes. The pricing of options and corporate liabilities.
Journalof Political Economy , 81(3):637–654, 1973.[7] Robert C Merton. Optimum consumption and portfolio rules in a continuous-time model. In
Stochastic Optimization Models in Finance , pages 621–661. Elsevier, 1975.[8] Robert C Merton. Option pricing when underlying stock returns are discontinuous.
Journalof Financial Economics , 3(1-2):125–144, 1976.[9] Ole Peters. Optimal leverage from non-ergodicity.
Quantitative Finance , 11(11):1593–1602,2011.[10] Gleb Oshanin and Gregory Schehr. Two stock options at the races: Black–scholes forecasts.
Quantitative Finance , 12(9):1325–1333, 2012.[11] HGE Hentschel and Itamar Procaccia. Fractal nature of turbulence as manifested in turbulentdiffusion.
Physical Review A , 27(2):1266, 1983.[12] M Sc Mario Heidern¨atsch.
On the diffusion in inhomogeneous systems . PhD thesis, Technis-chen Universit¨at Chemnitz, 2015.[13] E Baskin and A Iomin. Superdiffusion on a comb structure.
Physical Review Letters ,93(12):120603, 2004.[14] Nassim Nicholas Taleb.
The Black Swan: The impact of the highly improbable , volume 2.Random House, 2007. 3015] John Cox. Notes on option pricing i: Constant elasticity of variance diffusions.
Unpublishednote, Stanford University, Graduate School of Business , 1975.[16] Steven L Heston. A closed-form solution for options with stochastic volatility with applica-tions to bond and currency options.
The review of financial studies , 6(2):327–343, 1993.[17] Patrick S Hagan, Deep Kumar, Andrew S Lesniewski, and Diana E Woodward. Managingsmile risk.
The Best of Wilmott , 1:249–296, 2002.[18] Andrew Matacz. Financial modeling and option theory with the truncated l´evy process.
In-ternational Journal of Theoretical and Applied Finance , 3(01):143–160, 2000.[19] Lisa Borland. A theory of non-gaussian option pricing.
Quantitative Finance , 2(6):415–431,2002.[20] Lisa Borland and Jean-Philippe Bouchaud. A non-gaussian option pricing model with skew.
Quantitative Finance , 4(5):499–514, 2004.[21] L Moriconi. Delta hedged option valuation with underlying non-gaussian returns.
PhysicaA: Statistical Mechanics and its Applications , 380:343–350, 2007.[22] Daniel T Cassidy, Michael J Hamp, and Rachid Ouyed. Pricing european options with a logstudent’s t-distribution: A gosset formula.
Physica A: Statistical Mechanics and its Applica-tions , 389(24):5736–5748, 2010.[23] Lasko Basnarkov, Viktor Stojkoski, Zoran Utkovski, and Ljupco Kocarev. Option pricingwith heavy-tailed distributions of logarithmic returns. arXiv preprint arXiv:1807.01756 ,2018.[24] Marcin Magdziarz. Black-scholes formula in subdiffusive regime.
Journal of StatisticalPhysics , 136(3):553–564, 2009.[25] CN Angstmann, BI Henry, and AV McGann. Time-fractional geometric brownian motionfrom continuous time random walks.
Physica A , 526:121002, 2019.[26] Grzegorz Krzy˙zanowski, Marcin Magdziarz, and Łukasz Płociniczak. A weighted finite dif-ference method for subdiffusive black–scholes model.
Computers & Mathematics with Ap-plications , 80(5):653–670, 2020.[27] Enrico Scalas, Rudolf Gorenflo, and Francesco Mainardi. Fractional calculus and continuous-time finance.
Physica A: Statistical Mechanics and its Applications , 284(1-4):376–384, 2000.[28] Marco Raberto, Enrico Scalas, and Francesco Mainardi. Waiting-times and returns in high-frequency financial data: an empirical study.
Physica A: Statistical Mechanics and its Appli-cations , 314(1-4):749–755, 2002.[29] Ole Peters and William Klein. Ergodicity breaking in geometric brownian motion.
PhysicalReview Letters , 110(10):100603, 2013. 3130] Viktor Stojkoski, Zoran Utkovski, Lasko Basnarkov, and Ljupco Kocarev. Cooperation dy-namics in networked geometric brownian motion.
Physical Review E , 99(6):062312, 2019.[31] Jun Wang, Jin-Rong Liang, Long-Jin Lv, Wei-Yuan Qiu, and Fu-Yao Ren. Continuous timeblack–scholes equation with transaction costs in subdiffusive fractional brownian motionregime.
Physica A , 391(3):750–759, 2012.[32] Gulnur Karipova and Marcin Magdziarz. Pricing of basket options in subdiffusive fractionalblack–scholes model.
Chaos, Solitons & Fractals , 102:245–253, 2017.[33] Janusz Gajda and Agnieszka Wyłoma´nska. Geometric brownian motion with tempered stablewaiting times.
Journal of Statistical Physics , 148(2):296–305, 2012.[34] N Leibovich and E Barkai. Infinite ergodic theory for heterogeneous diffusion processes.
Physical Review E , 99(4):042138, 2019.[35] Robert C Merton. Theory of rational option pricing.
The Bell Journal of economics andmanagement science , pages 141–183, 1973.[36] John C Hull.
Options futures and other derivatives . Pearson Education India, 2003.[37] Steven G Kou. A jump-diffusion model for option pricing.
Management Science , 48(8):1086–1101, 2002.[38] John Cox. Notes on option pricing i: Constant elasticity of variance diffusions.
Unpublishednote, Stanford University, Graduate School of Business , 1975.[39] Steven L Heston. A closed-form solution for options with stochastic volatility with applica-tions to bond and currency options.
The Review of Financial Studies , 6(2):327–343, 1993.[40] Francesco Mainardi.
Fractional calculus and waves in linear viscoelasticity: an introductionto mathematical models . World Scientific, 2010.[41] Hans C Fogedby. Langevin equations for continuous time l´evy flights.
Physical Review E ,50(2):1657, 1994.[42] Chao Li.
Option pricing with generalized continuous time random walk models . PhD thesis,Queen Mary University of London, 2016.[43] Willliam Feller.
An introduction to probability theory and its applications, vol 2 . John Wiley& Sons, 2008.[44] Marcin Magdziarz, Aleksander Weron, and Karina Weron. Fractional fokker-planck dynam-ics: Stochastic representation and computer simulation.
Physical Review E , 75(1):016708,2007.[45] Marcin Magdziarz, Aleksander Weron, and Joseph Klafter. Equivalence of the fractionalfokker-planck and subordinated langevin equations: the case of a time-dependent force.
Phys-ical review letters , 101(21):210601, 2008. 3246] Roberto Garra and Roberto Garrappa. The prabhakar or three parameter mittag–leffler func-tion: Theory and application.
Communications in Nonlinear Science and Numerical Simula-tion , 56:314–329, 2018.[47] Didier Sornette, Anders Johansen, and Jean-Philippe Bouchaud. Stock market crashes, pre-cursors and replicas.
Journal de Physique I , 6(1):167–175, 1996.[48] Antonio Mura, Murad S Taqqu, and Francesco Mainardi. Non-markovian diffusion equationsand processes: analysis and simulations.
Physica A: Statistical Mechanics and its Applica-tions , 387(21):5033–5064, 2008.[49] Antonio Mura and Gianni Pagnini. Characterizations and simulations of a class of stochasticprocesses to model anomalous diffusion.
Journal of Physics A: Mathematical and Theoreti-cal , 41(28):285003, 2008.[50] Vittoria Sposini, Aleksei V Chechkin, Flavio Seno, Gianni Pagnini, and Ralf Metzler. Ran-dom diffusivity from stochastic equations: comparison of two models for brownian yet non-gaussian diffusion.
New Journal of Physics , 20(4):043044, 2018.[51] Marcin Magdziarz and Janusz Gajda. Anomalous dynamics of black–scholes model time-changed by inverse subordinators.
Acta Physica Polonica B , 43(5), 2012.[52] Ralf Metzler and Joseph Klafter. The random walk’s guide to anomalous diffusion: a frac-tional dynamics approach.
Physics Reports , 339(1):1–77, 2000.[53] E Barkai. Fractional fokker-planck equation, solution, and application.
Physical Review E ,63(4):046118, 2001.[54] Mark M Meerschaert, David A Benson, Hans-Peter Scheffler, and Boris Baeumer. Stochasticsolution of space-time fractional diffusion equations.
Physical Review E , 65(4):041103, 2002.[55] Johannes HP Schulz, Eli Barkai, and Ralf Metzler. Aging renewal theory and application torandom walks.
Physical Review X , 4(1):011028, 2014.[56] Ren´e L Schilling, Renming Song, and Zoran Vondracek.
Bernstein functions: theory andapplications , volume 37. Walter de Gruyter, 2012.[57] Trifce Sandev, Ralf Metzler, and Aleksei Chechkin. From continuous time random walks tothe generalized diffusion equation.
Fractional Calculus and Applied Analysis , 21(1):10–28,2018.[58] Trifce Sandev, Igor M Sokolov, Ralf Metzler, and Aleksei Chechkin. Beyond monofractionalkinetics.
Chaos, Solitons & Fractals , 102:210–217, 2017.[59] Tilak Raj Prabhakar. A singular integral equation with a generalized mittag leffler functionin the kernel.
Yokohama Mathematical Journal , 19:7–15, 1971.[60] Trifce Sandev, Aleksei V Chechkin, Nickolay Korabel, Holger Kantz, Igor M Sokolov, andRalf Metzler. Distributed-order diffusion equations and multifractality: Models and solutions.
Physical Review E , 92(4):042117, 2015. 3361] Andrey G Cherstvy, Deepak Vinod, Erez Aghion, Aleksei V Chechkin, and Ralf Metzler.Time averaging, ageing and delay analysis of financial time series.
New Journal of Physics ,19(6):063045, 2017.[62] Omar El Euch.
Quantitative Finance under rough volatility . PhD thesis, Sorbonne universit´e,2018.[63] Rohit Jain and Kizhakeyil L Sebastian. Diffusion in a crowded, rearranging environment.
The Journal of Physical Chemistry B , 120(16):3988–3992, 2016.[64] Aleksei V Chechkin, Flavio Seno, Ralf Metzler, and Igor M Sokolov. Brownian yet non-gaussian diffusion: from superstatistics to subordination of diffusing diffusivities.
PhysicalReview X , 7(2):021002, 2017.[65] Vittoria Sposini, Aleksei V Chechkin, Flavio Seno, Gianni Pagnini, and Ralf Metzler. Ran-dom diffusivity from stochastic equations: comparison of two models for brownian yet non-gaussian diffusion.
New Journal of Physics , 20(4):043044, 2018.[66] Wei Wang, Andrey G Cherstvy, Xianbin Liu, and Ralf Metzler. Anomalous diffusion andnonergodicity for heterogeneous diffusion processes with fractional gaussian noise.
PhysicalReview E , 102(1):012146, 2020.[67] Wei Wang, Flavio Seno, Igor M Sokolov, Aleksei V Chechkin, and Ralf Metzler. Unexpectedcrossovers in correlated random-diffusivity processes.
New Journal of Physics , 22(8):083041,2020.[68] Trifce Sandev, Alexander Iomin, and Kocarev Ljupco. Hitting times in turbulent diffusiondue to multiplicative noise.
Physical Review E , 102(4):042109, 2020.[69] Fritz Oberhettinger.
Tables of Mellin transforms . Springer Science & Business Media, 2012.[70] Christopher N Angstmann, Isaac C Donnelly, Bruce Ian Henry, TAM Langlands, and PeterStraka. Generalized continuous time random walks, master equations, and fractional fokker–planck equations.
SIAM Journal on Applied Mathematics , 75(4):1445–1468, 2015.[71] Christopher N Angstmann, Isaac C Donnelly, and Bruce I Henry. Continuous time randomwalks with reactions forcing and trapping.
Mathematical Modelling of Natural Phenomena ,8(2):17–27, 2013.[72] Trifce Sandev, Aleksei Chechkin, Holger Kantz, and Ralf Metzler. Diffusion and fokker-planck-smoluchowski equations with generalized memory kernel.
Fractional Calculus andApplied Analysis , 18(4):1006, 2015.[73] Arakaparampil M Mathai, Ram Kishore Saxena, and Hans J Haubold.