Next Article in Journal
The Meta-Dynamic Nature of Consciousness
Next Article in Special Issue
On (Non-)Monotonicity and Phase Diagram of Finitary Random Interlacement
Previous Article in Journal / Special Issue
A Continuous-Time Random Walk Extension of the Gillis Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Generalised Geometric Brownian Motion: Theory and Applications to Option Pricing

1
Faculty of Economics, Ss. Cyril and Methodius University, 1000 Skopje, Macedonia
2
Research Centre for Computer Science and Information Technologies, Macedonian Academy of Sciences and Arts, Bul. Krste Misirkov 2, 1000 Skopje, Macedonia
3
Institute of Physics & Astronomy, University of Potsdam, D-14776 Potsdam-Golm, Germany
4
Institute of Physics, Faculty of Natural Sciences and Mathematics, Ss. Cyril and Methodius University, Arhimedova 3, 1000 Skopje, Macedonia
5
Faculty of Computer Science and Engineering, Ss. Cyril and Methodius University, P.O. Box 393, 1000 Skopje, Macedonia
*
Author to whom correspondence should be addressed.
Submission received: 30 October 2020 / Revised: 11 December 2020 / Accepted: 16 December 2020 / Published: 18 December 2020
(This article belongs to the Special Issue New Trends in Random Walks)

Abstract

:
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.

1. 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 randomness in the system, whereas the ensemble growth rate is solely dependent on the drift. If only a single system is to be modelled, in the long run only the time-averaged growth rate is observed. This is naturally the case in financial market dynamics, for which only single time series exist, 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 diffusion, 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 a turbulent diffusion, the contaminant spreads very fast. For the case of Richardson diffusion, the position-dependent diffusion coefficient behaves as D ( x ) x 4 / 3 and the relative mean squared displacement (MSD) scales as x 2 ( t ) t 3 [11]. However, the fast spread of contaminants can be essentially increased due 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 of various characteristics observed in reality [14]. As a solution, several alternating theories have been proposed, among which are: stochastic volatility [15,16,17]; local volatility [18,19]; time-varying volatility [20,21], models utilising stochastic processes in which the noise follows a fat-tailed distribution [22,23,24,25,26,27]; and, generalisations of GBM based on subdiffusion [28,29,30]. In the first approach, the volatility is a stochastic process itself. In local volatility models, asset prices follow a stochastic differential equation whose diffusion coefficient is a function of the price. The time-varying volatility models assume that the latent volatility is predictable with respect to the information set. In contrast to stochastic volatility, in these models the conditional variance is a deterministic function of the model parameters and past data. Utilising stochastic processes in which the noise follows a fat-tailed distribution intuitively leads to the desired result. However, it has been acknowledged that several models that are described as a fat-tailed process can, in truth, be derived from stochastic and local volatility models [31,32,33,34]. The last, i.e., the subdiffusive approach, differently from the other views, assumes anomalous price dynamics. Concretely, the observation that the distribution of log returns is fat-tailed can be attributed to prolonged periods, in which the price of the asset exhibits approximately constant extreme values. These constant periods can be considered to be trapping of particles, as is done in physical systems that manifest anomalous diffusion (subdiffusion) [35,36]. While the resulting subdiffusive GBM (sGBM) is able to easily reproduce real-life properties, the literature lacks an extensive study in which the exact empirical characteristics of the subdiffusive model are presented.
The purpose of this paper is to propose a unifying framework for the application of subdiffusive GBM models in option pricing. We do this by providing a thorough investigation on the properties of the so-called generalised GBM (gGBM) [37]. gGBM is a stochastic process whose behaviour is critically determined by a memory kernel. By choosing the appropriate kernel, we recover the standard GBM and the typically used subdiffusive GBM models [28,29,30]. In order to understand the behaviour of gGBM under various kernels, we perform a detailed analysis of the moments, log-moments, and the expectation of the periodic log returns, and obtain the corresponding probability density functions by using the subordination approach. We show that the dynamics of the model can be easily adjusted in order to mimic periods of constant prices and/or fat-tailed observations of returns, thus corresponding to realistic scenarios. More importantly, it is known that gGBM leads to risk-neutral asset price dynamics, and, thus, it is adequate for their modelling [37]. We utilise this property of the model to investigate its capability to predict empirical option values. We find that the performance of a kernel ultimately depends on the parameters of the option, such as its maturity and moneyness. The first property describes the time that is left for the option to be exercised, whereas the second characteristic depicts the relative position of the current price with respect to the strike price of the option. At first sight, this conclusion appears intuitive—obviously the known information for the properties of the asset greatly impacts its price. However, 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. 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.

2. Background

2.1. 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
d x ( t ) = x ( t ) μ d t + σ d B ( t ) , x 0 = x ( 0 ) ,
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
x ( t ) = x 0 e ( μ σ 2 2 ) t + σ B ( t ) , x 0 = x ( 0 ) > 0 .
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),
t f ( x , t ) = μ x x f ( x , t ) + σ 2 2 2 x 2 x 2 f ( x , t ) ,
with initial condition f ( x , t = 0 ) = δ ( x x 0 ) . The solution of Equation (3) is the famed log-normal distribution
f ( x , t ) = 1 x 2 π σ 2 t × exp log x log x 0 μ ¯ t 2 2 σ 2 t .
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,
x ( t ) = x 0 e μ t ,
and
x 2 ( t ) = x 0 2 e ( σ 2 + 2 μ ) t ,
respectively. Thus, the variance becomes
x 2 ( t ) x ( t ) 2 = x 0 2 e 2 μ t e σ 2 t 1 .
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 3 ( t ) = x 0 3 e ( 3 σ 2 + 3 μ ) t ,
x 4 ( t ) = x 0 4 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,
g = ( y y ) 3 ( y 2 y 2 ) 3 / 2 ,
κ = ( y y ) 4 ( y 2 y 2 ) 2 3 .
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
g = e σ 2 t 1 e σ 2 t + 2 ,
κ = e 4 σ 2 t + 2 e 3 σ 2 t + 3 e 2 σ 2 t 6 .
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 ) d x , see Appendix A, Equation (A7). The mean value of the logarithm of x ( t ) becomes
log x ( t ) = log x 0 + μ ¯ t ,
from where for the expectation of the periodic log return with period Δ t , one finds
1 Δ t log x ( t + Δ t ) / x ( t ) Δ t 0 μ ¯ = d d t log x ( t ) .
The second log-moment is given by
log 2 x ( t ) = log 2 x 0 + 2 μ ¯ log x 0 + σ 2 t + μ ¯ 2 t 2 ,
which for the log-variance yields
log 2 x ( t ) log x ( t ) 2 = σ 2 t .
From Equations (A7), (14), and (16), we find the third and fourth log-moments,
log 3 x ( t ) = log 3 x 0 + 3 μ ¯ log 2 x 0 + σ 2 log x 0 t + 3 μ ¯ μ ¯ log x 0 + σ 2 t 2 + μ ¯ 3 t 3 ,
log 4 x ( t ) = log 4 x 0 + 2 2 μ ¯ log 3 x 0 + 3 σ 2 log 2 x 0 t + 3 2 μ ¯ μ ¯ log 2 x 0 + 2 σ 2 log x 0 + σ 4 t 2 + 4 μ ¯ 2 μ ¯ log x 0 + 3 σ 2 / 2 t 3 + μ ¯ 4 t 4 ,
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 ).

2.2. 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.,
C ( x , t ) = e r ( T t ) K ( x ( T ) K ) f ( x , T , | x 0 , 0 ) d x ,
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],
t + σ 2 x 2 2 2 x 2 r + r x x C BS ( x , t ) = 0 ,
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
t C BS ( x , t ) = σ 2 x 2 2 2 x 2 r + r x x C BS ( x , t ) .
with initial condition C BS ( x , t = 0 ) = max { x K , 0 } , x 0 , and boundary conditions C BS ( x = 0 , t ) = 0 , t 0 , and C ( x , t ) x .
The solution is
C BS ( x 0 , T , K , t ) = N ( d 1 ) x ( t ) N ( d 2 ) K e ( μ σ 2 2 ) ( T t )
d 1 = 1 σ T t log x ( t ) K + μ ( T t )
d 2 = d 1 σ T t ,
where N ( x ) = 1 2 π x e u 2 / 2 d u 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].

