VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA265 10.15559/24-VMSTA265 Research Article Subdiffusive option price model with Inverse Gaussian subordinator https://orcid.org/0000-0002-7652-8157 ShchestyukNataliyanataliya.shchestyuk@oru.seab∗ https://orcid.org/0000-0003-1804-8494 TyshchenkoSergiis.tyshchenko@ukma.edu.uab School of Business, Örebro University, Fakultetsgatan 1, Örebro, 701 82, Sweden Faculty of Informatics, National University of Kyiv-Mohyla Academy, Scovoroda Street, 2, Kyiv, 04070, Ukraine Corresponding author. 2025 1211202412213515216120246820243102024 © 2025 The Author(s). Published by VTeX2025 Open access article under the CC BY license.

The paper focuses on the option price subdiffusive model under the unusual behavior of the market, when the price may not be changed for some time, which is a quite common situation in modern illiquid financial markets or during global crises. In the model, the risk-free bond motion and classical geometrical Brownian motion (GBM) are time-changed by an inverted inverse Gaussian( IG) subordinator. We explore the correlation structure of the subdiffusive GBM stock returns process, discuss option pricing techniques based on the martingale option pricing method and the fractal Dupire equation, and demonstrate how it applies in the case of the IG subordinator.

 Option pricing subdiffusion models subordinator inverse subordinator time-changed process hitting time 91B24 91G20 91G60 Nataliya Shchestyuk acknowledges financial support from the project “Portfolio management for illiquid markets” (Dnr: 20220099) funded by the Knowledge Foundation. Introduction

Nowadays, modeling the financial market dynamics using diffusion processes has become an active research area in financial risk management, asset valuation, and derivatives pricing. However, the classical diffusion models, like the Black-Scholes-Merton (B-S) model and others based on Brownian motion (BM) have a normal distribution, independence of returns and are perpetually moving. These assumptions are inconsistent with empirical properties such as heavy-tailed and skewed marginal distributions, dependence on squared returns, and constant motionless periods (called also trapping events). The trapping events can be observed, for example, during global crises that negatively affect financial activity, and some types of risky assets have periods in their dynamics without changes. Such behavior is typical also for illiquid markets with a low number of transactions, interest rate markets, and commodity markets.

To overcome these difficulties, one can notice that the constant periods of stagnation in financial processes are analogous to the trapping events of the subdiffusive particle, therefore, the physical models of subdiffusion (a kind of anomalous diffusion, see [18]) can be successfully applied to describe financial data. Subdiffusion in the different physical systems arises from some memory effect of previous states, as a result of some fractal structure of the background space, or due to some nonlinear interactions inherent in the system, etc. To model subdiffusion, the time structure of the stochastic process is changed, and the time-changed process is no longer Markov. (Markov means each new step in the motion depends only on the present state and is independent of the previous states.)

The idea of the time-changed process was introduced in [3]. By time change we mean the replacement of the calendar (deterministic) time in a considered stochastic process S t using a new random clock. In our case, as a particular random clock, we consider the inverse H t of a subordinator G t . We recall that a subordinator is a nondecreasing Lévy process, i.e. it has stationary and independent increments and it is stochastically continuous. The inverse subordinator is used as a new random (“hitting”) time H t , thus as a studied model we obtain S H t . The distributional properties, asymptotic behavior, and the simulation procedures of the time-change process were considered for different types of subordinators and their inverses: α-stable, exponential, gamma, Pareto, Mittag-Leffler, and tempered distributions (see [22, 10, 14, 26, 2] and other).

However, the pricing of derivatives in this framework remains a complicated and under-researched problem. Nevertheless, more and more attention has been paid to this problem recently. In the papers [16, 17], the authors applied the subdiffusive mechanism of trapping events to describe financial data demonstrating periods of constant values and introduced the subdiffusive geometric BM and subdiffusive arithmetic BM. Ref. [16] showed that the considered models for α-subordinators are arbitrage-free but incomplete, and used the martingale option pricing method. A general result that implies that subdiffusive markets are generally incomplete is given in [24], Theorem 2. In [20], Theorem 5, a pseudodifferential equation that generalizes the Black and Scholes equation for the underlying time-changed geometric Brownian motion (GBM) was derived. The other technique for European option pricing was proposed in [7], where the derivative for the time in the Dupire equation is replaced by a Dzerbayshan–Caputo (D-C) derivative for jump-diffusion. The applications of anomalous diffusions for option pricing and the volatility term structure were considered in [9, 24]. Option pricing in time-changed Lévy models with compound Poisson jumps was explored in [8]. The theory in these papers was detailed for inverse α-stable ([16, 17, 7]), inverse tempered stable ([17, 24]), Mittag-Leffler ([20]) and inverted Poisson processes ([7]).

The paper considers the GBM model in a subdiffusion regime, time-changed by an inverted Inverse Gaussian ( IG) subordinator. We aim to show that the studied subdiffusive model demonstrates long-range dependence of stock returns and to discuss two different ways for option pricing. In addition, we propose a procedure for evaluating value-at-risk in the considered model. From this perspective, our study closes a gap which is of particular interest to investors.

The paper is organized as follows.

