Generalised Geometric Brownian Motion: Theory and Applications to Option Pricing

Classical option pricing schemes assume that the value of a financial asset follows a geometric Brownian motion (GBM). However, a growing body of studies suggest that a simple GBM trajectory is not an adequate representation for asset dynamics, due to irregularities found when comparing its properties with empirical distributions. As a solution, we investigate a generalisation of GBM where the introduction of a memory kernel critically determines the behaviour of the stochastic process. We find the general expressions for the moments, log-moments, and the expectation of the periodic log returns, and then obtain the corresponding probability density functions 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 the resulting generalised GBM (gGBM) in order to examine the empirical performance of a selected group of kernels in the pricing of European call options. Our results indicate that the performance of a kernel ultimately depends on the maturity of the option and its moneyness.


Introduction
Geometric Brownian motion (GBM) frequently features in mathematical modelling. The advantage of modelling through this process lies in its universality, as it represents an attractor of more complex models that exhibit non-ergodic dynamics [1][2][3]. As such, GBM has been used to underlie the dynamics of a diverse set of natural phenomena, including the distribution of incomes, body weights, rainfall, fragment sizes in rock crushing processes, etc. [4,5]. Nevertheless, perhaps the best-known application of GBM is in finance, and, in particular, in terms of the Black-Scholes (BS) model (or Black-Scholes-Merton model) [6][7][8] for the pricing of European options.
By construction, GBM is a simple continuous-time stochastic process in which the logarithm of the randomly varying quantity of interest follows a Brownian motion with drift. Its non-ergodicity is manifested in the difference between the growth rate that was observed in an individual trajectory and the ensemble average growth [9]. The time-averaged growth rate is dependent on both the drift and models that easily allow for such structural changes. We believe that the resolution to this issue lies in applying the concepts of time-averaging and ergodicity breaking to modelling financial time-series, and the gGBM framework offers a computationally inexpensive and efficiently tractable solution.
The paper is organised as follows. In Section 2, we provide an overview of GBM in the BS model and its use in option pricing. We also give detailed results for the so-called sGBM in terms of the fractional Fokker-Planck equation and its corresponding continuous time random walk (CTRW) model. In Section 3, we present gGBM and describe its properties using the subordination approach. In particular, we derive the corresponding Fokker-Planck equation with a memory kernel and obtain the respective moments and log-moments. The general function that is used in the Lévy exponent occurs as a memory kernel in the Fokker-Planck equation, which allows for us to recover the previously known results for GBM and sGBM. We consider generalisations of GBM and sGBM by introducing tempered sGBM, a mix of GBM and sGBM, as well as a mix of sGBMs. A numerical investigation of the properties of the model is given in Section 4 and an empirical example of application of the gGBM in option pricing is presented in Section 5. Section 6 summarises our findings. In the Appendices, we give detailed calculations as well as derivation of the Fokker-Planck equation for the gGBM within the CTRW theory.