2.3. 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])
t f α ( x , t ) = RL D t 1 α μ x x f α ( x , t ) + σ 2 2 2 x 2 x 2 f α ( x , t ) ,
where
RL D t ν f ( t ) = 1 Γ ( 1 ν ) d d t 0 t ( t t ) ν f ( t ) d t
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 L RL D t ν f ( t ) ( s ) = s ν L f ( t ) ( s ) RL I t 1 ν f ( 0 + ) , where RL I t ν f ( t ) = 1 Γ ( ν ) 0 t ( t t ) ν 1 f ( t ) d t is the Riemann–Liouville fractional integral. In order to avoid the somewhat unusual fractional initial condition, alternatively, we could alternatively use the integral version of the equation [46]
f α ( x , t ) f α ( x , 0 ) = RL D t α μ x x f α ( x , t ) + σ 2 2 2 x 2 x 2 f α ( x , t ) ,
where L RL D t α f ( t ) ( s ) = s α L { f ( t ) } ( s ) . 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],
d d u x ( u ) = μ x ( u ) + σ x ( u ) ξ ( u ) ,
d d u T ( u ) = ζ ( u ) .
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 ) = 0 u τ ( u ) d u , 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 s T ( 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 Θ t T ( u ) ,
where Θ ( z ) is the Heaviside theta function. The Laplace transform then yields
h ^ ( u , s ) = u 1 s 0 δ t T ( u ) e s t d t = u 1 s e s T ( u ) = u 1 s e s α u = s α 1 e s α u .
Hence, f α ( x , t ) = δ ( x X ( t ) ) = δ ( x X ( S ( t ) ) = 0 f ( x , u ) h ( u , t ) d t , from where one can easily arrive to the fractional Fokker–Planck Equation (26).
The mean value for sGBM is given by [29,48]
x ( t ) = x 0 E α μ t α ,
where E α ( z ) is the one parameter Mittag–Leffler (ML) function [35,45]
E α ( z ) = k = 0 z k Γ ( α k + 1 ) ,
with ( z C ; R ( α ) > 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 α ( a t α ) ( s ) = s α 1 s α a . The asymptotic behaviour of the mean is given by
x ( t ) x 0 1 + μ t α / Γ ( 1 + α ) e μ t α / Γ ( 1 + α ) , t 1 , α 1 e μ 1 / α t , t 1 .
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]
x 2 ( t ) = x 0 2 E α ( σ 2 + 2 μ ) t α x 2 ( 0 ) 1 + ( σ 2 + 2 μ ) t α / Γ ( 1 + α ) e ( σ 2 + 2 μ ) t α / Γ ( 1 + α ) , t 1 , α 1 e ( σ 2 + 2 μ ) 1 / α t , t 1 .
From here, one concludes that sGBM is an exponentially fast process. Moreover, the third and fourth moments, respectively, become
x 3 ( t ) = x 0 3 E α ( 3 σ 2 + 3 μ ) t α ,
x 4 ( t ) = x 0 4 E α ( 6 σ 2 + 4 μ ) t α .
The first log-moment has the form [29]
log x ( t ) = log x 0 + μ ¯ 0 t t α 1 Γ ( α ) d t = log x 0 + μ ¯ t α Γ ( 1 + α ) ,
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 log x ( t + Δ t ) / x ( t ) Δ t 0 μ ¯ t α 1 Γ ( α ) .
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]
log 2 x ( t ) = log 2 x 0 + 2 μ ¯ log x 0 + σ 2 t α Γ ( 1 + α ) + 2 μ ¯ 2 t 2 α Γ ( 1 + 2 α ) .
from where, for the log-variance, one finds [29]
log 2 x ( t ) log x ( t ) 2 = σ 2 t α Γ ( 1 + α ) + μ ¯ 2 2 Γ ( 1 + 2 α ) 1 Γ 2 ( 1 + α ) t 2 α ,
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
log 3 x ( t ) = log 3 x 0 + 3 μ ¯ log 2 x 0 + σ 2 log x 0 t α Γ ( α + 1 ) + 6 μ ¯ μ ¯ log x 0 + σ 2 t 2 α Γ ( 2 α + 1 ) + 6 μ ¯ 3 t 3 α Γ ( 3 α + 1 ) ,
log 4 x ( t ) = log 4 x 0 + 2 2 μ ¯ log 3 x 0 + 3 σ 2 log 2 x 0 t α Γ ( α + 1 ) + 6 2 μ ¯ μ ¯ log 2 x 0 + 2 σ 2 log x 0 + σ 4 t 2 α Γ ( 2 α + 1 ) + 24 μ ¯ 2 μ ¯ log x 0 + 3 σ 2 / 2 t 3 α Γ ( 3 α + 1 ) + 24 μ ¯ 4 t 4 α Γ ( 4 α + 1 ) .

3. 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 s T ( u ) = e u Ψ ^ ( s ) , with Ψ ^ ( s ) = 1 / η ^ ( s ) .

3.1. 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
X ( t ) = x S ( t ) ,
where X ( t ) is the generalised GBM (gGBM), S ( t ) = inf u > 0 : T ( u ) > t is the operational time, and T ( u ) is an infinite divisible process, i.e., a strictly increasing Lévy motion with [37]
e s T ( 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]
P ( x , t ) = 0 f ( x , u ) h ( u , t ) d u ,
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) reads
P ^ ( x , s ) = L P ( x , t ) = 0 e s t P ( x , t ) d t = 0 f ( x , u ) h ^ ( u , s ) d u ,
where h ^ ( u , s ) = L h ( u , t ) . By considering
h ^ ( u , s ) = Ψ ^ ( s ) s e u Ψ ^ ( s ) = 1 s η ^ ( s ) e u η ^ ( s ) ,
we then have
P ^ ( x , s ) = 1 s η ^ ( s ) 0 f ( x , u ) e u η ^ ( s ) d u = 1 s η ^ ( s ) f ^ x , 1 η ^ ( s ) .
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
s P ^ ( x , s ) P ( x , 0 ) = s η ^ ( s ) μ x x P ^ ( x , s ) + σ 2 2 2 x 2 x 2 P ^ ( x , s ) .
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)
t P ( x , t ) = t 0 t η ( t t ) μ x x P ( x , t ) + σ 2 2 2 x 2 x 2 P ( x , t ) d t ,
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],
P ^ ( x , s ) = 0 1 x 2 π σ 2 u × exp log x log x 0 μ ¯ u 2 2 σ 2 u 1 s η ^ ( s ) e u η ^ ( s ) d u = 1 / [ s η ^ ( s ) ] x μ ¯ 2 + 2 σ 2 / η ^ ( s ) exp log x log x 0 σ 2 μ ¯ 2 + 2 σ 2 / η ^ ( s ) μ ¯ , x > x 0 , 1 , x = x 0 , exp log x log x 0 σ 2 μ ¯ 2 + 2 σ 2 / η ^ ( s ) + μ ¯ , x < x 0 ,
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
0 t γ ( t t ) t P ( x , t ) d t = μ x x P ( x , t ) + σ 2 2 2 x 2 x 2 P ( x , t ) ,
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 γ ( t ) = L 1 s 1 η ^ 1 ( s ) = L 1 1 = δ ( t ) . For sGBM ( η ( t ) = t α 1 / Γ ( α ) , i.e., η ^ ( s ) = s α ) the memory kernel becomes γ ( t ) = L 1 s α 1 = t α / Γ ( 1 α ) , and, thus, Equation (53) reads
C D t α P ( x , t ) = μ x x P ( x , t ) + σ 2 2 2 x 2 x 2 P ( x , t ) ,
where
C D t ν f ( t ) = 1 Γ ( 1 ν ) 0 t ( t t ) ν d d t f ( t ) d t
is the Caputo fractional derivative of order 0 < ν < 1 [45]. 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 ν 1 f ( 0 + ) . We note that, with the appropriate restrictions for η ( t ) and γ ( t ) , both formulations are equivalent.
Remark 3.
For η ( t ) = t α 1 Γ ( α ) , 0 < α < 1 , gGBM corresponds to sGBM. Using the subordination approach one finds [28]
h ^ ( u , s ) = s α 1 e u s α = s α 1 H 0 , 1 1 , 0 u s α ( 0 , 1 ) ,
where H p , q m , n ( z ) is the Fox H-function, see Appendix D. By inverse Laplace transform (A42), one obtains [28]
h ( u , t ) = L 1 h ^ ( u , s ) = t α H 1 , 1 1 , 0 u t α ( 1 α , α ) ( 0 , 1 ) = 1 u H 1 , 1 1 , 0 u t α ( 1 , α ) ( 1 , 1 ) ,
where we applied property (A43). The solution in Laplace space then becomes
P ^ ( x , s ) = 0 1 x 2 π σ 2 u × exp log x log x 0 μ ¯ u 2 2 σ 2 u s α 1 e u s α d u = s α 1 x μ ¯ 2 + 2 σ 2 s α × exp log x log x 0 σ 2 μ ¯ 2 + 2 σ 2 s α μ ¯ , x > x 0 , 1 , x = x 0 , exp log x log x 0 σ 2 μ ¯ 2 + 2 σ 2 s α + μ ¯ , x < x 0 ,
which is obtained in Ref. [48] in a similar way. From here, we can plot the PDF using numerical inverse Laplace transform techniques.

3.2. 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]
C gBS ( x , t ) = e r ( S ( T ) t ) ( x ( S ( T ) ) K ) x = 0 C BS ( x , u ) h ( u , T ) d u ,
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 finds
C ^ gBS ( x , s ) = 1 s η ^ ( s ) C ^ BS ( x , 1 / η ^ ( s ) ) .
Therefore, from Equation (22), the corresponding equation for the option price becomes [48]
t C gBS ( x , t ) = t 0 t η ( t t ) σ 2 x 2 2 2 x 2 r + r x x C gBS ( x , t ) d t .

3.3. Calculation of Moments