For introducing a subdiffusive model for the dynamics of a financial market we assume that the market consists of at least three components: one riskless asset B t , one risky asset with price S t , and one derivative security, usually called call option with price C t . The dynamics of the first two components in the subdiffusive framework are presented in the second section. We consider the IG process as the subordinator (waiting time) and the inverse Inverse Gaussian IIG process as the inverse subordinator (hitting time) for subdiffusive GBM, describe their properties and features, demonstrate the simulation for them, and for the subdiffusive GBM stock returns processes. We explore the correlation structure of the subdiffusive GBM stock returns process and assume that the stock returns process has long-range dependence and it is presented in the squared returns. Finally, we mentioned the Fractional Fokker–Planck equation (FFPE) for IG subordinator as the usual approach for modeling subdiffusion in physics. This equation describes the probability density function w ( t ) of the studied subdiffusive stock process, and we discuss how it can be used for risk measuring.

The next section focuses on the third component of the model and discusses two option pricing techniques. The one technique is very common for option pricing and can be given by the discounted expectation of the payoff conditional on the natural filtration up to time t concerning an equivalent martingale measure. The other technique is based on the fractal Dupire equation and uses the Dzerbayshan–Caputo (D-C) derivative. Since the form of the D-C derivative depends upon the chosen inverted Lévy subordinator, we demonstrate how it applies to the IG subordinator. The fourth section contains the numerical illustration for real financial data.

 Subdiffusive GBM model with <inline-formula id="j_vmsta265_ineq_017"><alternatives><mml:math> <mml:mi mathvariant="italic">IG</mml:mi></mml:math><tex-math><![CDATA[$\mathit{IG}$]]></tex-math></alternatives></inline-formula> subordinator

To build the mathematical model for the dynamics of a financial market first, we need to assume what kinds of securities evolve on the market and to describe their dynamics. Assume that the market consists of at least one riskless asset, usually called bond B t , a risky asset with price S t , usually called the stock, and one derivative security, usually called call option, or put option, which will have a certain payoff at a specified date in the future, depending on the values taken by the stock up to that date.

The idea is to replace the calendar time t in risk-free bond motion and classical GBM by some stochastic process H t , which will represent a stochastic clock or a stochastic operation time.

The time-changed risk-free bond has a value at time t equal to d B H t B H t = r d H t , B 0 = 1 , and the movement of the underlying risk assets S t follows a subdiffusive geometric Brownian motion (GBM) d S H t S H t = μ + σ 2 2 d H t + σ d B H t , t > 0 , with the solution S t = S 0 e μ H t + σ B H t , t > 0 . In formula (3) the standard diffusive process S t = S H t is time-changed by some stochastic process H t , which is called the inverse subordinator (“hitting time”). In general, for a given subordinator G t , the inverse subordinator H t is defined as H t = inf τ > 0 : G τ ≥ t , and interpreted as the first time at which G t hits the barrier t. H t is positive and nondecreasing and has all the properties to be used as a stochastic clock. By construction, the inverted process may be constant. Therefore, any process subordinated by H t exhibits motionless periods.

The definition (4) of the inverse subordinator is based on the use of some other random process called a subordinator G t . The subordinator G t is generally a nondecreasing stochastic process with stationary independent increments (Lévy process), taking values in R + and having Laplace transform E e − u G t = e − t Ψ ( u ) , where Ψ ( u ) is called Lévy exponent, which can be written as Ψ ( u ) = b u + ∫ 0 + ∞ 1 − e − u x ν ˜ ( d x ) . Here, b ≥ 0 is the drift parameter. If for simplicity, following [10], we assume b = 0, then ν ˜ ( d x ) is an appropriate Lévy measure.

The subordinator G t in our framework is often called the “waiting” time. In this paper, we consider the Inverse Gaussian process as a subordinator.

The IG subordinator and its inverse

Inverse Gaussian ( IG) subordinator G t is a nondecreasing Lévy process, where the increments G t + s − G s follow the inverse Gaussian distribution ϱ ( δ t , γ ) with probabilities density function (PDF) [22] g ( x , t ) = δ t 2 π x 3 e δ γ t − ( δ 2 t 2 / x + γ 2 x ) / 2 , x > 0 , and with the Lévy measure ν ˜ ( d x ) = δ 2 π x 3 e − γ 2 x 2 d x , x > 0 , t > 0 . For γ = δ = 1 we have the standard IG distribution in the form f ( x , t ) = t 2 π x 3 e − ( x − t ) 2 2 x , x > 0 , t > 0 . The Laplace transform of the IG subordinator is E e − u G t = e − t δ ( 2 u + γ 2 − γ ) , therefore, the Laplace exponent for IG process is given by Ψ IG ( u ) = δ ( 2 u + γ 2 − γ ) .

The tail probability for the IG subordinator is of the form ([26]) P ( G t > x ) ∼ 2 π δ t γ 2 e γ δ t x − 3 / 2 e − ( γ 2 / 2 ) x , x → ∞ , where f ( x ) ∼ g ( x ) as x → 0 means lim x → 0 f ( x ) g ( x ) = 1. Notice, that the tail probability decreases exponentially for γ ≠ 0. Thus all moments of the process G t are finite.