Standard GBM
GBM has been applied in a variety of scientific fields [1,3,9,[38][39][40]. Mathematically, it is represented by the Langevin equation where x(t) is the particle position, µ is the drift, σ > 0 is the volatility, and B(t) represents a standard Brownian motion. The solution to Equation (1), in the Itô sense, is When the dynamics of the asset price follows a GBM, then a risk-neutral distribution (probability distribution that takes into account the risk of future price fluctuations) can be easily found by solving the corresponding Fokker-Planck equation to Equation (1), with initial condition f (x, t = 0) = δ(x − x 0 ). The solution of Equation (3) is the famed log-normal distribution whereμ = µ − σ 2 /2. We point out that this representation corresponds to the Itô interpretation of the multiplicative noise. There are also Stratonovich and Klimontovich-Hänggi interpretations, for which the corresponding Fokker-Planck equations are slightly different, see Refs. [12,41]. In finance math literature, the Itô convention is the standard interpretation. From the solution, it follows that the mean value and mean squared displacement (MSD) have exponential dependence on time, and x 2 (t) = x 2 0 e (σ 2 +2 µ)t , (6) respectively. Thus, the variance becomes The exact derivation of the GBM distribution and its moments is given in Appendix A. In the same way, one calculates the third and fourth moments, which are given by x 4 (t) = x 4 0 e (6σ 2 +4µ)t , respectively. The third and fourth moments are used to estimate the skewness g, and, respectively, excess kurtosis κ of the probability distribution of a random variable y, The skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable around its first moment, whereas the excess kurtosis evaluates the "tailedness" of the probability distribution. From these relations for the random variable x, one finds the skewness and excess kurtosis in the form Evidently, in GBM, the diffusion coefficient scales proportionally with the square of the position of the particle, i.e., D(x) = σ 2 x 2 /2, and, thus, the MSD has an exponential dependence on time. A more convenient measure instead of the MSD for geometric processes is the behaviour of the expectation of the logarithm of x(t), which, in asset pricing terms, represents the continuously compounded return of the asset. In the case of GBM, the expectation of the logarithm of the particle position has a linear dependence on time. This can be shown by calculation of the log-moments log n x(t) = ∞ 0 log n x P(x, t) dx, see Appendix A, Equation (A7). The mean value of the logarithm of x(t) becomes from where for the expectation of the periodic log return with period ∆t, one finds The second log-moment is given by which for the log-variance yields From Equations (A7), (14), and (16), we find the third and fourth log-moments, respectively. The logarithm of the process in GBM has both skewness and excess kurtosis of 0, which can be shown by using y → log x in Equations (10) and (11), and the previous results for the first four log-moments. This implies that there is no asymmetry and excess "tailedness" in GBM. However, real world return distributions are known to exhibit both positive asymmetry (g > 0) and fat-tailedness, i.e., positive excess kurtosis (κ > 0).

Black-Scholes Formula
As previously said, perhaps the best-known application of GBM is in finance and, in particular, the BS model for pricing of European options. Formally, a European option is a contract that gives the buyer (the owner or holder of the option) the right, but not the obligation, to buy, or sell an underlying asset or instrument x(T) at a specified strike price K on a specified date T. The seller has the corresponding obligation to fulfil 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 price is referred to as a call; an option that conveys the right of the owner to sell at a specific price is referred to as a put. Here, we are going to consider the valuation of call options, denoted as C(x, t), with the note that the derived results easily extend to put options.
In the modelling of financial assets, a standard assumption is that there is a risk-neutral distribution f (x, t) for the price of the asset. This measure is simply a probability distribution that takes into account the risk of future price fluctuations. Once a risk-neutral distribution is assigned, the value of the option is obtained by discounting the expectation of its value at the maturity T with respect to that distribution [6,42], i.e., where r is the risk-free rate of return and x 0 is the asset price at the beginning (t = 0). Notice that the integral is only calculated for the region of prices where the option has positive value, since, for asset price less than K, the option would not be exercised (i.e., its value is 0). Assuming that the asset price follows GBM dynamics, Equations (20) and (4) can be combined in order to derive an analytical formula for the value of the call option in the BS model. In particular, the European option C BS (x, t) (20) is a solution of the Black-Scholes equation, see, for example, [43], with initial condition C BS (x, T) = max{x − K, 0}, x ≥ 0, and boundary conditions C BS (x = 0, t) = 0, t ≥ T, and C BS (x → ∞, t) → x. By using t = 0 and T → t, one finds the equation with initial condition C BS (x, t = 0) = max{x − K, 0}, x ≥ 0, and boundary conditions The solution is where N(x) = 1 √ 2π x −∞ e −u 2 /2 du is the cumulative distribution function of the Gaussian distribution with zero mean and unit variance. Put simply, the two terms in the BS formula describe the current price of the asset weighted by the probability that the investor will exercise its option at time t and the discounted price of the strike price weighted by its exercise probability. The terms d 1,2 can be seen as measures of the moneyness of the option and N(d 1,2 ) as probabilities that the option will expire, while its value is in the money. The neat BS formulation has allowed for the model to be widely applied in both theoretical investigations and empirical implementations. However, the BS model has failed to adequately reproduce a plethora of real world properties. 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 a plethora of cases. In order to deal with this problem, extensions of the BS model have emerged, which include a combination of the GBM with jumps [8,44] or with stochastic volatility [15,16].

Subdiffusive GBM
One of the reasons why the standard GBM is not able to explain empirical data is because it fails to reproduce periods of constant prices that appear on markets with low number of transactions. The price in these constant periods can be described as a trapped particle, which, in physical systems, manifests anomalous diffusion (subdiffusion) [35,36]. To deal with this problem, the so-called subdiffusive GBM (sGBM) has been developed by Magdziarz [28], by using the subordination approach. The corresponding equation for the sGBM becomes the following fractional Fokker-Planck equation [28] (see also [29]) where is the Riemann-Liouville fractional derivative of order 0 < ν < 1 [45]. The Laplace transform of the Riemann-Liouville fractional derivative of a given function reads In order to avoid the somewhat unusual fractional initial condition, alternatively, we could alternatively use the integral version of the equation [46] f where . In Ref. [29], the time fractional Fokker-Planck Equation (26) for sGBM is derived within the CTRW theory for a particle on a geometric lattice in the presence of a logarithmic potential.
Here we note that the fractional Fokker-Planck Equation (26) can be obtained using the Langevin equation approach [47], i.e., by considering a CTRW model that was described by a coupled Langevin equations [48], Therefore, x(t) is parametrised in terms of the number of steps u, and the connection to the physical time t is given by T (u) = u 0 τ(u ) du, where τ(u) is a total of individual waiting times τ for each step. In mathematical terms, this is called subordination [49][50][51]. The noise ξ(u) is a white noise with zero mean and correlation ξ(u)ξ(u ) = 2δ(u − u ), while ζ(u) is one-sided α-stable Lévy noise with the stable index 0 < α < 1. The inverse process S(t) of the one-sided α-stable Levy process T (u) with a characteristic function e −sT (u) = e −s α u is given by S(t) = inf {u > 0 : T (u) > t}, i.e., it represents a collection of first passage times [47]. The CTRW is defined by the subordinated process X (t) = x(S (t)).
The PDF h(u, t) of the inverse process S(t) can be found from the relation [47] h(u, t) = − ∂ ∂u where Θ(z) is the Heaviside theta function. The Laplace transform then yieldŝ dt, from where one can easily arrive to the fractional Fokker-Planck Equation (26).
The mean value for sGBM is given by [29,48] where E α (z) is the one parameter Mittag-Leffler (ML) function [35,45] with (z ∈ C; (α) > 0), and Γ(·) is the Gamma function. The ML function is a generalisation of the exponential function, since E 1 (z) = e z . The Laplace transform of the one parameter ML function reads L {E α (at α )} (s) = s α−1 s α −a . The asymptotic behaviour of the mean is given by For the short time limit, we use the first two terms from the series expansion of the ML function (34), while, for the long time limit, we apply its asymptotic expansion formula E α (z) ∼ 1 α e z 1/α , z 1 [45,52]. Here, we note that the asymptotic behaviour of the ML function with negative argument has a power-law form, i.e., E α (−z α ) ∼ z −α Γ(1−α) for z 1 and 0 < α < 2 [45,52]. The MSD also is given through the one parameter ML function [29,48] From here, one concludes that sGBM is an exponentially fast process. Moreover, the third and fourth moments, respectively, become The first log-moment has the form [29] log which gives a power-law dependence with respect to time of the expectation of the log return with period ∆t, i.e., [29] 1 ∆t Such models have been used, for example, to explain the dynamics of an asset before a market crash [53]. The second log-moment becomes [29] from where, for the log-variance, one finds [29] which, in the long time, scales as t 2α (0 < α < 1), contrary to the linear scaling t for regular GBM (α = 1). For the third and fourth log-moments, we obtain

Generalised GBM
In this section, we consider a generalisation of GBM, under which the standard and subdiffusive GBM arise as special cases, by using the subordination approach. Here, we present analytical expressions for the first four moments and log-moments of the process for a variety of special cases. They are thoroughly analysed in the numerical experiments section.
The continuous time random walk approach to the corresponding Fokker-Planck equation is given in Appendix B in detail. The same Fokker-Planck equation can be obtained using the coupled Langevin equations approach [47], as given in Equations (29) and (30), where the waiting times are given by e −sT (u) = e −uΨ(s) , withΨ(s) = 1/η(s).

Subordination Approach
The generalisation of GBM that we consider is the model introduced by Magdziarz and Gajda [37] in the form of a stochastic process and T (u) is an infinite divisible process, i.e., a strictly increasing Lévy motion with [37] e −sT (u) = e −uΨ(s) , andΨ(s) is the Lévy exponent [28,37,39]. Here we considerΨ(s) = 1/η(s). The current process should not be confused with the generalised grey Brownian motion, see Ref. [54][55][56].
Next we find the PDF of gGBM which subordinates the processes from the time scale t (physical time) to the GBM on a time scale u (operational time). Specifically, the PDF P(x, t) of a given random process X (t) can be represented as [28,46,[57][58][59] 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 Laplace space, Equation (46) readŝ we then haveP By Laplace transform of the Fokker-Planck Equation (3) for the GBM, and using relation (49), one finds that the PDF P(x, s) satisfies After applying an inverse Laplace transform, we arrive at the generalised Fokker-Planck equation (see Refs. [37,48], where one-sided α-stable waiting times are considered in detail) where η(t) is a so-called memory kernel. One observes that, for η(t) = 1, we arrive at the Fokker-Planck Equation (3) for the GBM, and for η(t) = t α−1 Γ(α) at the time fractional Fokker-Planck Equation (26) for the sGBM. From Equations (47) and (48), we find for the PDF in the Laplace domain, see also Ref. [48],

Remark 1.
Here, we note that there are restrictions on the choice of the memory kernel η(t) since the PDF (46) should be non-negative. From the subordination integral it follows that the subordination function h(u, t) should be non-negative, which, according to the Bernstein theorem, means that its Laplace transform (48) should be a completely monotone function [60]. Therefore, the PDF (46) will be non-negative if 1/[sη(s)] is a completely monotone function, and 1/η(s) is a Bernstein function, see Refs. [61,62].

Remark 2.
We note that Equation (50) can be written in an equivalent form as where the memory kernel γ(t) is connected to η(t) in Laplace space as γ(s) = 1/[sη(s)], see Ref. [61]. From this relation, we find that for GBM (η(t) = 1, i.e.,η(s) = 1/s) the memory kernel γ(t) is given by , and, thus, Equation (53) reads where is the Caputo fractional derivative of order 0 < ν < 1 [45]. The Laplace transform of the Caputo derivative of a given function reads . We note that, with the appropriate restrictions for η(t) and γ(t), both formulations are equivalent.
, gGBM corresponds to sGBM. Using the subordination approach one finds [28] where H m,n p,q (z) is the Fox H-function, see Appendix D. By inverse Laplace transform (A42), one obtains [28] h where we applied property (A43). The solution in Laplace space then becomeŝ which is obtained in Ref. [48] in a similar way. From here, we can plot the PDF using numerical inverse Laplace transform techniques.

Generalised BS Formula
If we consider that the asset price follows a gGBM, then the generalised BS (gBS) formula for the option price is [37] where C BS (x, t) is taken from the BS Formula (25), and h(x, T) is the subordination function defined by Equation (48) in the Laplace domain. By Laplace transform one findŝ Therefore, from Equation (22), the corresponding equation for the option price becomes [48]

Calculation of Moments
The nth moment X n (t) = ∞ 0 x n P(x, t) dx can be calculated by multiplying both sides of Equation (51) by x n and integration over x, see Appendix C. In the Laplace domain, this results in From this result, we reproduce the normalisation condition x 0 (t) = x 0 0 = 1. The general results for the first four moments in terms of the memory kernel become The log-moments log n x(t) = ∞ 0 log n x P(x, t) dx, can also be calculated exactly through the memory kernel, see Appendix C. The normalisation condition is satisfied, i.e., log 0 x(t) = 1, while the log-mean reads From here, we find, for the expectation of the periodic log return with period ∆t where I(t) = η(t) dt, i.e., I (t) = η(t). Therefore, the expectation of the periodic log returns behaves as the rate of the first log-moment, which is proportional to the memory kernel η(t). Moreover, for the second log-moment, we find from where the log-variance becomes The general results for third and fourth moments are given in Appendix C. From all of these general formulas, one can easily recover the previous results for the standard GBM (η(t) = 1, i.e.,η(s) = 1/s) and sGBM (η(t) = t α−1 /Γ(α), i.e.,η(s) = s −α , 0 < α < 1).

Exponentially Truncated Subdiffusive GBM
As an example for another memory kernel in gGBM, we consider a power-law memory kernel with exponential truncation, where τ is a characteristic crossover time scale, 0 < α < 1. Such forms are important in many real-world applications, in which the scale-free nature of the waiting time dynamics is broken at macroscopic times t τ [61]. Therefore,η where we use the shift rule of the Laplace transform, L e −at f (t) =F(s + a), forF(s) = L { f (t)}.
The mean value reads, and the MSD is The third and fourth moments become respectively. From the general result for the log-mean, we find that where (z, β ∈ C; (α) > 0). Here, γ(a, z) = z 0 t a−1 e −t dt = Γ(a)e −z z a E 1,a+1 (z) is the incomplete gamma function, and is the two parameter ML function [45]. Note that, the Laplace transform of the two parameter ML function reads L t β−1 E α,β (at α ) (s) = s α−β s α −a . The asymptotic expansion formula for the two parameter ML function is E α,β (z) ∼ 1 α e z 1/α z (1−β)/α , z 1 [45,52], while the asymptotic behaviour for negative arguments is given by power-law decay, E α,β (−z α ) ∼ z −α Γ(β−α) , z 1 [45,52]. For the expectation of the periodic log return with period ∆t, we find This leads to a long run log return of 0, whereas on the short time scale the same observable behaves in the same way as sGBM. As such, the model can be used in order to model early herd behaviour, where the price of an asset grows simply as a consequence of investors following trends (short run behaviour), which last until the trade of the asset becomes congested (long run behaviour). The second log-moment is from where the log-variance becomes Here, we note that, for t/τ 1, the obtained results correspond to those that were obtained for sGBM, as it should be since the exponential truncation has no influence on the process. We observe that on the long run the log-variance becomes constant, i.e., it is equal to σ 2 τ α +μ τ 2α . In a similar way, for the third and fourth log-moments, we find log 4 x(t) = log 4 x 0 + 2 2μ log 3 x 0 + 3σ 2 log 2 x 0 τ α γ(α, t/τ) Γ(α) The subordination function, in this case, is given bŷ from where one can analyse the PDF P(x, t).

Combined Standard and Subdiffusive GBM
As another application, let us consider the combination of GBM and sGBM, represented by the memory kernel where 0 < α < 1, w 1 + w 2 = 1, andη This case combines both motions governed by Equations (3) and (26). In this case, in a jump picture, normal GBM steps occur with weight w 2 , while power-law waiting time steps are realised with weight w 1 .
The mean value for this case is given by is the three parameter ML function [63], and (γ) n = Γ(γ + n)/Γ(γ) is the Pochhammer symbol. The Laplace transform of the three parameter ML function reads L t β−1 E γ α,β (at α ) (s) = s αγ−β (s α −a) γ . From here, we see that, for w 1 = 0 and w 2 = 1, only the term for n = 0 in Equation (88) survives, which yields the result for standard GBM as it should be. The opposite case, with w 1 = 1 and w 2 = 0, yields as it should be for the sGBM. For the second, third, and fourth moments, we find Following the same procedure as previously, for the log-mean we find log x(t) = log x 0 +μ w 1 t α Γ(α + 1) and for the expectation of the periodic log return with period ∆t, This model introduces subdiffusive and trapping asset dynamics on short time scales (i.e., then the part multiplied with w 1 is much bigger), whereas, on the long run, we recover the standard GBM dynamics. The second log-moment yields from where the log-variance becomes Similarly to the behaviour of the first log moment, in the log variance, for short time scales, the sGBM dynamics dominates. However, we observe that, on the long run, the dynamics is a combination of the two kernels, since the dominant term is w 1 w 2 t α+1 . Moreover, the third and fourth log-moments read log 4 x(t) = log 4 x 0 + 2 2μ log 3 x 0 + 3σ 2 log 2 x 0 w 1 t α Γ(α + 1) The subordination function for this case is given bŷ where the Lévy exponent isΨ(s) = w 1 s −1 + w 2 s −α −1 .

Mix of Subdiffusive GBMs
We may further analyse the case of a mix of two sGBM with different power-law memory functions, where 0 < α 1 < α 2 < 1, w 1 + w 2 = 1, and This situation corresponds to the case of two different groups of periods of constant prices. For physical systems, this situation means that the process is a sum of two random processes with different waiting times [64], as represented by the memory kernel (101). Therefore, for the mean, we find while, for the second, third, and fourth moments, respectively, we obtain x 4 (t) = x 4 0 ∞ ∑ n=0 w n 1 (6σ 2 + 4µ) n t α 1 n E n+1 α 2 ,α 1 n+1 w 2 (6σ 2 + 4µ)t α 2 .
Similarly, the log-mean yields The expectation of the log return with period ∆t, then becomes Because 0 < α 1 < α 2 < 1, on short times, the part of first sGBM dominates, whereas, on long times, it is the characteristic of the second sGBM that determines the dynamics. The second log-moment becomes and for the log-variance we find In this case, for short times, the kernel with the smaller exponent dominates the variance. Interestingly, for long times, this observable is determined by the magnitude of the larger exponent, which is opposite from the previous kernel examples. Moreover, the third and fourth log-moments read For the mix of subdiffusive GBMs, the subordination function becomeŝ where the Lévy exponent isΨ(s) = [w 1 s −α 1 + w 2 s −α 2 ] −1 , see Refs. [64,65]. Figure 1a provides an intuitive illustration of the gGBM dynamics under various choices for the kernel. As argued, for standard GBM we observe smooth dynamics without periods of constant prices, whereas there is more turbulence in the asset price dynamics in the gGBM case. The periods of constant prices that are 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 longer periods of constant prices. In Figure 1b,c we plot, respectively, the numerical approximations for the first moment and the MSD for GBM, sGBM, a mix of GBM and sGBM, and a mix of sGBMs. One can easily notice the nonlinear behaviour in the generalisations of GBM. For long times, all gGBMs give an exponential dependence of the first moment and the MSD on time but with smaller slope than the one of GBM. Finally, Figure 1d shows the empirical PDF for the logarithmic return at t = 1 year. For each of the studied generalisations of GBM, the PDF is characterised with fatter tails (which should increase as the α parameters decrease), which means that it is more prone to producing values that fall far from the average. The fat-tailed property can be observed in greater detail in Figure 2, where we plot the skewness g and excess kurtosis κ in gGBM as a function of α (for sGBM and the mix of GBM-sGBM) or α 1 (for the mix of sGBM) for the logarithmic return at t = 1 year. As is the case with general fat-tailed distributions, all of the generalisations exhibit positive skewness and excess kurtosis. This is exactly what makes the gGBM framework useful in understanding the statistical behaviour of the asset price dynamics. Moreover, the figure indicates that there is an inverse relationship between these two statistics and α (α 1 ), i.e., as α (α 1 ) increases, g and κ decrease. For small α (α 1 ), sGBM is the model with the largest skewness and excess kurtosis, followed by the mix of sGBM and the mix of GBM-sGBM. For large α (α 1 ), the mix of GBM-sGBM remains the model with the lowest g and κ, but, now, the mix of sGBM becomes the model with the largest values of the two statistics. The rationale for this observation is that the GBM-sGBM model already includes a GBM term, which greatly induces the fat tails, whereas the operational time in the mix of sGBM becomes dominated by the value of the second sGBM, which has a fixed value for α 2 , thereby resulting in larger skewness and excess kurtosis. Finally, we investigate the dependence of two standard quantities that are relevant in an option pricing scheme: (i) at-the-money (ATM) implied volatility and (ii) ATM volatility skew with respect to the gGBM parameters. Both of the quantities are a part of the moneyness property of the option. Moneyness describes the relative position of the current price of the asset (x 0 ) with respect to the strike price of the option (K). An option whose strike price is equal to the current price of the asset is said to be at-the-money; if the strike price is larger than the current price, the option is "out of the money"; and, if the strike price is smaller than the current price, the option is described to be "in the money".

Numerical Experiments
Moneyness is usually examined by plotting the implied volatility of an option for an underlying asset as a function of its strike price. This plot is formally known as the volatility smile. Its name is derived from the usual pattern, which suggests that the implied volatility is the lowest for options that are ATM, i.e., the plot looks like a parabola (smile). The ATM volatility skew is the first derivative of the volatility smile at the ATM point [66]. A larger volatility skew implies that the implied volatility increases faster for options that are near the money (options whose strike price is near the current asset price), and vice versa. Figure 3 displays the functional relationship between the two studied quantities and the parameter α (for sGBM and the mix of GBM-sGBM) or α 1 (for the mix of sGBMs). For each gGBM model, we get that lower ATM implied volatilities can be attributed to lower diffusion rates, thereby suggesting that the gGBM framework can be used to describe the empirical at-the-money observations. Identically, there is a direct relationship between ATM volatility skews and α (α 1 ). In other words, in the gGBM framework, the slope of the volatility smile for at-the-money options behaves in the same manner as the magnitude of the implied volatility. This leads to lower near-the-money implied volatilities. We assume that r = 0.02. Moreover, for the mix GBM-sGBM case, we set w 1 = w 2 = 0.5, and for the mix of sGBM case α 2 = 0.8 and w 1 = w 2 = 0.5.

Empirical Example
In order to illustrate the power of the gGBM framework in the evaluation of options, we utilise empirical data of American options for two companies: Tesla (TSLA) and Apple (AAPL). By definition, the dynamics of American options differ from European, as they allow for exercising of the option at any time before the option expires. Nevertheless, as given in Ref. [67], one can rely on the fact that American options on non-dividend-paying stocks have the same value as their European counterpart. This relation has allowed for the empirical examination of a pricing scheme of European options to be widely done via data for American ones.
For our analysis, we use freely available data from the Nasdaq's Options Trading Centre. This dataset offers daily data for options of all companies quoted on the NASDAQ stock market. However, the options for most companies have small sample size. Therefore, we restricted the empirical analysis to Tesla and Apple, whose options are more frequently traded. In our estimations, the risk-free rate of return r is taken simply as the three-Month Treasury Bill Secondary Market Rate at the date of observation. The volatility parameter σ, on the other hand, was inferred from the values of the options on the market as the value that produces the minimum root mean squared pricing error (RMSE) in their fit.
For both assets, we evaluate the predictive performance of the same sGBM, a mix of GBM-sGBM, and a mix of sGBM models that were used in the numerical analysis of the previous section. The predicted option prices by the model were inferred via a Monte Carlo estimation of (59) (see Refs [28,48]).

Moneyness in TSLA:
Let us now turn our attention to Figure 4, where we use TSLA data gathered on 1st March 2018 on options which expire on 16th March 2018 to examine the dependence of the sGBM, the mix of GBM-sGBM and the mix of sGBM models on α in predicting the option price. We discover that, in general, the TSLA asset price dynamics are best described as a sGBM, and the minimal error occurs around α = 0.2. For large α, the mix of sGBM becomes the best performer, because it inherently includes a process with a lower subdiffusion and, thus, it has a close resemblance to the sGBM process with low α than the other models. For every α, the mix of GBM-sGBM has the worst performance since it includes a term with a normal diffusion. In the inset plot of Figure 4, we provide a more detailed information on the predictive properties of sGBM by examining its relation with the moneyness of the option. In it, we vary the diffusion parameter α, and plot the absolute difference in the estimated option price C g and the observed option price as 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 the TSLA price, a transition occurs and α = 1 becomes the worst predictor of the option price, whereas the lower the subdiffusion parameter, the better prediction we obtain. For options that are out of the money, it appears that the performance of the prediction for the option price does not depend on α. These findings are in agreement with the discussion presented in the previous section regarding the ATM implied volatility and ATM volatility skew.
Maturity in Apple (AAPL): Next, we use AAPL data that were gathered for at-the-money options on 28th February 2018 and examine how the maturity T affects the performance of the same models in predicting the option price. We focus solely on data for at-the-money options in order to remove the potential bias in the prediction error that arises from the potential moneyness effect, which as we saw in the above paragraph might arise (i.e., we discovered that the error rate with respect to α (α 1 ) depends on the moneyness of the option).
For this purpose, in Table 1, we show the minimum RMSE and the corresponding optimal α (α 1 ) of the AAPL option price prediction for seven different maturity periods. We observe that for the shortest-term maturity (T = 0.006 years, or two days), a subdiffusive model (in particular sGBM or a mix of GBM-sGBM) is the model that best describes the option price. As the maturity increases, the optimal α (α 1 ) also increases, reaching a maximal value of 1 for options maturing at T = 0.101 years (one month) and T = 0.14 years (35 working days). However, for the options with the longest maturity, again a model with a very low subdiffusion rate is an optimal fit. Hence, the RMSE of the prediction depends on both the maturity of the options and choice of α (α 1 ). Table 1. Minimum prediction error and optimal α (α 1 ) for ATM AAPL options. For the mix GBM-sGBM case we set w 1 = w 2 = 0.5, and for the mix of sGBM case α 2 = 0.8 and w 1 = w 2 = 0.5. To provide a better depiction on the role of α (α 1 ), in the predictive performance of sGBM, a mix of GBM-sGBM and a mix of sGBM, in Figure 5 we depict the RMSE of the option price prediction as a function of the parameter α (α 1 ). We observe that, in general, only for the options with the longest maturity, one can observe a clear extremum, whereas, for every other maturity, it appears that every α (α 1 ) represents an adequate fit. Increasing this parameter will only lead to marginal and insignificant improvements in the error rate. As a consequence, one might even argue that different gGBM kernels can lead to similar outcomes in the pricing of options, an interesting finding as such. Evidently, the performance of a kernel ultimately depends on the physical properties of the option. On the first sight, this conclusion appears intuitive-obviously the known information for the properties of the asset greatly impacts its price, the observation that a slight change in the known information may drastically change the dynamics suggests that there is a need in the option pricing literature for models that easily allow for such structural changes. In this aspect, we believe that the generalised GBM approach offers a computationally inexpensive and efficiently tractable solution to this issue. Consequently, we stress that a significant improvement of the description of the data in the gGBM framework can be achieved with comparatively few additional parameters.

Conclusions
We investigated the potential of GBM extensions that are based on subdiffusion to model and predict the price of options. By assuming that the price of the asset underlying the option undergoes a subdiffusive 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 gGBM instance is critically determined by a memory kernel. The advantage of gGBM comes in the flavour of allowing various forms for the functional form of the kernel. Depending on its choice, we may end up with asset price dynamics whose behaviour significantly varies on the short time in comparison to its long run characteristics. This, in turn, may induce observations of the properties of the asset price that more closely mimic realistic behaviour than standard GBM.
We explored the ability of gGBM to fit and predict real option values. Our empirical analysis confirmed the characteristics of gGBM, as we discovered that the performance of a certain choice of memory kernel is uniquely determined by the parameters of the option, such as its maturity and its moneyness. Because 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 in the analysis of a single time-series (or a set of few), which is characterised with non-ergodic dynamics. The non-ergodicity creates non-equilibrium dynamics which, consequently, makes studies of the ensemble behaviour irrelevant. This leads to the introduction of novel strategies for analysing financial data [9,68].
In line with our conclusions, we believe that the next step in uncovering the properties of gGBM is demonstrating the ergodicity breaking of the process. Because multiplicative processes are frequently present in nature, this will not only extend the framework of gGBM in analysing financial data, but will also provide an avenue for applying the model in other scientific domains. Another fruitful research direction would be to incorporate the properties of gGBM in a wider framework for financial modelling, which includes the concept of "rough volatility", where the instantaneous volatility is driven by a (rough) fractional Brownian motion [69]. Building an explanatory model for the volatility in terms of gGBM would bring novel insights regarding the theoretical and empirical characteristics of the asset prices. We also leave for future analysis the problem of gGBM with stochastic volatility, which can be treated in the framework of the Fokker-Planck equation for gGBM with time varying volatility σ(t), in analogy of the diffusing-diffusivity models for heterogeneous media [56,[70][71][72][73].  we obtain the solution of the Fokker-Planck equation for GBM Here we used that We used the properties of the inverse Mellin transform [75], Therefore, from the solution (A3) we conclude that the solution of the Fokker-Planck equation is a log-normal distribution.
The nth moment x n (t) = ∞ 0 x n P(x, t) dx of the solution of Equation (3) can be obtained by multiplying the both sides of the equation with x n and integration over x. Thus, one has ∂ ∂t x n (t) = σ 2 2 n(n − 1) + µ n x n (t) , from where the nth moment becomes x n (t) = x n (0) e (σ 2 n(n−1)/2+µ n)t .
For n = 0 one observes that the solution of the Fokker-Planck equation for GBM is normalised, i.e., x 0 (t) = 1. The mean value (n = 1) and the MSD have exponential dependence on time, x(t) = x(0) e µ t and x 2 (t) = x 2 (0) e (σ 2 +2 µ)t , respectively, and thus, the variance becomes The log-moments log n x = ∞ 0 log n xP(x, t) dx can be obtained by multiplying the both sides of Equation (3) with log n x and integration over x. Therefore, one finds the following equation (see Ref. [45] for details)
For n = 2 we obtain the second log-moment ∂ ∂t log 2 x(t) = 2 µ − σ 2 2 log x(t) + σ 2 log 0 x(t) (A10) which is given by Therefore, for the log-variance one finds linear dependence on time

Appendix B. Derivation of the Fokker-Planck Equation for gGBM from CTRW Theory
We use the approach given in Refs. [58,76,77]. Let us consider a CTRW for a particle at position x i which can move right to the position x i+1 = h x i or to left at position x i−1 = 1 h x i , h > 0. For the CTRW on a geometric lattice we use h = 1 + u, and at the end we will find the diffusion limit u → 0. The probability density function (PDF) for the particle to jump to right is p r (x, t), and for jump to left p l (x, t). The total probability is p r (x, t) + p l (x, t) = 1.
We consider a multiplicative jump length PDF on a geometric lattice [58], and a waiting time PDF ψ(t), related to the survival probability by By substitution in the master equation [58] ∂ ∂t where K(t) = L −1 ψ (s)/φ(s) , one finds From this result we obtain the normalisation condition, ∂ ∂t x 0 (t) = 0, i.e., x 0 (t) = x 0 (0) = 1. For n = 1, we find the equation for the mean value ∂ ∂t and its Laplace pair In terms of the memory kernel γ(t), Equation (A21) reads We note that for the standard case with η(t) = 1 (η(s) = 1/s) we recover the previously obtained results for the GBM. For n = 2 we obtain the equation for the second moment, or the MSD, and its Laplace pair or x 2 (s) =γ (s) sγ(s) − (σ 2 + 2µ) x 2 (0) .
We also calculate the log-moments log n x(t) = ∞ 0 log n x P(x, t) dx, which satisfy the following integral equation