The nth moment X n ( t ) = 0 x n P ( x , t ) d x 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
X ^ n ( s ) = x 0 n s 1 1 η ^ ( s ) σ 2 2 n ( n 1 ) + μ n .
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
X ^ ( s ) = x 0 s 1 1 μ η ^ ( s ) ,
X ^ 2 ( s ) = x 0 2 s 1 1 ( σ 2 + 2 μ ) η ^ ( s ) ,
X ^ 3 ( s ) = x 0 3 s 1 1 3 η ^ ( s ) σ 2 + μ ,
X ^ 4 ( s ) = x 0 4 s 1 1 4 η ^ ( s ) 3 σ 2 2 + μ .
The log-moments log n x ( t ) = 0 log n x P ( x , t ) d x , 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
log x ( t ) = log x 0 + μ ¯ 0 t η ( t ) d t .
From here, we find, for the expectation of the periodic log return with period Δ t
1 Δ t log x ( t + Δ t ) / x ( t ) = μ ¯ 1 Δ t t t + Δ t η ( t ) d t = μ ¯ I ( t + Δ t ) I ( t ) Δ t Δ t 0 μ ¯ η ( t ) ,
where I ( t ) = η ( t ) d t , i.e., I ( t ) = η ( t ) . Therefore, the expectation of the periodic log returns behaves as the rate of the first log-moment,
1 Δ t log x ( t + Δ t ) / x ( t ) Δ t 0 d d t log x ( t ) ,
which is proportional to the memory kernel η ( t ) . Moreover, for the second log-moment, we find
log 2 x ( t ) = log 2 x 0 + 0 t η ( t t ) 2 μ ¯ log x 0 + μ ¯ 0 t η ( t ) d t + σ 2 d t ,
from where the log-variance becomes
log 2 x ( t ) log x ( t ) 2 = σ 2 0 t η ( t ) d t + μ ¯ 2 2 0 t η ( t t ) 0 t η ( t ) d t d t 0 t η ( t ) d t 2 .
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 ).

3.4. Exponentially Truncated Subdiffusive GBM