The qth moments for the IG process are found in the paper [11]: E ( G t q ) = 2 π δ δ γ q − 1 / 2 t q + 1 / 2 e δ γ t K q − 1 / 2 ( δ γ t ) , where K q ( ω ) is the modified Bessel function of the third kind with index q. Moreover, in [11] it was shown that if t → ∞, then E ( G t q ) ∼ δ γ q t q . For the standard distribution, we have E ( G ( t ) ) ∼ t, var ( G ( t ) ) ∼ t 2 .

The algorithm of the simulation of the IG process G t for time points t 1 = 1 n , t 2 = 2 n , … , t n = 1 was proposed in some literature (see, for example, [26, 22]). Since the process G t has independent and stationary increments, F i = G t i − G t i − 1 = G d t ∼ ϱ ( d t , 1 ) for i = 1 , 2 , … n and d t = 1 n , so following [26] we generate n i.i.d. IG variables F i assuming γ = δ = 1. For this, we simulate a standard normal random variable N and a uniform [ 0 , 1 ] random variable U. Then assign X = N 2 and Y = d t + X 2 − 1 2 4 d t X + X 2 . According this algorithm, if U ≤ d t d t + Y then return Y; otherwise return ( d t ) 2 Y . After assigning G t 0 = 0 and G t i = ∑ j = 1 i F j , i = 1 , 2 , … , n, we obtain G t 1 , G t 2 , … , G t n , which are simulated values of the IG process at times t 1 , t 2 , … , t n , respectively. The simulated trajectories of the standard IG process are shown in Figure 1.

Simulation of the IG process trajectories for γ = δ = 1

Now we discuss the properties of the inverse subordinator (hitting time) H t defined by (4). The inverse to the inverse Gaussian ( IIG) process is not Lévy and has monotonically increasing continuous sample paths. Moreover, the sample paths of the IIG process are constant over the intervals where G t has jumped. It follows from the fact that the trajectories of G t are strictly increasing with jumps.

The distribution of H t is not infinitely divisible ([26]).

The density function h ( x , t ) of H t can be put in the integral form ([25]): h ( x , t ) = δ π e δ γ x − γ 2 2 ∫ 0 ∞ e − t y y + γ 2 2 ( γ sin ( δ x 2 y ) + 2 y cos ( δ x 2 y ) ) d y . It is worth noticing that if γ = 0 then Laplace exponent (8) has the form Ψ γ = 0 ( u ) = ( δ 2 ) u 1 / 2 , and it means that we have the stable distribution TS ( 1 / 2 , 0 , δ ). In the case γ ≠ 0 we have the tempered stable distribution TS ( 1 / 2 , γ 2 , δ ). For the particular case γ = 0, δ = 1 we have 1 / 2 −stable distribution S ( α ), α = 1 / 2, and the density function for the inverse IG process is equal h ( x , t ) = 2 π t e − x 2 2 t . It should be mentioned that the expression (10) has been extended to general inverse tempered stable subordinators in [13], and the case of a stable subordinator can be deduced as a limit. So, the explicit expression (10) for the hitting time allows us to find the closed form for a tempered 1 / 2-stable subordinator and 1 / 2-stable process as particular cases. Moreover, the closed form of the density distribution of hitting time will be used to compute option pricing in the next section.

The qth moments for the IIG process may be numerically evaluated for known t by using the density function h ( x , t ). However, an explicit expression for the first and second order was obtained by using the Laplace transformation in [25].

Moreover, in [11] you can find the following result regarding the asymptotic behavior of E ( H t ) and Var ( H t ). If t → ∞, then E ( H t ) ∼ γ δ t , γ > 0 , 1 δ 2 t π t , γ = 0 , Var ( H t ) ∼ γ δ 2 t 2 . For the standard distribution we have E ( H t ) ∼ t, v a r ( H t ) ∼ t 2 .

To simulate the approximate trajectory of the inverse subordinator H t , we define H t Δ as follows [26]: H t Δ = [ min { n ∈ N : G Δ n > t } − 1 ] Δ , n = 1 , 2 , … , where Δ is the step length and G Δ n is the value of the Inverse Gaussian process G t at Δ n. The simulation of the trajectory H t for γ = δ = 1 is demonstrated in Figure 2.

Simulation of the inverse to the IG process trajectories for γ = δ = 1

For modeling of the stochastic subdiffusive GBM we propose the iterative scheme S t + 1 = S t + μ S t H t Δ + σ S t H t Δ ε t , where ε is a white noise with the normal standard distribution, and H t Δ are given in (13). The trajectories of the subdiffusion GBM with the inverse to the IG process are demonstrated in Figure 3.

Simulation of the subdiffusive GBM time-changed by an inverse IG subordinator

 Correlation structure for log returns of time-changed GBM

Exploring the dependence structure for stock processes or their stock returns is an important issue in the study of diffusive and subdiffusive models (see, for example, [9, 26, 14] and [4]).

First, we compute log S H t = log S 0 + μ H t + σ B H t ( 1 ) , where S H t is a stochastic process given by (3).

