1 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 ($\mathit{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 $\mathit{IG}$ process as the subordinator (waiting time) and the inverse Inverse Gaussian $\mathit{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 $\mathit{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 $\mathit{IG}$ subordinator. The fourth section contains the numerical illustration for real financial data.
2 Subdiffusive GBM model with $\mathit{IG}$ 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
and the movement of the underlying risk assets ${S_{t}}$ follows a subdiffusive geometric Brownian motion (GBM)
with the solution
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
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.
(2)
\[ \frac{d{S_{{H_{t}}}}}{{S_{{H_{t}}}}}=\left(\mu +\frac{{\sigma ^{2}}}{2}\right)d{H_{t}}+\sigma d{B_{{H_{t}}}},\hspace{1em}t\gt 0,\]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
where $\Psi (u)$ is called Lévy exponent, which can be written as
Here, $b\ge 0$ is the drift parameter. If for simplicity, following [10], we assume $b=0$, then $\tilde{\nu }(dx)$ 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.
2.1 The IG subordinator and its inverse
Inverse Gaussian ($\mathit{IG}$) subordinator ${G_{t}}$ is a nondecreasing Lévy process, where the increments ${G_{t+s}}-{G_{s}}$ follow the inverse Gaussian distribution $\varrho (\delta t,\gamma )$ with probabilities density function (PDF) [22]
For $\gamma =\delta =1$ we have the standard $\mathit{IG}$ distribution in the form
\[ g(x,t)=\frac{\delta t}{\sqrt{2\pi {x^{3}}}}{e^{\delta \gamma t-({\delta ^{2}}{t^{2}}/x+{\gamma ^{2}}x)/2}},\hspace{1em}x\gt 0,\]
and with the Lévy measure
(7)
\[ \tilde{\nu }(dx)=\frac{\delta }{\sqrt{2\pi {x^{3}}}}{e^{\left(-\frac{{\gamma ^{2}}x}{2}\right)}}dx,\hspace{1em}x\gt 0,\hspace{2.5pt}t\gt 0.\]
\[ f(x,t)=\frac{t}{\sqrt{2\pi {x^{3}}}}{e^{\left(-\frac{{(x-t)^{2}}}{2x}\right)}},\hspace{1em}x\gt 0,\hspace{2.5pt}t\gt 0.\]
The Laplace transform of the $\mathit{IG}$ subordinator is
therefore, the Laplace exponent for $\mathit{IG}$ process is given by
The tail probability for the $\mathit{IG}$ subordinator is of the form ([26])
\[ P({G_{t}}\gt x)\sim \sqrt{\frac{2}{\pi }}\frac{\delta t}{{\gamma ^{2}}}{e^{\gamma \delta t}}{x^{-3/2}}{e^{-({\gamma ^{2}}/2)x}},\hspace{1em}x\to \infty ,\]
where $f(x)\sim g(x)$ as $x\to 0$ means ${\lim \nolimits_{x\to 0}}\frac{f(x)}{g(x)}=1$. Notice, that the tail probability decreases exponentially for $\gamma \ne 0$. Thus all moments of the process ${G_{t}}$ are finite.The qth moments for the $\mathit{IG}$ process are found in the paper [11]:
\[ E({G_{t}^{q}})=\sqrt{\frac{2}{\pi }\delta }{\left(\frac{\delta }{\gamma }\right)^{q-1/2}}{t^{q+1/2}}{e^{\delta \gamma t}}{K_{q-1/2}}(\delta \gamma t),\]
where ${K_{q}}(\omega )$ is the modified Bessel function of the third kind with index q. Moreover, in [11] it was shown that if $t\to \infty $, then
For the standard distribution, we have $E(G(t))\sim t$, $\mathit{var}(G(t))\sim {t^{2}}$.The algorithm of the simulation of the $\mathit{IG}$ process ${G_{t}}$ for time points ${t_{1}}=\frac{1}{n},{t_{2}}=\frac{2}{n},\dots ,{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_{dt}}\sim \varrho (dt,1)$ for $i=1,2,\dots n$ and $dt=\frac{1}{n}$, so following [26] we generate n i.i.d. $\mathit{IG}$ variables ${F_{i}}$ assuming $\gamma =\delta =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=dt+\frac{X}{2}-\frac{1}{2}\sqrt{4dtX+{X^{2}}}$. According this algorithm, if $U\le \frac{dt}{dt+Y}$ then return Y; otherwise return $\frac{{(dt)^{2}}}{Y}$. After assigning ${G_{{t_{0}}}}=0$ and ${G_{{t_{i}}}}={\textstyle\sum _{j=1}^{i}}{F_{j}}$, $i=1,2,\dots ,n$, we obtain ${G_{{t_{1}}}},{G_{{t_{2}}}},\dots ,{G_{{t_{n}}}}$, which are simulated values of the $\mathit{IG}$ process at times ${t_{1}},{t_{2}},\dots ,{t_{n}}$, respectively. The simulated trajectories of the standard $\mathit{IG}$ process are shown in Figure 1.
Now we discuss the properties of the inverse subordinator (hitting time) ${H_{t}}$ defined by (4). The inverse to the inverse Gaussian ($\mathit{IIG}$) process is not Lévy and has monotonically increasing continuous sample paths. Moreover, the sample paths of the $\mathit{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]):
It is worth noticing that if $\gamma =0$ then Laplace exponent (8) has the form ${\Psi _{\gamma =0}}(u)=(\delta \sqrt{2}){u^{1/2}}$, and it means that we have the stable distribution $\mathit{TS}(1/2,0,\delta )$. In the case $\gamma \ne 0$ we have the tempered stable distribution $\mathit{TS}(1/2,{\gamma ^{2}},\delta )$. For the particular case $\gamma =0$, $\delta =1$ we have $1/2-$stable distribution $S(\alpha )$, $\alpha =1/2$, and the density function for the inverse $\mathit{IG}$ process is equal
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.
(10)
\[ h(x,t)=\frac{\delta }{\pi }{e^{\delta \gamma x-\frac{{\gamma ^{2}}}{2}}}{\int _{0}^{\infty }}\frac{{e^{-ty}}}{y+\frac{{\gamma ^{2}}}{2}}(\gamma \sin (\delta x\sqrt{2y})+\sqrt{2y}\cos (\delta x\sqrt{2y}))dy.\]The qth moments for the $\mathit{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 $\mathit{Var}({H_{t}})$. If $t\to \infty $, then
For the standard distribution we have $E({H_{t}})\sim t$, $var({H_{t}})\sim {t^{2}}$.
(11)
\[ E({H_{t}})\sim \left\{\begin{array}{c@{\hskip10.0pt}c}\left(\frac{\gamma }{\delta }\right)t,& \gamma \gt 0,\\ {} \left(\frac{1}{\delta }\sqrt{\frac{2t}{\pi }}\right)t,& \gamma =0,\end{array}\right.\]To simulate the approximate trajectory of the inverse subordinator ${H_{t}}$, we define ${H_{t}^{\Delta }}$ as follows [26]:
where Δ is the step length and ${G_{\Delta n}}$ is the value of the Inverse Gaussian process ${G_{t}}$ at $\Delta n$. The simulation of the trajectory ${H_{t}}$ for $\gamma =\delta =1$ is demonstrated in Figure 2.
(13)
\[ {H_{t}^{\Delta }}=[\min \{n\in N:{G_{\Delta n}}\gt t\}-1]\Delta ,\hspace{1em}n=1,2,\dots ,\]For modeling of the stochastic subdiffusive GBM we propose the iterative scheme
where ε is a white noise with the normal standard distribution, and ${H_{t}^{\Delta }}$ are given in (13). The trajectories of the subdiffusion GBM with the inverse to the $\mathit{IG}$ process are demonstrated in Figure 3.
(14)
\[ {S_{t+1}}={S_{t}}+\mu {S_{t}}{H_{t}^{\Delta }}+\sigma {S_{t}}\sqrt{{H_{t}^{\Delta }}}{\varepsilon _{t}},\]2.2 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
where ${S_{{H_{t}}}}$ is a stochastic process given by (3).
Then denote ${\tau _{t}}={H_{t}}-{H_{t-1}}$ and obtain the stock returns in the form
using the scaling law of the Brownian motion.
Proposition 1.
Let ${X_{t}}$ be a stochastic process given by (15), then for any integer $k\ge 0$:
-
2. There is a long-range dependence in the squared returns:in particular, ${\mathrm{Cov}_{\mu =0}}({X_{t}^{2}},{X_{t+k}^{2}})={\sigma ^{4}}\mathrm{Cov}({\tau _{t}},{\tau _{t+k}})$.
(17)
\[\begin{aligned}{}& \mathrm{Cov}({X_{t}^{2}},{X_{t+k}^{2}})\\ {} & \hspace{1em}={\sigma ^{4}}\mathrm{Cov}({\tau _{t}},{\tau _{t+k}})+{\mu ^{4}}\mathrm{Cov}({\tau _{t}^{2}},{\tau _{t+k}^{2}})+2{\mu ^{2}}{\sigma ^{2}}\mathrm{Cov}({\tau _{t}^{2}},{\tau _{t+k}});\end{aligned}\]
Proof.
For the covariance function of the log returns process ${X_{t}}$ (15), we obtain
\[\begin{aligned}{}& \mathrm{Cov}({X_{t}},{X_{t+k}})=\mathrm{Cov}(\mu {\tau _{t}}+\sigma \sqrt{{\tau _{t}}}{B_{1}^{(1)}},\mu {\tau _{t+k}}\sigma \sqrt{{\tau _{t+k}}}{B_{1}^{(2)}})\\ {} & \hspace{1em}=\mathrm{E}\left[\left(\mu {\tau _{t}}+\sigma \sqrt{{\tau _{t}}}{B_{1}^{(1)}}-\mu E[{\tau _{t}}]\right)\left(\mu {\tau _{t+k}}+\sigma \sqrt{{\tau _{t+k}}}{B_{1}^{(2)}}-\mu E[{\tau _{t+k}}]\right)\right]\\ {} & \hspace{1em}={\mu ^{2}}\left(\mathrm{E}[{\tau _{t}}{\tau _{t+k}}]-\mathrm{E}[{\tau _{t}}]\mathrm{E}[{\tau _{t+k}}]\right)={\mu ^{2}}\mathrm{Cov}({\tau _{t}},{\tau _{t+k}}),\end{aligned}\]
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
\[\begin{aligned}{}& \mathrm{Cov}({X_{t}^{2}},{X_{t+k}^{2}})=\mathrm{Cov}\left({(\mu {\tau _{t}}+\sigma \sqrt{{\tau _{t}}}{B_{1}^{(1)}})^{2}},{(\mu {\tau _{t+k}}+\sigma \sqrt{{\tau _{t+k}}}{B_{1}^{(2)}})^{2}}\right)\\ {} & \hspace{1em}={\sigma ^{4}}\mathrm{Cov}({\tau _{t}},{\tau _{t+k}})+{\mu ^{4}}\mathrm{Cov}({\tau _{t}^{2}},{\tau _{t+k}^{2}})\\ {} & \hspace{2em}+{\mu ^{2}}{\sigma ^{2}}\mathrm{Cov}({\tau _{t}^{2}},{\tau _{t+k}})+{\mu ^{2}}{\sigma ^{2}}\mathrm{Cov}({\tau _{t+k}},{\tau _{t}^{2}}),\end{aligned}\]
where from the last expression we immediately get the result. □Thus, if the increments ${\tau _{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.
2.3 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:
where ${\Phi _{t}}$ is the integro-differential operator defined as
with the memory kernel $M(t)$ defined via its Laplace transform
where $\Psi (u)$ is the Laplace exponent of the $\mathit{IG}$ subordinator and can be written as in (8). This implies that the memory kernel $M(t)$ can be expressed as
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).
(18)
\[ \frac{\partial w}{\partial t}={\Phi _{t}}\left[-\mu \frac{\partial }{\partial x}+\frac{{\sigma ^{2}}}{2}\frac{{\partial ^{2}}}{\partial {x^{2}}}\right]w(x,t),\]Thus, the possibility of numerical computing of the probability density function $w(t)$ for the studied subdiffusive stock process (with $\mathit{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 $(\Omega ,\mathcal{F},P)$ be the probability space. The value-at-risk of level α, $0\lt \alpha \le 1$, is a probability functional defined as α-quantile of the profit (loss) function $Y\in L(\Omega )$:
where W is the distribution function of Y, ${W^{-1}}$ is the quantile function of α, $0\lt \alpha \le 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
and the value-at-risk of level α, $0\lt \alpha \le 1$, for the random variable Y is
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
3 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
where
both are functions of five parameters $T,K,{S_{0}},r,\sigma $, and $\Phi (\cdot )$ is the standard normal cumulative distribution function.
(24)
\[ {d_{1}}=\frac{\log \frac{{S_{0}}}{K}+rT+\frac{1}{2}{\sigma ^{2}}T}{\sigma \sqrt{T}},\hspace{1em}{d_{2}}=\frac{\log \frac{{S_{0}}}{K}+rT-\frac{1}{2}{\sigma ^{2}}T}{\sigma \sqrt{T}}\]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 $\mathit{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 $(\Omega ,\mathcal{F},P)$. Here, Ω is the sample space, $\mathcal{F}$ contains all statements that can be made about the behavior of prices, and P is the “objective” probability measure. We denote by $({\mathcal{F}_{t}}),t\in [0,T]$, the information about the history of asset prices ${S_{H}}$ up to time t. $({\mathcal{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.
Proposition 2.
The European call option price ${C_{H}}={C_{H}}\left(S,K,T,\sigma \right)$ for a time-changed version of the B-S model with $\mathit{IG}$ subordinator satisfies
where
and $C(T)=C(S,K,T,\sigma )$ is given by (23).
(25)
\[ {C_{H}}={\int _{0}^{\infty }}\hspace{-0.1667em}\hspace{-0.1667em}{\int _{0}^{\infty }}\hspace{-0.1667em}\hspace{-0.1667em}C(x)\frac{\delta }{\pi }{e^{\delta \gamma x-\frac{{\gamma ^{2}}}{2}}}\frac{{e^{-Ty}}}{y+\frac{{\gamma ^{2}}}{2}}\left(\gamma \sin (\rho (x,y))+\sqrt{2y}\cos (\rho (x,y))\right)dydx,\]Proof.
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
was introduced, where $\tau =\frac{\mu +{\sigma ^{2}}}{\sigma }$ and $A\in \mathcal{F}$, and it was proved that subdiffusive GBM (3) with α-stable subordinator is a martingale with respect to $\mathbb{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.
(26)
\[ \mathbb{Q}(A)={\int _{A}}\exp \left\{-\tau {B_{{H_{T}}}}-\frac{{\tau ^{2}}}{2}{H_{T}}\right\}d\mathbb{P}\]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 $\mathbb{Q}$ is not unique, and therefore the subdiffusive markets are generally incomplete. Particularly it applies to inverse-tempered stable subordinators. Thus, the subdiffusive $\mathit{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 $\mathit{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 $\mathbb{Q}$ (26):
Here, ${h_{\Psi }}(x,T)$ is the PDF of ${H_{T}}$ for the subordinator with Lévy exponent Ψ and $C(S,K,T,\sigma )$ is given by (23).
\[ {C_{H}}\left(S,K,T,\sigma \right)=\left\langle C\left(S,K,{H_{T}},\sigma \right)\right\rangle ={E^{\mathbb{Q}}}\left({e^{-r{H_{T}}}}{({S_{{H_{T}}}}-K)^{+}}|{\mathcal{F}_{0}}\right).\]
Therefore, conditioning on ${H_{T}}$, we obtain
(27)
\[ {C_{H}}\left(S,K,T,\sigma \right)={\int _{0}^{\infty }}C\left(S,K,x,\sigma \right){h_{\Psi }}(x,T)dx.\]Remark 1.
The presented option pricing subdiffusive model with $IG(\delta ,\gamma )$ 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 $\mathit{IG}$ subordinator with parameters δ and γ can be considered as tempered stable Lévy subordinator $\mathit{TS}(1/2,{\gamma ^{2}},\delta )$, 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 ${\gamma ^{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 $\gamma =0$ we get a regime where prices exhibit maximum staleness at all temporal scales. Therefore, $\mathit{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).
Remark 2.
In the above equation (27) we can evaluate the subdiffusive call price $C(\cdot )$ 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)
where $C(S,K,T,\sigma )$ is taken from the Black–Scholes option pricing formula (23).
(28)
\[ {C_{H}}(S,K,T,r,\sigma )=\langle C(S,K,{H_{T}},\sigma )\rangle =\frac{1}{n}{\sum \limits_{i=1}^{n}}C(S,K,{H_{T}^{(i)}},\sigma ),\]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 $\Psi (.)$, which is often called the Bernstein function.
The next proposition is the application of this PIDE to B-S subdiffusion with $\mathit{IG}$ subordinators.
Proposition 3.
The European call option price ${C_{H}}(T,k)={C_{H}}\left(S,K,T,\sigma \right)$ for a time-changed version of the B-S model with $\mathit{IG}$ subordinator is a solution of a fractional PIDE equation
\[\begin{aligned}{}& {\int _{0}^{T}}\delta \frac{\partial }{\partial t}{C_{H}}(T-s,k)\left(\gamma (\Phi \left(\gamma \sqrt{s}\right)-1)+\frac{{e^{-\frac{s{\gamma ^{2}}}{2}}}}{\sqrt{2\pi s}}\right)ds\\ {} & \hspace{1em}=-\frac{r}{2}\frac{\partial }{\partial k}{C_{H}}(T,k)+\frac{{\sigma ^{2}}}{4}\frac{{\partial ^{2}}}{\partial {k^{2}}}{C_{H}}(T,k),\end{aligned}\]
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}}}})^{+}}$.
Proof.
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
where ${_{}^{\Psi }}Du(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])
The function Ψ is the Lévy exponent for a given subordinator ${G_{t}}$, and we use its Lévy–Khintchine representation
where $\tilde{\nu }$ is the Lévy measure for this subordinator.
(29)
\[ {_{}^{\Psi }}{DC_{H}}(T,k)=-r\frac{\partial }{\partial k}{C_{H}}(T,k)+\frac{{\sigma ^{2}}}{2}\frac{{\partial ^{2}}}{\partial {k^{2}}}{C_{H}}(T,k),\](30)
\[ {_{}^{\Psi }}Du(t)=b\frac{d}{dt}u(t)+{\int _{0}^{t}}\frac{\partial }{\partial t}u(t-s)\nu (s)ds.\]For $\mathit{IG}$ Lévy subordinator the Lévy measure $\tilde{\nu }(s)$ is defined by (7). The integral kernel $\nu (s)$ is the integral of $\tilde{\nu }$ over $(s,\infty )$:
Thus, if we substitute (31) in (29), the proposition follows. □
\[ \nu (s)={\int _{s}^{+\infty }}\hspace{-0.1667em}\hspace{-0.1667em}\frac{\delta }{\sqrt{2\pi {x^{3}}}}{e^{\left(-\frac{{\gamma ^{2}}x}{2}\right)}}dx=\frac{\delta \gamma }{2\sqrt{\pi }}{\int _{\frac{s{\gamma ^{2}}}{2}}^{+\infty }}\frac{1}{\sqrt{{z^{3}}}}{e^{-z}}dz=\frac{\delta \gamma }{2\sqrt{\pi }}\Gamma (-1/2,\frac{s{\gamma ^{2}}}{2}),\]
where $\Gamma (a,b)$ is the upper incomplete Gamma function. If we use the fact that $\Gamma (a+1,b)=a\Gamma (a,b)+{b^{a}}{e^{-b}}$ and expression for $\Gamma (1/2,b)$, then D-C derivative (30) for $\mathit{IG}$ subordinator can be written as
\[ {_{}^{\Psi }}Du(t)=\delta \gamma {\int _{0}^{t}}\frac{\partial }{\partial t}u(t-s)\left(\mathit{erf}\left(\frac{\gamma \sqrt{2s}}{2}\right)-1+\frac{2{e^{-\frac{s{\gamma ^{2}}}{2}}}}{\gamma \sqrt{2\pi s}}\right)ds,\]
where $\mathit{erf}(.)$ is the error function, which can be expressed in terms of the standard normal cumulative distribution function $\Phi (\cdot )$:
(31)
\[ {_{}^{\Psi }}Du(t)=2\delta {\int _{0}^{t}}\frac{\partial }{\partial t}u(t-s)\left(\gamma (\Phi \left(\gamma s\right)-1)+\frac{{e^{-\frac{s{\gamma ^{2}}}{2}}}}{\sqrt{2\pi s}}\right)ds.\]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 $\mathit{IG}$ subordinator.
4 Numerical illustration
We illustrate the call option pricing for the subdiffusive $\mathit{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$, $\mu =-0.0015$.
To simplify the calculations, we put parameters $\gamma =\delta =1$ because this assumption allows getting the desired result regarding the asymptotic behavior of the hitting time. In this case, when $t\to \infty $, it follows that $E({H_{t}})\sim t$, $\mathit{var}({H_{t}})\sim {t^{2}}$ (11)–(12). So, we estimate only the σ parameters based on the least squares technique. We obtain that the value $\sigma =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.
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.
Table 1.
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):
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.
Table 2.
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 |
Table 3.
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.
5 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 $\mathit{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.