As an example for another memory kernel in gGBM, we consider a power-law memory kernel with exponential truncation,
η ( t ) = t α 1 Γ ( α ) e t τ ,
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,
η ^ ( s ) = ( s + τ 1 ) α ,
where we use the shift rule of the Laplace transform, L e a t f ( t ) = F ^ ( s + a ) , for F ^ ( s ) = L f ( t ) .
The mean value reads,
x ( t ) = x 0 L 1 s 1 1 μ ( s + τ 1 ) α ( t ) = x 0 L 1 s + τ 1 s ( s + τ 1 ) α 1 ( s + τ 1 ) α μ ( t ) = x 0 e t / τ E α μ t α + τ 1 0 t e t / τ E α μ t α d t ,
and the MSD is
x 2 ( t ) = x 0 2 L 1 s 1 1 ( σ 2 + 2 μ ) ( s + τ 1 ) α ( t ) = x 0 2 L 1 s + τ 1 s ( s + τ 1 ) α 1 ( s + τ 1 ) α ( σ 2 + 2 μ ) ( t ) = x 0 2 e t / τ E α ( σ 2 + 2 μ ) t α + τ 1 0 t e t / τ E α ( σ 2 + 2 μ ) t α d t .
The third and fourth moments become
x 3 ( t ) = x 0 3 L 1 s 1 1 3 ( s + τ 1 ) α σ 2 + μ ( t ) = x 0 3 L 1 s + τ 1 s ( s + τ 1 ) α 1 ( s + τ 1 ) α 3 σ 2 + μ ( t ) = x 0 3 e t / τ E α ( 3 σ 2 + 3 μ ) t α + τ 1 0 t e t / τ E α ( 3 σ 2 + 3 μ ) t α d t ,
x 4 ( t ) = x 0 4 L 1 s 1 1 ( s + τ 1 ) α 6 σ 2 + 4 μ ( t ) = x 0 4 L 1 s + τ 1 s ( s + τ 1 ) α 1 ( s + τ 1 ) α 6 σ 2 + 4 μ ( t ) = x 0 4 e t / τ E α ( 6 σ 2 + 4 μ ) t α + τ 1 0 t e t / τ E α ( 6 σ 2 + 4 μ ) t α d t ,
respectively.
From the general result for the log-mean, we find that
log x ( t ) = log x 0 + μ ¯ τ α γ ( α , t / τ ) Γ ( α ) = log x 0 + μ ¯ e t / τ t α E 1 , α + 1 ( t / τ ) ,
where ( z , β C ; R ( α ) > 0 ). Here, γ ( a , z ) = 0 z t a 1 e t d t = Γ ( a ) e z z a E 1 , a + 1 ( z ) is the incomplete gamma function, and
E α , β ( z ) = k = 0 z k Γ ( α k + β ) .
is the two parameter ML function [45]. Note that, the Laplace transform of the two parameter ML function reads L t β 1 E α , β ( a t α ) ( 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
1 Δ t log x ( t + Δ t ) / x ( t ) Δ t 0 μ ¯ t α 1 Γ ( α ) e t / τ = d d t log x ( t ) .
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
log 2 x ( t ) = log 2 x 0 + 2 μ ¯ log x 0 + σ 2 τ α γ ( α , t / τ ) Γ ( α ) + 2 μ ¯ 2 τ 2 α γ ( 2 α , t / τ ) Γ ( 2 α ) ,
from where the log-variance becomes
log 2 x ( t ) log x ( t ) 2 = σ 2 τ α γ ( α , t / τ ) Γ ( α ) + μ ¯ 2 τ 2 α 2 γ ( 2 α , t / τ ) Γ ( 2 α ) γ 2 ( α , t / τ ) Γ 2 ( α ) .
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 3 x ( t ) = log 3 x 0 + 3 μ ¯ log 2 x 0 + σ 2 log x 0 τ α γ ( α , t / τ ) Γ ( α ) + 6 μ ¯ μ ¯ log x 0 + σ 2 τ 2 α γ ( 2 α , t / τ ) Γ ( 2 α ) + 6 μ ¯ 3 τ 3 α γ ( 3 α , t / τ ) Γ ( 3 α ) ,
log 4 x ( t ) = log 4 x 0 + 2 2 μ ¯ log 3 x 0 + 3 σ 2 log 2 x 0 τ α γ ( α , t / τ ) Γ ( α ) + 6 2 μ ¯ μ ¯ log 2 x 0 + 2 σ 2 log x 0 + σ 4 τ 2 α γ ( 2 α , t / τ ) Γ ( 2 α ) + 24 μ ¯ 2 μ ¯ log x 0 + 3 σ 2 / 2 τ 3 α γ ( 3 α , t / τ ) Γ ( 3 α ) + 24 μ ¯ 4 τ 4 α γ ( 4 α , t / τ ) Γ ( 4 α ) .
The subordination function, in this case, is given by
h ^ ( u , s ) = ( s + τ 1 ) α s e u ( s + τ 1 ) α = 1 + ( s τ ) 1 ( s + τ 1 ) α 1 e u ( s + τ 1 ) α , h ( u , t ) = e t / τ H 1 , 1 1 , 0 u t α ( 1 , α ) ( 1 , 1 ) + τ 1 0 t e t / τ H 1 , 1 1 , 0 u t α ( 1 , α ) ( 1 , 1 ) d t ,
from where one can analyse the PDF P ( x , t ) .

3.5. Combined Standard and Subdiffusive GBM

As another application, let us consider the combination of GBM and sGBM, represented by the memory kernel
η ( t ) = w 1 t α 1 Γ ( α ) + w 2 ,
where 0 < α < 1 , w 1 + w 2 = 1 , and
η ^ ( s ) = w 1 s α + w 2 s 1 .
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
x ( t ) = x 0 L 1 s 1 1 μ w 1 s α + w 2 s 1 ( t ) = x 0 n = 0 w 1 n μ n t α n E 1 , α n + 1 n + 1 w 2 μ t = x 0 n = 0 w 1 n μ n t α n Γ ( α n + 1 ) 1 F 1 n + 1 ; α n + 1 ; w 2 μ t ,
where 1 F 1 a ; b ; z = k = 0 ( a ) k ( b ) k z k k ! is the Kummer confluent hypergeometric function,
E α , β γ ( z ) = n = 0 ( γ ) n Γ ( α n + β ) z n n ! ,
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 α , β γ ( a t α ) ( 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
x ( t ) = x 0 n = 0 μ n t α n Γ ( α n + 1 ) = E α μ t α ,
as it should be for the sGBM. For the second, third, and fourth moments, we find
x 2 ( t ) = x 0 2 n = 0 w 1 n ( σ 2 + 2 μ ) n t α n E 1 , α n + 1 n + 1 w 2 ( σ 2 + 2 μ ) t ,
x 3 ( t ) = x 0 3 n = 0 w 1 n ( 3 σ 2 + 3 μ ) n t α n E 1 , α n + 1 n + 1 w 2 ( 3 σ 2 + 3 μ ) t ,
x 4 ( t ) = x 0 4 n = 0 w 1 n ( 6 σ 2 + 4 μ ) n t α n E 1 , α n + 1 n + 1 w 2 ( 6 σ 2 + 4 μ ) t .
Following the same procedure as previously, for the log-mean we find
log x ( t ) = log x 0 + μ ¯ w 1 t α Γ ( α + 1 ) + w 2 t .
and for the expectation of the periodic log return with period Δ t ,
1 Δ t log x ( t + Δ t ) / x ( t ) Δ t 0 μ ¯ w 1 t α 1 Γ ( α ) + w 2 = d d t log x ( 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
log 2 x ( t ) = log 2 x 0 + 2 μ ¯ log x 0 + σ 2 w 1 t α Γ ( α + 1 ) + w 2 t + 2 μ ¯ 2 w 1 2 t 2 α Γ ( 2 α + 1 ) + 2 w 1 w 2 t α + 1 Γ ( α + 2 ) + w 2 2 t 2 2 ,
from where the log-variance becomes
log 2 x ( t ) log x ( t ) 2 = σ 2 w 1 t α Γ ( α + 1 ) + w 2 t + μ ¯ 2 w 1 2 t 2 α 2 Γ ( 2 α + 1 ) 1 Γ 2 ( α + 1 ) + 2 μ ¯ 2 w 1 w 2 t α + 1 2 Γ ( α + 2 ) 1 Γ ( α + 1 ) .
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 3 x ( t ) = log 3 x 0 + 3 μ ¯ log 2 x 0 + σ 2 log x 0 w 1 t α Γ ( α + 1 ) + w 2 t + 6 μ ¯ μ ¯ log x 0 + σ 2 w 1 2 t 2 α Γ ( 2 α + 1 ) + 2 w 1 w 2 t α + 1 Γ ( α + 2 ) + w 2 2 t 2 2 + 6 μ ¯ 3 w 1 3 t 3 α Γ ( 3 α + 1 ) + 3 w 1 2 w 2 t 2 α + 1 Γ ( 2 α + 2 ) + 3 w 1 w 2 2 t α + 2 Γ ( α + 3 ) + w 2 3 t 3 6 ,
log 4 x ( t ) = log 4 x 0 + 2 2 μ ¯ log 3 x 0 + 3 σ 2 log 2 x 0 w 1 t α Γ ( α + 1 ) + w 2 t + 6 2 μ ¯ μ ¯ log 2 x 0 + 2 σ 2 log x 0 + σ 4 w 1 2 t 2 α Γ ( 2 α + 1 ) + 2 w 1 w 2 t α + 1 Γ ( α + 2 ) + w 2 2 t 2 2 + 24 μ ¯ 2 μ ¯ log x 0 + 3 σ 2 / 2 w 1 3 t 3 α Γ ( 3 α + 1 ) + 3 w 1 2 w 2 t 2 α + 1 Γ ( 2 α + 2 ) + 3 w 1 w 2 2 t α + 2 Γ ( α + 3 ) + w 2 3 t 3 6 + 24 μ ¯ 4 w 1 4 t 4 α Γ ( 4 α + 1 ) + 4 w 1 3 w 2 t 3 α + 1 Γ ( 3 α + 2 ) + 6 w 1 2 w 2 2 t 2 α + 2 Γ ( 2 α + 2 ) + 4 w 1 w 2 3 t α + 3 Γ ( α + 4 ) + w 2 4 t 4 24 .
The subordination function for this case is given by
h ^ ( u , s ) = 1 w 1 + w 2 s 1 α e u w 1 s 1 + w 2 s α ,
where the Lévy exponent is Ψ ^ ( s ) = w 1 s 1 + w 2 s α 1 .

3.6. Mix of Subdiffusive GBMs

We may further analyse the case of a mix of two sGBM with different power-law memory functions,
η ( t ) = w 1 t α 1 1 Γ ( α 1 ) + w 2 t α 2 1 Γ ( α 2 ) ,
where 0 < α 1 < α 2 < 1 , w 1 + w 2 = 1 , and
η ^ ( s ) = w 1 s α 1 + w 2 s α 2 .
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
x ( t ) = x 0 L 1 s 1 1 μ w 1 s α 1 + w 2 s α 2 ( t ) = x 0 n = 0 w 1 n μ n t α 1 n E α 2 , α 1 n + 1 n + 1 w 2 μ t α 2 ,
while, for the second, third, and fourth moments, respectively, we obtain
x 2 ( t ) = x 0 2 n = 0 w 1 n ( σ 2 + 2 μ ) n t α 1 n E α 2 , α 1 n + 1 n + 1 w 2 ( σ 2 + 2 μ ) t α 2 ,
x 3 ( t ) = x 0 3 n = 0 w 1 n ( 3 σ 2 + 3 μ ) n t α 1 n E α 2 , α 1 n + 1 n + 1 w 2 ( 3 σ 2 + 3 μ ) t α 2 ,
x 4 ( t ) = x 0 4 n = 0 w 1 n ( 6 σ 2 + 4 μ ) n t α 1 n E α 2 , α 1 n + 1 n + 1 w 2 ( 6 σ 2 + 4 μ ) t α 2 .
Similarly, the log-mean yields
log x ( t ) = log x 0 + μ ¯ w 1 t α 1 Γ ( α 1 + 1 ) + w 2 t α 2 Γ ( α 2 + 1 ) .
The expectation of the log return with period Δ t , then becomes
1 Δ t log x ( t + Δ t ) / x ( t ) Δ t 0 μ ¯ w 1 t α 1 1 Γ ( α 1 ) + w 2 t α 2 1 Γ ( α 2 ) = d d t log x ( t ) .
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
log 2 x ( t ) = log 2 x 0 + 2 μ ¯ log x 0 + σ 2 w 1 t α 1 Γ ( α 1 + 1 ) + w 2 t α 2 Γ ( α 2 + 1 ) + 2 μ ¯ 2 w 1 2 t 2 α 1 Γ ( 2 α 1 + 1 ) + 2 w 1 w 2 t α 1 + α 2 Γ ( α 1 + α 2 + 1 ) + w 2 2 t 2 α 2 Γ ( 2 α 2 + 1 ) ,
and for the log-variance we find
log 2 x ( t ) log x ( t ) 2 = σ 2 w 1 t α 1 Γ ( α 1 + 1 ) + w 2 t α 2 Γ ( α 2 + 1 ) + 2 μ ¯ 2 w 1 w 2 t α 1 + α 2 2 Γ ( α 1 + α 2 + 1 ) 1 Γ ( α 1 + 1 ) Γ ( α 2 + 1 ) + μ ¯ 2 w 1 2 t 2 α 1 2 Γ ( 2 α 1 + 1 ) 1 Γ 2 ( α 1 + 1 ) + μ ¯ 2 w 2 2 t 2 α 2 2 Γ ( 2 α 2 + 1 ) 1 Γ 2 ( α 2 + 1 ) .
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
log 3 x ( t ) = log 3 x 0 + 3 μ ¯ log 2 x 0 + σ 2 log x 0 w 1 t α 1 Γ ( α 1 + 1 ) + w 2 t α 2 Γ ( α 2 + 1 ) + 6 μ ¯ μ ¯ log x 0 + σ 2 w 1 2 t 2 α 1 Γ ( 2 α 1 + 1 ) + 2 w 1 w 2 t α 1 + α 2 Γ ( α 1 + α 2 + 1 ) + w 2 2 t 2 α 2 Γ ( 2 α 2 + 1 ) + 6 μ ¯ 3 w 1 3 t 3 α 1 Γ ( 3 α 1 + 1 ) + 3 w 1 2 w 2 t 2 α 1 + α 2 Γ ( 2 α 1 + α 2 + 1 ) + 3 w 1 w 2 2 t α 1 + 2 α 2 Γ ( α 1 + 2 α 2 + 1 ) + w 2 3 t 3 α 2 Γ ( 3 α 2 + 1 ) ,
log 4 x ( t ) = log 4 x 0 + 2 2 μ ¯ log 3 x 0 + 3 σ 2 log 2 x 0 w 1 t α 1 Γ ( α 1 + 1 ) + w 2 t α 2 Γ ( α 2 + 1 ) + 6 2 μ ¯ μ ¯ log 2 x 0 + 2 σ 2 log x 0 + σ 4 w 1 2 t 2 α 1 Γ ( 2 α 1 + 1 ) + 2 w 1 w 2 t α 1 + α 2 Γ ( α 1 + α 2 + 1 ) + w 2 2 t 2 α 2 Γ ( 2 α 2 + 1 ) + 24 μ ¯ 2 μ ¯ log x 0 + 3 σ 2 / 2 w 1 3 t 3 α 1 Γ ( 3 α 1 + 1 ) + 3 w 1 2 w 2 t 2 α 1 + α 2 Γ ( 2 α 1 + α 2 + 1 ) + 3 w 1 w 2 2 t α 1 + 2 α 2 Γ ( α 1 + 2 α 2 + 1 ) + w 2 3 t 3 α 2 Γ ( 3 α 2 + 1 ) + 24 μ ¯ 4 w 1 4 t 4 α 1 Γ ( 4 α 1 + 1 ) + 4 w 1 3 w 2 t 3 α 1 + α 2 Γ ( 3 α 1 + α 2 + 1 ) + 6 w 1 2 w 2 2 t 2 α 1 + 2 α 2 Γ ( 2 α 1 + 2 α 2 + 1 ) + 4 w 1 w 2 3 t α 1 + 3 α 2 Γ ( α 1 + 3 α 2 + 1 ) + w 2 4 t 4 α 2 Γ ( 4 α 2 + 1 ) .
For the mix of subdiffusive GBMs, the subordination function becomes
h ^ ( u , s ) = 1 w 1 s 1 α 1 + w 2 s 1 α 2 e u w 1 s α 1 + w 2 s α 2 ,
where the Lévy exponent is Ψ ^ ( s ) = w 1 s α 1 + w 2 s α 2 1 , see Refs. [64,65].

4. Numerical Experiments

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”.
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.

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 ).
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.

6. 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].

Author Contributions

All authors conceived this research. V.S., T.S., and L.B. performed the analysis and prepared the figures. All authors wrote and reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge funding from the German Science Foundation (DFG grant ME 1535/6-1). T.S. was supported by the Alexander von Humboldt Foundation. R.M. acknowledges an Alexander von Humboldt Polish Honorary Research Scholarship from the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP).