Then denote τ t = H t − H t − 1 and obtain the stock returns in the form X t = log S H t S H t − 1 = μ τ t + σ τ t B 1 ( 1 ) , using the scaling law of the Brownian motion.

Let X t be a stochastic process given by (15), then for any integer k ≥ 0:

Cov ( X t , X t + k ) = μ 2 Cov ( τ t , τ t + k ) ; and stock returns are uncorrelated for μ = 0: Cov μ = 0 ( X t , X t + k ) = 0.

There is a long-range dependence in the squared returns: Cov ( X t 2 , X t + k 2 ) = σ 4 Cov ( τ t , τ t + k ) + μ 4 Cov ( τ t 2 , τ t + k 2 ) + 2 μ 2 σ 2 Cov ( τ t 2 , τ t + k ) ; in particular, Cov μ = 0 ( X t 2 , X t + k 2 ) = σ 4 Cov ( τ t , τ t + k ).

For the covariance function of the log returns process X t (15), we obtain Cov ( X t , X t + k ) = Cov ( μ τ t + σ τ t B 1 ( 1 ) , μ τ t + k σ τ t + k B 1 ( 2 ) ) = E μ τ t + σ τ t B 1 ( 1 ) − μ E [ τ t ] μ τ t + k + σ τ t + k B 1 ( 2 ) − μ E [ τ t + k ] = μ 2 E [ τ t τ t + k ] − E [ τ t ] E [ τ t + k ] = μ 2 Cov ( τ t , τ t + k ) , where the second and third equalities follow from the independence between the two stochastic processes.

The covariance function of the process X t 2 can be computed as Cov ( X t 2 , X t + k 2 ) = Cov ( μ τ t + σ τ t B 1 ( 1 ) ) 2 , ( μ τ t + k + σ τ t + k B 1 ( 2 ) ) 2 = σ 4 Cov ( τ t , τ t + k ) + μ 4 Cov ( τ t 2 , τ t + k 2 ) + μ 2 σ 2 Cov ( τ t 2 , τ t + k ) + μ 2 σ 2 Cov ( τ t + k , τ t 2 ) , where from the last expression we immediately get the result.  □

Thus, if the increments τ t of the subordinated process H t exhibit long-range dependence, the same holds for the log returns of X t 2 and we have long-range dependence for stock process (3) in the squared returns.

 Fractional Fokker–Planck equation and risk measuring for subdiffusion

The usual model of subdiffusion in physics is the celebrated Fractional Fokker–Planck equation (see, for example, [18]). This equation was derived from the continuous time random walk scheme with heavy-tailed waiting times and describes the probability density function w ( t ) of the subdiffusive studied stock process: ∂ w ∂ t = Φ t − μ ∂ ∂ x + σ 2 2 ∂ 2 ∂ x 2 w ( x , t ) , where Φ t is the integro-differential operator defined as Φ t f ( t ) = d d t ∫ 0 t M ( t − y ) f ( y ) d y , with the memory kernel M ( t ) defined via its Laplace transform M ˜ ( u ) = ∫ 0 ∞ e − u t M ( t ) d t = 1 Ψ ( u ) , where Ψ ( u ) is the Laplace exponent of the IG subordinator and can be written as in (8). This implies that the memory kernel M ( t ) can be expressed as M ( t ) = L − 1 1 δ ( 2 u + γ 2 − γ ) , where L − 1 ( f ) is the inverse Laplace transform of the function f ( t ), which can be computed as in [10]. Thus, the formula (18) allows us to find, at least in some particular cases of parameters γ, δ, closed-form formulas for the PDF of the studied subdiffusive stock processes. In general, approximated solutions w ( t ) of (18) can be derived by the finite element method for FFPE (see, for example, [6]) or by the Monte Carlo techniques based on the simulation algorithm of the time-changed stock process (see the section above).

Thus, the possibility of numerical computing of the probability density function w ( t ) for the studied subdiffusive stock process (with IG subordinator) opens the way to evaluate the value-at-risk (VaR) in this model.

The value-at-risk is a quite useful tool for investors and can be used for understanding the past and making medium-term and strategic decisions for the future. On the other side, we can apply VaR for checking of the model performance. For this, we can use the most important criterion of a risk management system, namely, to check if the regulatory requirements are fulfilled.

VaR can be defined as α-quantile of the profit (loss) function.

Let ( Ω , F , P ) be the probability space. The value-at-risk of level α, 0 < α ≤ 1, is a probability functional defined as α-quantile of the profit (loss) function Y ∈ L ( Ω ): VaR α ( Y ) = W − 1 ( α ) = inf { y ∈ R : α ≤ W ( y ) } , where W is the distribution function of Y, W − 1 is the quantile function of α, 0 < α ≤ 1.

Let the time horizon coincide with the time to maturity, then the loss (profit) function of the call option with strike price K is Y = Y ( S ) = | S − K | + − c 0 , and the value-at-risk of level α, 0 < α ≤ 1, for the random variable Y is VaR α ( Y ) = VaR α ( | S − K | + ) − c 0 due to the translation-equivariant property of the probability functional VaR.