Acknowledgments

T.S. acknowledges Andrey Cherstvy for fruitful discussions and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
GBMGeometric Brownian motion
sGBMSubdiffusive geometric Brownian motion
gGBMGeneralised geometric Brownian motion
BSBlack-Scholes
CTRWContinuous time random walk
MSDMean squared displacement
MLMittag-Leffler
ATMat-the-money
RMSERoot mean squared error
TSLATesla
AAPLApple

Appendix A. Solution of the Fokker-Planck Equation for Standard GBM

The solution of Equation (3) can be found by using the Laplace-Mellin transform method [74]. The Laplace transform is defined by
F ^ ( s ) = L f ( t ) ( s ) = 0 f ( t ) e s t d t ,
while the Mellin transform as [75]
F ˜ ( q ) = M f ( x ) ( q ) = 0 x q 1 f ( x ) d x .
The inverse Mellin transform then reads
f ( x ) = M 1 F ˜ ( q ) ( x ) = 1 2 π ı c ı c + ı x q F ˜ ( q ) d q .
Therefore, by performing Laplace transform in respect to t and Mellin transform in respect to x in Equation (3), we have
F ^ ˜ ( q , s ) = x 0 q 1 × 1 s σ 2 2 ( q 1 ) ( q 2 ) + μ ( q 1 ) ,
where we use M δ ( x x 0 ) ( q ) = x 0 q 1 . Then the inverse Laplace transform yields
F ˜ ( q , t ) = x 0 q 1 × exp σ 2 2 q + 1 2 2 μ σ 2 3 2 t μ σ 2 2 2 2 σ 2 t ,
where we use L 1 1 s a ( t ) = e a t . Applying the inverse Mellin transform and looking for the solution in the form of the convolution integral of two functions [75], M h ( x ) ( q ) = H ˜ ( q ) and M g ( x ) ( q ) = G ˜ ( q ) ,
M 1 H ˜ ( q ) G ˜ ( q ) ( x ) = 0 h ( r ) g ( x / r ) d r r ,
we obtain the solution of the Fokker-Planck equation for GBM
f ( x , t ) = 0 δ ( r x 0 ) × exp log x r μ σ 2 2 t 2 2 σ 2 t ( x / r ) 2 π σ 2 t d r r = 1 x 2 π σ 2 t × exp log x x 0 μ σ 2 2 t 2 2 σ 2 t .
Here we used that
h ( x ) = M 1 x 0 q 1 ( x ) = δ ( x x 0 )
and
g ( x ) = M 1 exp σ 2 2 q + 1 2 2 μ σ 2 3 2 t μ σ 2 2 2 2 σ 2 t ( x ) = 1 x 2 π σ 2 t × exp log x μ σ 2 2 t 2 2 σ 2 t .
We used the properties of the inverse Mellin transform [75], M 1 f ( q + a ) ( x ) = x a M 1 f ( q ) and
M 1 exp α q 2 ( x ) = 1 4 π α e x 2 4 α .
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 ) d x 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
x 2 ( t ) x ( t ) 2 = x 2 ( 0 ) e 2 μ t e σ 2 t 1 .
The log-moments log n x = 0 log n x P ( x , t ) d x 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)
t log n x ( t ) = μ σ 2 2 n log n 1 x ( t ) + σ 2 2 n ( n 1 ) log n 2 x ( t ) .
From here it follows that t log 0 x ( t ) = 0 , i.e., log 0 x ( t ) = log 0 x ( 0 ) = 1 . The case n = 1 yields the mean value of the logarithm of x ( t ) ,
t log x ( t ) = μ σ 2 2 log 0 x ( t ) = 1
i.e.,
log x ( t ) = log x ( 0 ) + μ σ 2 2 t .
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 )
which is given by
log 2 x ( t ) = log 2 x ( 0 ) + μ σ 2 2 2 t 2 + 2 μ σ 2 2 log x ( 0 ) t + σ 2 t .
Therefore, for the log-variance one finds linear dependence on time
log 2 x ( t ) log x ( t ) 2 = σ 2 t .

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],
λ ( x i , t , x j ) = p r ( x j , t ) δ ( x i [ 1 + u ] x j ) + p l ( x j , t ) δ ( x i x j / [ 1 + u ] ) ,
and a waiting time PDF ψ ( t ) , related to the survival probability by
ϕ ( t ) = 1 0 t ψ ( t ) d t , i . e . , ϕ ^ ( s ) = 1 ψ ^ ( s ) s .
By substitution in the master equation [58]
t ρ ( x i , t ) = j λ ( x i , t , x j ) 0 t K ( t t ) ρ ( x j , t ) d t 0 t K ( t t ) ρ ( x i , t ) d t ,
where K ( t ) = L 1 ψ ^ ( s ) / ϕ ^ ( s ) , one finds
t ρ ( x i , t ) = p r x i 1 + u , t 0 t K ( t t ) ρ x i 1 + u , t d t + p l [ 1 + u ] x i , t 0 t K ( t t ) ρ [ 1 + u ] x i , t d t 0 t K ( t t ) ρ ( x i , t ) d t .
We consider generalised waiting time PDF, which in the Laplace space has the form [61,78]
ψ ^ ( s ) = 1 1 + τ η / η ^ ( s ) ,
where τ η is a time parameter, which depends on η ( t ) . Therefore,
ϕ ^ ( s ) = τ η / η ^ ( s ) s 1 + τ η / η ^ ( s ) ,
and
K ^ ( s ) = 1 τ η s × η ^ ( s ) ,
from where we find that 0 t K ( t t ) f ( t ) d t 1 τ η d d t 0 t η ( t t ) f ( t ) d t . From Equation (A14) then we obtain
t ρ ( x i , t ) = 1 τ η p r x i 1 + u , t t 0 t η ( t t ) ρ x i 1 + u , t d t + 1 τ η p l [ 1 + u ] x i , t t 0 t η ( t t ) ρ [ 1 + u ] x i , t d t 1 τ η t 0 t η ( t t ) ρ ( x i , t ) d t .
Let us now consider the diffusion limit ( u 0 and τ η 0 ) of Equation (A15). From the normalisation condition of the PDF ρ ( x , t ) given by i ρ ( x i , t ) = 1 and by using position-dependent lattice spacing Δ x i = u x i , one finds i ρ u ( x i , t ) u x i Δ x i = 1 , such that lim Δ x i 0 i ρ u ( x i , t ) u x i Δ x i = 1 [58]. By defining the function P u ( x , t ) = ρ u ( x , t ) / [ u x ] , one concludes that P ( x , t ) = lim u 0 P u ( x , t ) is normalised, i.e., 0 P ( x , t ) d x = 1 . By introducing B u ( x i , t ) = p r ( x i , t ) p l ( x i , t ) and b 0 ( x , t ) = lim u 0 u B u ( x , t ) , in the diffusion limit u 0 and τ η 0 , where we assume that B 0 ( x , t ) = lim u 0 B u ( x , t ) = 0 [58], we arrive to the following Fokker-Planck equation
t P ( x , t ) = D t 0 t η ( t t ) x x 2 x 1 k B T x 2 F ( x ) P ( x , t ) d t ,
where D = lim u , τ η 0 u 2 / [ 2 τ η ] , F ( x ) = V ( x ) = k B T [ 2 b 0 ( x ) 1 ] / x , and P ( x , t ) exp V ( x ) k B T is obtained from the long time steady state Boltzmann distribution [58]. 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 = σ 2 / 2 and v = 2 μ / D the Fokker-Planck equation becomes
t P ( x , t ) = t 0 t η ( t t ) x σ 2 x 2 2 x + [ σ 2 μ ] x P ( x , t ) d t ,
which can be rewritten in the form of Equation (51).

Appendix C. General Results for nth Moment

If we multiply both sides of Equation (51) by x n , and integrate over x we find the nth moment x n ( t ) = 0 x n P ( x , t ) d x ,
t x n ( t ) = σ 2 2 n ( n 1 ) + μ n μ d d t 0 t η ( t t ) x n ( t ) d t ,
from where in the Laplace space it reads
x ^ n ( s ) = s 1 1 η ^ ( s ) σ 2 2 n ( n 1 ) + μ n x n ( 0 ) .
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 x ( t ) = μ d d t 0 t η ( t t ) x ( t ) d t ,
and its Laplace pair
x ^ ( s ) = s 1 1 μ η ^ ( s ) x ( 0 ) .
In terms of the memory kernel γ ( t ) , Equation (A21) reads
x ^ ( s ) = γ ^ ( s ) s γ ^ ( s ) μ x ( 0 ) .
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,
t x 2 ( t ) = ( σ 2 + 2 μ ) d d t 0 t η ( t t ) x 2 ( t ) d t ,
and its Laplace pair
x ^ 2 ( s ) = s 1 1 ( σ 2 + 2 μ ) η ^ ( s ) x 2 ( 0 ) ,
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 ) d x , which satisfy the following integral equation
t log n x ( t ) = t 0 t η ( t t ) μ σ 2 2 n log n 1 x ( t ) + σ 2 2 n ( n 1 ) log n 2 x ( t ) d t .
Thus, we find t log 0 x ( t ) = 0 , i.e., log 0 x ( t ) = log 0 x ( 0 ) = 1 . The mean value ( n = 1 ) becomes
t log x ( t ) = μ σ 2 2 t 0 t η ( t t ) log 0 x ( t ) = 1 d t
from where it follows
log x ( t ) = log x ( 0 ) + μ σ 2 2 0 t η ( t ) d t .
For the expectation of the periodic log return with period Δ t , we find
1 Δ t log x ( t + Δ t ) / x ( t ) = μ σ 2 2 1 Δ t t t + Δ t η ( t ) d t = μ σ 2 2 I ( t + Δ t ) I ( t ) Δ t Δ t 0 μ σ 2 2 η ( t ) ,
where I ( t ) = η ( t ) d t , i.e., I ( t ) = η ( t ) . Therefore, the expectation of the periodic log returns behaves as the rate of the first log-moment,
1 Δ t log x ( t + Δ t ) / x ( t ) Δ t 0 d d t log x ( t ) .
For n = 2 we obtain the second log-moment
t log 2 x ( t ) = t 0 t η ( t t ) 2 μ σ 2 2 log x ( t ) + σ 2 log 0 x ( t ) d t
i.e.,
log 2 x ( t ) = log 2 x ( 0 ) + 2 μ σ 2 2 log x ( 0 ) + σ 2 0 t η ( t ) d t + 2 μ σ 2 2 2 0 t η ( t t ) 0 t η ( t ) d t d t ,
and the log-variance becomes
log 2 x ( t ) log x ( t ) 2 = σ 2 0 t η ( t ) d t + μ σ 2 2 2 2 0 t η ( t t ) 0 t η ( t ) d t d t 0 t η ( t ) d t 2 .
Following the same procedure, we find the third ( n = 3 ) and fourth ( n = 4 ) log-moments,
log 3 x ( t ) = log 3 x ( 0 ) + 3 μ σ 2 2 log 2 x ( 0 ) + σ 2 log x ( 0 ) I 1 ( t ) + 6 μ σ 2 2 μ σ 2 2 log x ( 0 ) + σ 2 I 2 ( t ) + 6 μ σ 2 2 3 I 3 ,
log 4 x ( t ) = log 4 x ( 0 ) + 2 2 μ σ 2 2 log 3 x ( 0 ) + 3 σ 2 log 2 x ( 0 ) I 1 ( t ) + 6 2 μ σ 2 2 μ σ 2 2 log 2 x ( 0 ) + 2 σ 2 log x ( 0 ) + σ 4 I 2 ( t ) + 24 μ σ 2 2 2 μ σ 2 2 log x ( 0 ) + 3 σ 2 / 2 I 3 ( t ) + 24 μ σ 2 2 4 I 4 ( t ) ,
where
I 1 ( t ) = 0 t η ( t 1 ) d t 1 ,
I 2 ( t ) = 0 t η ( t t 1 ) 0 t 1 η ( t 2 ) d t 2 d t 1 ,
I 3 ( t ) = 0 t η ( t t 1 ) 0 t 1 η ( t 1 t 2 ) 0 t 2 η ( t 3 ) d t 3 d t 2 d t 1 ,
I 4 ( t ) = 0 t η ( t t 1 ) 0 t 1 η ( t 1 t 2 ) 0 t 2 η ( t 2 t 3 ) 0 t 3 η ( t 4 ) d t 4 d t 3 d t 2 d t 1 .

Appendix D. Fox H-Function

The Fox H-function is defined by [79]
H p , q m , n z ( a 1 , A 1 ) , , ( a p , A p ) ( b 1 , B 1 ) , , ( b q , B q ) = H p , q m , n z ( a p , A p ) ( b q , B q ) = 1 2 π ı Ω θ ( s ) z s d s ,
where θ ( s ) is given by θ ( s ) = j = 1 m Γ ( b j B j s ) j = 1 n Γ ( 1 a j + A j s ) j = m + 1 q Γ ( 1 b j + B j s ) j = n + 1 p Γ ( a j A j s ) , 0 n p , 1 m q , a i , b j C , A i , B j R + , i = 1 , . . . , p , j = 1 , . . . , q . The contour Ω starting at c ı and ending at c + ı separates the poles of the function Γ ( b j + B j s ) , j = 1 , . . . , m from those of the function Γ ( 1 a i A i s ) , i = 1 , . . . , n .
A special case of the Fox H-function is the exponential function [79],
e z = H 0 , 1 1 , 0 z ( 0 , 1 ) .
The inverse Laplace transform of the Fox H-function reads [79]
L 1 s ρ H p , q m , n a s σ ( a p , A p ) ( b q , B q ) ( t ) = t ρ 1 H p + 1 , q m , n a t σ ( a p , A p ) , ( ρ , σ ) ( b q , B q ) .
The Fox H-functions have the following property [79]
z k H p , q m , n z ( a p , A p ) ( b q , B q ) = H p , q m , n z ( a p + k A p , A p ) ( b q + k B q , B q ) .