The cumulative distribution function (CDF) for Y ( S ) = | S − K | + is given [21] by W Y ( y ) = ∫ − ∞ y w S ( u + K ) d u , y ≥ 0 , 0 , y < 0 .

Thus, if the time horizon coincides with the time to maturity, one can find the value-at-risk of level α, 0 < α ≤ 1, as (21)–(22), where w ( t ) is a solution of (18) with memory kernel (19).

 Option pricing

Even if subdiffusions can be successfully applied for modeling illiquidity, option pricing in this framework remains insufficiently researched.

Recall that in the classical GBM model the fair price of the European call option is given by the Black–Scholes formula C ( S , K , T , r , σ ) = N ( d 1 ) S − N ( d 2 ) K e − r T , where d 1 = log S 0 K + r T + 1 2 σ 2 T σ T , d 2 = log S 0 K + r T − 1 2 σ 2 T σ T both are functions of five parameters T , K , S 0 , r , σ, and Φ ( · ) is the standard normal cumulative distribution function.

Consider a time-changed version of the B-S model, where the price of the bond evolves as (1) and the underlying risky assets follow (2) with IG subordinator. Now we discuss two approaches to option pricing in the proposed model with inverse Gaussian subordinator.

Let the evolution of this market up to time horizon T be contained in the probability space ( Ω , F , P ). Here, Ω is the sample space, F contains all statements that can be made about the behavior of prices, and P is the “objective” probability measure. We denote by ( F t ) , t ∈ [ 0 , T ], the information about the history of asset prices S H up to time t. ( F t ) is also called filtration and is interpreted as the background information that is available for the investor. The more time proceeds the more information is revealed to the investor.

The European call option price C H = C H S , K , T , σ for a time-changed version of the B-S model with IG subordinator satisfies C H = ∫ 0 ∞ ∫ 0 ∞ C ( x ) δ π e δ γ x − γ 2 2 e − T y y + γ 2 2 γ sin ( ρ ( x , y ) ) + 2 y cos ( ρ ( x , y ) ) d y d x , where ρ ( x , y ) = δ x 2 y , and C ( T ) = C ( S , K , T , σ ) is given by (23).

For the subdiffusion market described above the usual requirement for fair option pricing is that arbitrage opportunities do not exist. For this, it is enough to prove the existence of the equivalent martingale measure. In [16], the measure Q ( A ) = ∫ A exp − τ B H T − τ 2 2 H T d P was introduced, where τ = μ + σ 2 σ and A ∈ F, and it was proved that subdiffusive GBM (3) with α-stable subordinator is a martingale with respect to Q. The Theorem 1 in Magdziarz and Schilling (2015) extend the results for the α-stable case to any case of the Lévy distribution. Consequently, using the Fundamental theorem of asset pricing (see [22]) we can state that the subdiffusion market model with any Lévy subordinator is arbitrage-free.

The second question is the completeness of the subdiffusive market model. In the mathematical theory of arbitrage, the model is complete if and only if there is a unique martingale measure. In Theorem 2 of [24] the whole family of possible transformations of equivalent martingale measures was constructed, which shows that the measure Q is not unique, and therefore the subdiffusive markets are generally incomplete. Particularly it applies to inverse-tempered stable subordinators. Thus, the subdiffusive IIG market is incomplete, and it is natural because it takes into account the presence of a market risk of trade duration.

Then, for arbitrage-free and incomplete IIG market, we apply a common technique for time-changed processes (see [16, 17, 4, 8]) and find the corresponding fair price of the European call option written on the time-changed asset as the discounted expected payoff under measure Q (26): C H S , K , T , σ = C S , K , H T , σ = E Q e − r H T ( S H T − K ) + | F 0 . Therefore, conditioning on H T , we obtain C H S , K , T , σ = ∫ 0 ∞ C S , K , x , σ h Ψ ( x , T ) d x . Here, h Ψ ( x , T ) is the PDF of H T for the subordinator with Lévy exponent Ψ and C ( S , K , T , σ ) is given by (23).

If we note that h Ψ ( x , T ) is PDF of I G G process for the considered model (10), then this ends the proof of formula (25).  □

The presented option pricing subdiffusive model with I G ( δ , γ ) subordinator has an interesting interpretation of its parameters in light of the market price of risk. Thanks to the theorem about the incompleteness of the subdiffusive market, the existence of a market price of duration risk naturally arises (see [24]). Once S t is calibrated to liquid market prices, the market price of duration risk will be reflected in the martingale density parameters, and these same features will be thus reflected in the (risk-adjusted) frequency and duration of trade pauses incorporated in hitting time H t .

Since IG subordinator with parameters δ and γ can be considered as tempered stable Lévy subordinator TS ( 1 / 2 , γ 2 , δ ), we can refer to the paper [24], where the market price of risk was explored for the class of tempered subdiffusive models. From this paper, we can state that the parameter γ 2 expresses belief in the speed at which the granular price evolution of the asset will revert to a fully “liquid” state, which is well approximated by a standard Lévy-driven diffusion. So, γ captures the risk of latency to liquidity, while absolute trade duration risk remains constant and is defined by 1 / 2. Moreover, if γ = 0 we get a regime where prices exhibit maximum staleness at all temporal scales. Therefore, IG models offer a unified framework in which ultrashort-term option pricing accounts for the microstructural duration effect, whereas options at typical maturities are priced according to the usual Lévy paradigm: the cut-off between the two states is encoded in γ (see [24] for more details).

In the above equation (27) we can evaluate the subdiffusive call price C ( · ) by computing the integral numerically. An alternative consists of calculating the price by Monte Carlo simulations. First, one simulates the trajectories for the inverse subordinator on the interval [ 0 , T ] by the approximation scheme (13). Then, one obtains the fair price as an estimation of the expected value for simulated prices where the inverse subordinator stands for calendar time T in (27) C H ( S , K , T , r , σ ) = ⟨ C ( S , K , H T , σ ) ⟩ = 1 n ∑ i = 1 n C ( S , K , H T ( i ) , σ ) , where C ( S , K , T , σ ) is taken from the Black–Scholes option pricing formula (23).

It is worth noticing, that the application of the Monte Carlo method for option pricing in subdiffusive models can be seen in the papers [16, 17] for α-stable subordinator.

The second approach to option pricing in the considered model is much more interesting and based on a fractional version of what is called Dupire’s equation. Dupire has established a forward partial differential equation for call options with local volatility. The fractional Dupire’s equation was proposed by [7] and it is integro-differential equation in partial derivatives. (PIDE). This equation was presented in a very general form and valid for all invertible Lévy subordinators. In the PIDE the derivative with respect to time was replaced by a convolution-type derivative, called the Dzerbayshan–Caputo (D-C) derivative. D-C derivative is a kind of fractional derivative, which is more advantageous than its classical counterparts due to capturing the past history. The Dzerbayshan–Caputo derivative depends upon the chosen kind of subordinator and its Lévy exponent Ψ ( . ), which is often called the Bernstein function.

The next proposition is the application of this PIDE to B-S subdiffusion with IG subordinators.

The European call option price C H ( T , k ) = C H S , K , T , σ for a time-changed version of the B-S model with IG subordinator is a solution of a fractional PIDE equation ∫ 0 T δ ∂ ∂ t C H ( T − s , k ) γ ( Φ γ s − 1 ) + e − s γ 2 2 2 π s d s = − r 2 ∂ ∂ k C H ( T , k ) + σ 2 4 ∂ 2 ∂ k 2 C H ( T , k ) , with the initial condition C H ( 0 , K ) = ( S H 0 − K ) + , and k = ln K. The European put option price is a solution of the same fractional PIDE equation but with initial condition P H ( 0 , K ) = ( K − P H 0 ) + .

In [7], a fractional PIDE equation for option pricing was presented in the fractional jump-diffusion setting. In the Black and Scholes (B-S) regime, when the Brownian volatility is constant and there are no jumps, the fractional Dupire equation has the form Ψ D C H ( T , k ) = − r ∂ ∂ k C H ( T , k ) + σ 2 2 ∂ 2 ∂ k 2 C H ( T , k ) , where Ψ D u ( t ) is the convolution-type derivative, called the Dzerbayshan–Caputo (D-C) derivative. The generalized D-C derivative according to the function Ψ is defined as (see, for example, [23]) Ψ D u ( t ) = b d d t u ( t ) + ∫ 0 t ∂ ∂ t u ( t − s ) ν ( s ) d s . The function Ψ is the Lévy exponent for a given subordinator G t , and we use its Lévy–Khintchine representation Ψ ( x ) = ∫ 0 + ∞ ( 1 − e − s z ) ν ˜ ( d z ) , where ν ˜ is the Lévy measure for this subordinator.

For IG Lévy subordinator the Lévy measure ν ˜ ( s ) is defined by (7). The integral kernel ν ( s ) is the integral of ν ˜ over ( s , ∞ ): ν ( s ) = ∫ s + ∞ δ 2 π x 3 e − γ 2 x 2 d x = δ γ 2 π ∫ s γ 2 2 + ∞ 1 z 3 e − z d z = δ γ 2 π Γ ( − 1 / 2 , s γ 2 2 ) , where Γ ( a , b ) is the upper incomplete Gamma function. If we use the fact that Γ ( a + 1 , b ) = a Γ ( a , b ) + b a e − b and expression for Γ ( 1 / 2 , b ), then D-C derivative (30) for IG subordinator can be written as Ψ D u ( t ) = δ γ ∫ 0 t ∂ ∂ t u ( t − s ) erf γ 2 s 2 − 1 + 2 e − s γ 2 2 γ 2 π s d s , where erf ( . ) is the error function, which can be expressed in terms of the standard normal cumulative distribution function Φ ( · ): Ψ D u ( t ) = 2 δ ∫ 0 t ∂ ∂ t u ( t − s ) γ ( Φ γ s − 1 ) + e − s γ 2 2 2 π s d s . Thus, if we substitute (31) in (29), the proposition follows.  □

We can notice that the derivative we obtain is a tempered D-C derivative (see [19, 20]). In general, these kinds of convolution operators were introduced in [12] and then further explored in [23, 1] and [5]. To solve (29), we need to extend a finite difference approach explored in [6] to the fractional kind of Dupire equation with IG subordinator.

 Numerical illustration