References

  1. Stojkoski, V.; Utkovski, Z.; Basnarkov, L.; Kocarev, L. Cooperation dynamics in networked geometric Brownian motion. Phys. Rev. E 2019, 99, 062312. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Stojkoski, V.; Karbevski, M.; Utkovski, Z.; Basnarkov, L.; Kocarev, L. Evolution of cooperation in populations with heterogeneous multiplicative resource dynamics. arXiv 2019, arXiv:1912.09205. [Google Scholar]
  3. Peters, O.; Klein, W. Ergodicity breaking in geometric Brownian motion. Phys. Rev. Lett. 2013, 110, 100603. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Aitchison, J.; Brown, J.A. The Lognormal Distribution with Special Reference to Its Uses in Economics; Cambridge University Press: Cambridge, UK, 1957. [Google Scholar]
  5. Redner, S. Random multiplicative processes: An elementary tutorial. Am. J. Phys. 1990, 58, 267–273. [Google Scholar] [CrossRef]
  6. Black, F.; Scholes, M. The pricing of options and corporate liabilities. J. Polit. Econ. 1973, 81, 637–654. [Google Scholar] [CrossRef] [Green Version]
  7. Merton, R.C. Optimum consumption and portfolio rules in a continuous-time model. In Stochastic Optimization Models in Finance; Elsevier: Amsterdam, The Netherlands, 1975; pp. 621–661. [Google Scholar]
  8. Merton, R.C. Option pricing when underlying stock returns are discontinuous. J. Financ. Econ. 1976, 3, 125–144. [Google Scholar] [CrossRef] [Green Version]
  9. Peters, O. Optimal leverage from non-ergodicity. Quant. Financ. 2011, 11, 1593–1602. [Google Scholar] [CrossRef] [Green Version]
  10. Oshanin, G.; Schehr, G. Two stock options at the races: Black–Scholes forecasts. Quant. Financ. 2012, 12, 1325–1333. [Google Scholar] [CrossRef] [Green Version]
  11. Hentschel, H.; Procaccia, I. Fractal nature of turbulence as manifested in turbulent diffusion. Phys. Rev. A 1983, 27, 1266. [Google Scholar] [CrossRef]
  12. Heidernätsch, M.S.M. On the Diffusion in Inhomogeneous Systems. Ph.D. Thesis, Technischen Universität Chemnitz, Chemnitz, Germany, 2015. [Google Scholar]
  13. Baskin, E.; Iomin, A. Superdiffusion on a comb structure. Phys. Rev. Lett. 2004, 93, 120603. [Google Scholar] [CrossRef]
  14. Taleb, N.N. The Black Swan: The Impact of the Highly Improbable; Random House: New York, NY, USA, 2007; Volume 2. [Google Scholar]
  15. Cox, J. Notes on Option Pricing I: Constant Elasticity of Variance Diffusions; Unpublished Note; Stanford University, Graduate School of Business: Stanford, CA, USA, 1975. [Google Scholar]
  16. Heston, S.L. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Rev. Financ. Stud. 1993, 6, 327–343. [Google Scholar] [CrossRef] [Green Version]
  17. Hagan, P.S.; Kumar, D.; Lesniewski, A.S.; Woodward, D.E. Managing smile risk. Best Wilmott 2002, 1, 249–296. [Google Scholar]
  18. Dupire, B. Pricing with a smile. Risk 1994, 7, 18–20. [Google Scholar]
  19. Derman, E.; Kani, I. The volatility smile and its implied tree. Goldman Sachs Quant. Strat. Res. Notes 1994, 2, 45–60. [Google Scholar]
  20. Engle, R.F.; Bollerslev, T. Modelling the persistence of conditional variances. Econ. Rev. 1986, 5, 1–50. [Google Scholar] [CrossRef]
  21. Bollerslev, T. Generalized autoregressive conditional heteroskedasticity. J. Econ. 1986, 31, 307–327. [Google Scholar] [CrossRef] [Green Version]
  22. Matacz, A. Financial modeling and option theory with the truncated Lévy process. Int. J. Theor. Appl. Financ. 2000, 3, 143–160. [Google Scholar] [CrossRef]
  23. Borland, L. A theory of non-Gaussian option pricing. Quant. Financ. 2002, 2, 415–431. [Google Scholar]
  24. Borland, L.; Bouchaud, J.P. A non-Gaussian option pricing model with skew. Quant. Financ. 2004, 4, 499–514. [Google Scholar] [CrossRef]
  25. Moriconi, L. Delta hedged option valuation with underlying non-Gaussian returns. Phys. A Stat. Mech. Appl. 2007, 380, 343–350. [Google Scholar] [CrossRef] [Green Version]
  26. Cassidy, D.T.; Hamp, M.J.; Ouyed, R. Pricing European options with a log Student’s t-distribution: A Gosset formula. Phys. A Stat. Mech. Appl. 2010, 389, 5736–5748. [Google Scholar] [CrossRef] [Green Version]
  27. Basnarkov, L.; Stojkoski, V.; Utkovski, Z.; Kocarev, L. Option Pricing with Heavy-Tailed Distributions of Logarithmic Returns. arXiv 2018, arXiv:1807.01756. [Google Scholar] [CrossRef] [Green Version]
  28. Magdziarz, M. Black-Scholes formula in subdiffusive regime. J. Stat. Phys. 2009, 136, 553–564. [Google Scholar] [CrossRef]
  29. Angstmann, C.; Henry, B.; McGann, A. Time-fractional geometric Brownian motion from continuous time random walks. Physica A 2019, 526, 121002. [Google Scholar] [CrossRef]
  30. Krzyżanowski, G.; Magdziarz, M.; Płociniczak, Ł. A weighted finite difference method for subdiffusive Black–Scholes model. Comput. Math. Appl. 2020, 80, 653–670. [Google Scholar]
  31. McCauley, J.L.; Gunaratne, G.H.; Bassler, K.E. Hurst exponents, Markov processes, and fractional Brownian motion. Phys. A Stat. Mech. Appl. 2007, 379, 1–9. [Google Scholar] [CrossRef] [Green Version]
  32. Vellekoop, M.; Nieuwenhuis, H. On option pricing models in the presence of heavy tails. Quant. Financ. 2007, 7, 563–573. [Google Scholar] [CrossRef] [Green Version]
  33. Bormetti, G.; Delpini, D. Exact moment scaling from multiplicative noise. Phys. Rev. E 2010, 81, 032102. [Google Scholar] [CrossRef] [Green Version]
  34. Delpini, D.; Bormetti, G. Minimal model of financial stylized facts. Phys. Rev. E 2011, 83, 041111. [Google Scholar] [CrossRef] [Green Version]
  35. Scalas, E.; Gorenflo, R.; Mainardi, F. Fractional calculus and continuous-time finance. Phys. A Stat. Mech. Appl. 2000, 284, 376–384. [Google Scholar] [CrossRef] [Green Version]
  36. Raberto, M.; Scalas, E.; Mainardi, F. Waiting-times and returns in high-frequency financial data: An empirical study. Phys. A Stat. Mech. Appl. 2002, 314, 749–755. [Google Scholar] [CrossRef] [Green Version]
  37. Magdziarz, M.; Gajda, J. Anomalous dynamics of Black–Scholes model time-changed by inverse subordinators. Acta Phys. Polonica B 2012, 43, 5. [Google Scholar]
  38. Wang, J.; Liang, J.R.; Lv, L.J.; Qiu, W.Y.; Ren, F.Y. Continuous time Black–Scholes equation with transaction costs in subdiffusive fractional Brownian motion regime. Physica A 2012, 391, 750–759. [Google Scholar] [CrossRef]
  39. Karipova, G.; Magdziarz, M. Pricing of basket options in subdiffusive fractional Black–Scholes model. Chaos Solitons Fractals 2017, 102, 245–253. [Google Scholar] [CrossRef]
  40. Gajda, J.; Wyłomańska, A. Geometric Brownian motion with tempered stable waiting times. J. Stat. Phys. 2012, 148, 296–305. [Google Scholar] [CrossRef]
  41. Leibovich, N.; Barkai, E. Infinite ergodic theory for heterogeneous diffusion processes. Phys. Rev. E 2019, 99, 042138. [Google Scholar] [CrossRef] [Green Version]
  42. Merton, R.C. Theory of rational option pricing. Bell J. Econ. Manag. Sci. 1973, 4, 141–183. [Google Scholar] [CrossRef] [Green Version]
  43. Hull, J.C. Options Futures and Other Derivatives; Pearson: Prentice Hall, NJ, USA, 2003. [Google Scholar]
  44. Kou, S.G. A jump-diffusion model for option pricing. Manag. Sci. 2002, 48, 1086–1101. [Google Scholar] [CrossRef] [Green Version]
  45. Mainardi, F. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models; World Scientific: Sigapore, 2010. [Google Scholar]
  46. Metzler, R.; Klafter, J. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep. 2000, 339, 1–77. [Google Scholar] [CrossRef]
  47. Fogedby, H.C. Langevin equations for continuous time Lévy flights. Phys. Rev. E 1994, 50, 1657. [Google Scholar] [CrossRef] [Green Version]
  48. Li, C. Option Pricing with Generalized Continuous Time Random Walk Models. Ph.D. Thesis, Queen Mary University of London, London, UK, 2016. [Google Scholar]
  49. Feller, W. An Introduction to Probability Theory and Its Applications; John Wiley & Sons.: Hoboken, NJ, USA, 2008; Volume 2. [Google Scholar]
  50. Magdziarz, M.; Weron, A.; Weron, K. Fractional Fokker-Planck dynamics: Stochastic representation and computer simulation. Phys. Rev. E 2007, 75, 016708. [Google Scholar] [CrossRef] [PubMed]
  51. Magdziarz, M.; Weron, A.; Klafter, J. Equivalence of the fractional Fokker-Planck and subordinated Langevin equations: The case of a time-dependent force. Phys. Rev. Lett. 2008, 101, 210601. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Garra, R.; Garrappa, R. The Prabhakar or three parameter Mittag–Leffler function: Theory and application. Commun. Nonlinear Sci. Numer. Simul. 2018, 56, 314–329. [Google Scholar] [CrossRef] [Green Version]
  53. Sornette, D.; Johansen, A.; Bouchaud, J.P. Stock market crashes, precursors and replicas. J. Phys. I 1996, 6, 167–175. [Google Scholar] [CrossRef] [Green Version]
  54. Mura, A.; Taqqu, M.S.; Mainardi, F. Non-Markovian diffusion equations and processes: Analysis and simulations. Phys. A Stat. Mech. Appl. 2008, 387, 5033–5064. [Google Scholar] [CrossRef] [Green Version]
  55. Mura, A.; Pagnini, G. Characterizations and simulations of a class of stochastic processes to model anomalous diffusion. J. Phys. A Math. Theor. 2008, 41, 285003. [Google Scholar] [CrossRef] [Green Version]
  56. Sposini, V.; Chechkin, A.V.; Seno, F.; Pagnini, G.; Metzler, R. Random diffusivity from stochastic equations: Comparison of two models for Brownian yet non-Gaussian diffusion. New J. Phys. 2018, 20, 043044. [Google Scholar] [CrossRef]
  57. Barkai, E. Fractional Fokker-Planck equation, solution, and application. Phys. Rev. E 2001, 63, 046118. [Google Scholar] [CrossRef] [Green Version]
  58. Meerschaert, M.M.; Benson, D.A.; Scheffler, H.P.; Baeumer, B. Stochastic solution of space-time fractional diffusion equations. Phys. Rev. E 2002, 65, 041103. [Google Scholar] [CrossRef] [Green Version]
  59. Schulz, J.H.; Barkai, E.; Metzler, R. Aging renewal theory and application to random walks. Phys. Rev. X 2014, 4, 011028. [Google Scholar] [CrossRef] [Green Version]
  60. Schilling, R.L.; Song, R.; Vondracek, Z. Bernstein Functions: Theory and Applications; Walter de Gruyter: Berlin, Germany, 2012; Volume 37. [Google Scholar]
  61. Sandev, T.; Metzler, R.; Chechkin, A. From continuous time random walks to the generalized diffusion equation. Fract. Calculus Appl. Anal. 2018, 21, 10–28. [Google Scholar] [CrossRef] [Green Version]
  62. Sandev, T.; Sokolov, I.M.; Metzler, R.; Chechkin, A. Beyond monofractional kinetics. Chaos Solitons Fractals 2017, 102, 210–217. [Google Scholar] [CrossRef]
  63. Prabhakar, T.R. A singular integral equation with a generalized Mittag Leffler function in the kernel. Yokohama Math. J. 1971, 19, 7–15. [Google Scholar]
  64. Sandev, T.; Chechkin, A.V.; Korabel, N.; Kantz, H.; Sokolov, I.M.; Metzler, R. Distributed-order diffusion equations and multifractality: Models and solutions. Phys. Rev. E 2015, 92, 042117. [Google Scholar] [CrossRef] [Green Version]
  65. Orzeł, S.; Mydlarczyk, W.; Jurlewicz, A. Accelerating subdiffusions governed by multiple-order time-fractional diffusion equations: Stochastic representation by a subordinated Brownian motion and computer simulations. Phys. Rev. E 2013, 87, 032110. [Google Scholar] [CrossRef] [Green Version]
  66. Bayer, C.; Friz, P.; Gatheral, J. Pricing under rough volatility. Quant. Financ. 2016, 16, 887–904. [Google Scholar] [CrossRef]
  67. Hull, J.; Treepongkaruna, S.; Colwell, D.; Heaney, R.; Pitt, D. Fundamentals of Futures and Options Markets; Pearson: London, UK, 2013. [Google Scholar]
  68. Cherstvy, A.G.; Vinod, D.; Aghion, E.; Chechkin, A.V.; Metzler, R. Time averaging, ageing and delay analysis of financial time series. New J. Phys. 2017, 19, 063045. [Google Scholar] [CrossRef] [Green Version]
  69. El Euch, O. Quantitative Finance under Rough Volatility. Ph.D. Thesis, Sorbonne Université, Paris, France, 2018. [Google Scholar]
  70. Jain, R.; Sebastian, K.L. Diffusion in a crowded, rearranging environment. J. Phys. Chem. B 2016, 120, 3988–3992. [Google Scholar] [CrossRef]
  71. Chechkin, A.V.; Seno, F.; Metzler, R.; Sokolov, I.M. Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities. Phys. Rev. X 2017, 7, 021002. [Google Scholar] [CrossRef] [Green Version]
  72. Wang, W.; Cherstvy, A.G.; Liu, X.; Metzler, R. Anomalous diffusion and nonergodicity for heterogeneous diffusion processes with fractional Gaussian noise. Phys. Rev. E 2020, 102, 012146. [Google Scholar] [CrossRef]
  73. Wang, W.; Seno, F.; Sokolov, I.M.; Chechkin, A.V.; Metzler, R. Unexpected crossovers in correlated random-diffusivity processes. New J. Phys. 2020, 22, 083041. [Google Scholar] [CrossRef]
  74. Sandev, T.; Iomin, A.; Ljupco, K. Hitting times in turbulent diffusion due to multiplicative noise. Phys. Rev. E 2020, 102, 042109. [Google Scholar] [CrossRef]
  75. Oberhettinger, F. Tables of Mellin Transforms; Springer Science & Business Media: Berlin, Germnay, 2012. [Google Scholar]
  76. Angstmann, C.N.; Donnelly, I.C.; Henry, B.I.; Langlands, T.; Straka, P. Generalized continuous time random walks, master equations, and fractional Fokker–Planck equations. SIAM J. Appl. Math. 2015, 75, 1445–1468. [Google Scholar] [CrossRef] [Green Version]
  77. Angstmann, C.N.; Donnelly, I.C.; Henry, B.I. Continuous time random walks with reactions forcing and trapping. Math. Model. Natural Phenomena 2013, 8, 17–27. [Google Scholar] [CrossRef] [Green Version]
  78. Sandev, T.; Chechkin, A.; Kantz, H.; Metzler, R. Diffusion and Fokker-Planck-Smoluchowski equations with generalized memory kernel. Fract. Calculus Appl. Anal. 2015, 18, 1006. [Google Scholar]
  79. Mathai, A.M.; Saxena, R.K.; Haubold, H.J. The H-Function: Theory and Applications; Springer Science & Business Media: Berlin, Germany, 2009. [Google Scholar]