We illustrate the call option pricing for the subdiffusive IIG market.

We choose Airbnb, Inc. Class A Common Stock listed on the NASDAQ stock exchange and use their daily returns during the last two years for model calibration. All the prices and relative data we take on 24 June 2022. We fixed the strike price as 100 and observed the options prices with time to maturity ranging from the 1st of July 2022 to the 20th of January 2023. Market parameters are: S 0 = 103.51, K = 100, r = 0.168, μ = − 0.0015.

To simplify the calculations, we put parameters γ = δ = 1 because this assumption allows getting the desired result regarding the asymptotic behavior of the hitting time. In this case, when t → ∞, it follows that E ( H t ) ∼ t, var ( H t ) ∼ t 2 (11)–(12). So, we estimate only the σ parameters based on the least squares technique. We obtain that the value σ = 0.3 is the best to minimize the square error over the difference between the real option quotes and the estimated ones.

We estimated the prices (28) of call options using Monte Carlo methods based on the above-described simulation procedure for H T .

The results are presented in Table 1 and in graphic shape, where we compare the B-S subdiffusive European call options with the classical one and with the market price.

Simulated prices for the diffusive and subdiffusive B-S models

As we can see from the graphics in Figure 4, the diffusive option pricing model shows better results in the short-term period, while the subdiffusive model is more effective in the long-term perspective.

Call option prices for strike price K = 100 and varying times of maturity

1 Jul 8 Jul 15 Jul 22 Jul 29 Jul 19 Aug 16 Sep 21 Oct 16 Dec 20 Jan
Market 5.60 6.88 8.10 9.25 9.95 12.75 14.72 16.15 19.50 21.30
B-S diff 6.02 7.47 8.65 9.70 10.64 13.09 15.83 18.77 22.83 25.08
B-S subdiff 7.40 8.37 9.51 9.91 11.32 12.72 14.69 15.61 18.11 18.91

To compare numerical results we use absolute relative percentage (ARPE) and root mean squared error (RMSE): ARPE = | C ( T k ) − C market ( T k ) | C market ( T k ) RMSE = 1 n Σ i = 1 n ( C ( T i ) − C market ( T i ) C market ( T i ) ) 2

It is worth mentioning, that in econometrics, the root mean squared error (RMSE) (32) is a key criterion for model selection. The mean squared error indicates the mean squared deviation between the forecast and the outcome. It sums the squared bias and the variance of the estimator. The advantage of the ARPE (33) relative to the RMSE measure is that it gives a percentage value of the pricing error. Therefore, if we use both these errors it provides more insight into the economic significance of performance differences.

Upon dividing the RMSE analysis into two distinct periods, a noteworthy observation emerges: during the initial short period, the diffusive model outperforms, while in the subsequent long period, the subdiffusive model exhibits superior performance.

The ARP errors for B-S diffusion, B-S subdiffusion to the market price

1 Jul 8 Jul 15 Jul 22 Jul 29 Jul 19 Aug 16 Sep 21 Oct 16 Dec 20 Jan Mean
B-S 0.08 0.09 0.07 0.05 0.07 0.03 0.08 0.16 0.17 0.18 0.10
B-S Subdiffusion 0.32 0.22 0.17 0.07 0.14 0.00 0.00 0.03 0.07 0.11 0.11

The RMS errors for diffusion and subdiffusion regarding the market price

1 Jul – 29 Jul 19 Aug – 20 Jan Overall
B-S 0.07 0.14 0.11
B-S Subdiffusion 0.20 0.06 0.15

In the framework of the paper, we just illustrate the application of the studied model and compare option pricing results in a situation when strike price K was fixed (in the money), while time to maturity T was changing.

For detailed model performance we need to examine the ARP pricing errors of the proposed option pricing models in more detail (see [15]) and consider the pricing errors as a regression on the time to maturity T (in years), the moneyness of the option, and a binary variable that is set to unity, if the option is a call and to zero in the case of a put. This can indicate a level of explanatory value of moneyness, maturity, and the put-call dummy in the model.

However, this model performance and application of a finite difference approach for solving the fractional Dupire equation is a future work beyond the scope of the present paper.

 Conclusion

The paper developed a subdiffusive model with the following features of stock returns, which are quite well documented in the financial and econometric literature: (i) the stochastic processes have continuous paths with motionless periods of time; (ii) the returns processes are uncorrelated; (iii) dependence is presented in squared returns; (iv) the hitting time is defined by the inverse to IG subordinator which is not infinitely divisible. For the model, two option pricing techniques were discussed, a procedure for evaluating value-at-risk was proposed. The results of comparing the studied model with the classical one show that the classical B-S model demonstrates better results in the short-term period, while the subdiffusive model is more effective in the long-term perspective.

Thanks to the proposed approaches, the investor gets tools that allow him to take into account the market’s illiquidity.

 References Ascione, G.: Abstract cauchy problems for the generalized fractional calculus. Nonlinear Anal. 209, 112339 (2021) MR4236481. https://doi.org/10.1016/j.na.2021.112339 Beghin, L., Macci, C., Martinucci, B.: Random time-changes and asymptotic results for a class of continuous-time markov chains on integers with alternating rates. Mod. Stoch. Theory Appl. 8(1), 63–91 (2020) MR4235564. https://doi.org/10.15559/20-vmsta169 Bochner, S.: Diffusion equation and stochastic processes. Proc. Natl. Acad. Sci. 35(7), 368–370 (1949) MR0030151. https://doi.org/10.1073/pnas.35.7.368 Castelli, F., Leonenko, N., Shchestyuk, N.: Student-like models for risky asset with dependence. Stoch. Anal. Appl. 35(3), 452–464 (2017) MR3609373. https://doi.org/10.1080/07362994.2016.1266945 Chen, Z.-Q.: Time fractional equations and probabilistic representation. Chaos Solitons Fractals 102, 168–174 (2017) MR3672008. https://doi.org/10.1016/j.chaos.2017.04.029 Deng, W.: Finite element method for the space and time fractional fokker–planck equation. SIAM J. Numer. Anal. 47(1), 204–226 (2009) MR2452858. https://doi.org/10.1137/080714130 Hainaut, D., Leonenko, N.: Option pricing in illiquid markets: A fractional jump–diffusion approach. J. Comput. Appl. Math. 381, 112995 (2021) MR4107797. https://doi.org/10.1016/j.cam.2020.112995 Ivanov, R.V., Ano, K.: Option pricing in time-changed lévy models with compound poisson jumps. Mod. Stoch. Theory Appl. 6(1), 81–107 (2018) MR3935428. https://doi.org/10.15559/18-vmsta124 Jacquier, A., Torricelli, L.: Anomalous diffusions in option prices: connecting trade duration and the volatility term structure. SIAM J. Financ. Math. 11(4), 1137–1167 (2020) MR4173215. https://doi.org/10.1137/19M1289832 Janczura, J., Wyłomańska, A.: Anomalous diffusion models: different types of subordinator distribution. arXiv preprint arXiv:1110.2868 (2011) Jorgensen, B.: Statistical Properties of the Generalized Inverse Gaussian Distribution vol. 9. Springer (2012) MR0648107 Kochubei, A.N.: General fractional calculus, evolution equations, and renewal processes. Integral Equ. Oper. Theory 71(4), 583–600 (2011) MR2854867. https://doi.org/10.1007/s00020-011-1918-8 Kumar, A., Vellaisamy, P.: Inverse tempered stable subordinators. Stat. Probab. Lett. 103, 134–141 (2015) MR3350873. https://doi.org/10.1016/j.spl.2015.04.010 Kumar, A., Wyłomańska, A., Połoczański, R., Sundar, S.: Fractional brownian motion time-changed by gamma and inverse gamma process. Physica A 468, 648–667 (2017) MR3580705. https://doi.org/10.1016/j.physa.2016.10.060 Lehar, A., Scheicher, M., Schittenkopf, C.: Garch vs. stochastic volatility: Option pricing and risk management. J. Bank. Finance 26(2-3), 323–345 (2002) Magdziarz, M.: Black-scholes formula in subdiffusive regime. J. Stat. Phys. 136, 553–564 (2009) MR2529684. https://doi.org/10.1007/s10955-009-9791-4 Magdziarz, M., Orzeł, S., Weron, A.: Option pricing in subdiffusive bachelier model. J. Stat. Phys. 145, 187–203 (2011) MR2841938. https://doi.org/10.1007/s10955-011-0310-z Metzler, R., Klafter, J.: The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339(1), 1–77 (2000) MR1809268. https://doi.org/10.1016/S0370-1573(00)00070-3 Sabzikar, F., Meerschaert, M.M., Chen, J.: Tempered fractional calculus. J. Comput. Phys. 293, 14–28 (2015) MR3342453. https://doi.org/10.1016/j.jcp.2014.04.024 Scalas, E., Toaldo, B.: Limit theorems for prices of options written on semi-markov processes. Theory Probab. Math. Stat. 105, 3–33 (2021) MR4421361. https://doi.org/10.1090/tpms Shchestyuk, N., Tyshchenko, S.: Option pricing and stochastic optimization. In: International Conference on Stochastic Processes and Algebraic Structures, pp. 651–665 (2019). Springer MR4607864. https://doi.org/10.1007/978-3-031-17820-7_28 Tankov, P.: Financial Modelling with Jump Processes. Chapman and Hall/CRC (2003) Toaldo, B.: Convolution-type derivatives, hitting-times of subordinators and time-changed c 0-semigroups. Potential Anal. 42, 115–140 (2015) MR3297989. https://doi.org/10.1007/s11118-014-9426-5 Torricelli, L.: Trade duration risk in subdiffusive financial models. Physica A 541, 123694 (2020) MR4043912. https://doi.org/10.1016/j.physa.2019.123694 Vellaisamy, P., Kumar, A.: The hitting time of an inverse gaussian process. Preprint (2011) Wyłomańska, A., Kumar, A., Połoczański, R., Vellaisamy, P.: Inverse gaussian and its inverse process as the subordinators of fractional brownian motion. Phys. Rev. E 94(4), 042128 (2016)