Figure 1. Generalised geometric Brownian motion (gGBM) properties. (a) An example for simulated individual trajectories of gGBM for different memory kernels: standard GBM (blue solid line), subdiffusive GBM (sGBM) (red dashed line), a mix of standard GBM and sGBM (yellow dotted line), and a mix of sGBM (violet dot-dashed line). (b) Numerical estimation for the first moment in GBM, sGBM, a mix of standard GBM and sGBM and a mix of sGBM as a function of time. (c) Same as (b), only for the second moment. (d) Empirical PDF for the logarithmic return at t = 1 year estimated from 1000 realisations of gGBM. (ad) In the simulations, μ = 0.03 and σ 2 = 0.02 . Moreover, for the sGBM case we set α = 0.8 , for the mix GBM-sGBM case, we set α = 0.8 and w 1 = w 2 = 0.5 , and for the mix of sGBM case α 1 = 0.8 , α 2 = 0.6 , and w 1 = w 2 = 0.5 .
Figure 1. Generalised geometric Brownian motion (gGBM) properties. (a) An example for simulated individual trajectories of gGBM for different memory kernels: standard GBM (blue solid line), subdiffusive GBM (sGBM) (red dashed line), a mix of standard GBM and sGBM (yellow dotted line), and a mix of sGBM (violet dot-dashed line). (b) Numerical estimation for the first moment in GBM, sGBM, a mix of standard GBM and sGBM and a mix of sGBM as a function of time. (c) Same as (b), only for the second moment. (d) Empirical PDF for the logarithmic return at t = 1 year estimated from 1000 realisations of gGBM. (ad) In the simulations, μ = 0.03 and σ 2 = 0.02 . Moreover, for the sGBM case we set α = 0.8 , for the mix GBM-sGBM case, we set α = 0.8 and w 1 = w 2 = 0.5 , and for the mix of sGBM case α 1 = 0.8 , α 2 = 0.6 , and w 1 = w 2 = 0.5 .
Entropy 22 01432 g001
Figure 2. gGBM skewness and excess kurtosis. (a) Skewness for the distribution of logarithmic return for sGBM, a mix of standard GBM and sGBM and a mix of sGBM at t = 1 year as a function of α (for sGBM and the mix of GBM-sGBM) or α 1 (for the mix of sGBMs). (b) Same as (a), only for the excess kurtosis. (a,b) The PDF for the logarithmic return at t = 1 year is estimated from 1000 realisations of gGBM. In the simulations, μ = 0.03 and σ 2 = 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 .
Figure 2. gGBM skewness and excess kurtosis. (a) Skewness for the distribution of logarithmic return for sGBM, a mix of standard GBM and sGBM and a mix of sGBM at t = 1 year as a function of α (for sGBM and the mix of GBM-sGBM) or α 1 (for the mix of sGBMs). (b) Same as (a), only for the excess kurtosis. (a,b) The PDF for the logarithmic return at t = 1 year is estimated from 1000 realisations of gGBM. In the simulations, μ = 0.03 and σ 2 = 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 .
Entropy 22 01432 g002
Figure 3. gGBM at-the-money (ATM) implied volatility and volatility skew. (a) ATM implied volatility for sGBM, a mix of standard GBM and sGBM and a mix of sGBM for an option with K = x 0 = 1 , T = 0.083 years (one month) as a function of α (for sGBM and the mix of GBM-sGBM) or α 1 (for the mix of sGBMs). (b) Same as (a), only for the ATM volatility skew. (a,b) 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 .
Figure 3. gGBM at-the-money (ATM) implied volatility and volatility skew. (a) ATM implied volatility for sGBM, a mix of standard GBM and sGBM and a mix of sGBM for an option with K = x 0 = 1 , T = 0.083 years (one month) as a function of α (for sGBM and the mix of GBM-sGBM) or α 1 (for the mix of sGBMs). (b) Same as (a), only for the ATM volatility skew. (a,b) 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 .
Entropy 22 01432 g003
Figure 4. Moneyness in Tesla (TSLA). Root mean squared error (RMSE) of the predicted option price as a function of α ( α 1 ) for sGBM, a mix of GBM-sGBM and a mix of sGBM. The inset plot gives the difference between the predicted TSLA option price C g and its real value C as a function of the strike price of the option for various choices of α . The data is taken on 1st March 2020 and describe the value of TSLA options that expire on 16th March 2020. 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 .
Figure 4. Moneyness in Tesla (TSLA). Root mean squared error (RMSE) of the predicted option price as a function of α ( α 1 ) for sGBM, a mix of GBM-sGBM and a mix of sGBM. The inset plot gives the difference between the predicted TSLA option price C g and its real value C as a function of the strike price of the option for various choices of α . The data is taken on 1st March 2020 and describe the value of TSLA options that expire on 16th March 2020. 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 .
Entropy 22 01432 g004
Figure 5. Maturity in at-the-money AAPL options. Root mean squared error (RMSE) of the prediction of the AAPL at-the-money option price with data taken on 28th February 2018 as a function of α for various maturity periods T (measured in years). (a) The gGBM model is sGBM. (b) The gGBM model is a mix of GBM-sGBM. (c) The gGBM model is a mix of sGBM. 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 .
Figure 5. Maturity in at-the-money AAPL options. Root mean squared error (RMSE) of the prediction of the AAPL at-the-money option price with data taken on 28th February 2018 as a function of α for various maturity periods T (measured in years). (a) The gGBM model is sGBM. (b) The gGBM model is a mix of GBM-sGBM. (c) The gGBM model is a mix of sGBM. 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 .
Entropy 22 01432 g005
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 .
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 .
MaturitysGBMGBM-sGBMmix of sGBM
(in Years) min α RMSE min α RMSE min α 1 RMSE
0.0060.890.280.840.280.880.31
0.0250.960.390.930.390.930.48
0.0440.970.460.930.460.970.6
0.0630.990.520.980.520.980.71
0.1011.000.720.990.720.991.00
0.1401.000.610.990.600.990.91
0.8880.334.520.356.380.385.47
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Stojkoski, V.; Sandev, T.; Basnarkov, L.; Kocarev, L.; Metzler, R. Generalised Geometric Brownian Motion: Theory and Applications to Option Pricing. Entropy 2020, 22, 1432. https://doi.org/10.3390/e22121432

AMA Style

Stojkoski V, Sandev T, Basnarkov L, Kocarev L, Metzler R. Generalised Geometric Brownian Motion: Theory and Applications to Option Pricing. Entropy. 2020; 22(12):1432. https://doi.org/10.3390/e22121432

Chicago/Turabian Style

Stojkoski, Viktor, Trifce Sandev, Lasko Basnarkov, Ljupco Kocarev, and Ralf Metzler. 2020. "Generalised Geometric Brownian Motion: Theory and Applications to Option Pricing" Entropy 22, no. 12: 1432. https://doi.org/10.3390/e22121432

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop