Pricing, Risk and Volatility in Subordinated Market Models

Abstract

We consider several market models, where time is subordinated to a stochastic process. These models are based on various time changes in the Lévy processes driving asset returns, or on fractional extensions of the diffusion equation; they were introduced to capture complex phenomena such as volatility clustering or long memory. After recalling recent results on option pricing in subordinated market models, we establish several analytical formulas for market sensitivities and portfolio performance in this class of models, and discuss some useful approximations when options are not far from the money. We also provide some tools for volatility modelling and delta hedging, as well as comparisons with numerical Fourier techniques.

1. Introduction

In this opening section, we provide a general introduction to the class of subordinated market models; we also present the key points investigated in the paper, as well as the work’s overall structure.

1.1. Time Subordination in Financial Modelling

Among the most striking patterns that are observable in financial time series are the phenomena of regime switching, clustering, and long memory or autocorrelation (see e.g., Cont (2007) and references therein). Such stylized facts have been evidenced for several decades, Mandelbrot famously remarking in Mandelbrot (1963) that large price changes tend to cluster together (“large changes tend to be followed by large changes, of either sign, and small changes tend to be followed by small changes”), thus creating periods of market turbulence (high volatility) alternating with periods of relative calm (low volatility). These empirical observations can be described, among other approaches, by agent based models focusing on economic interpretation, such as Lux and Marchesi (2000); Niu and Wang (2013), by tools from statistical mechanics and econophysics (Krawiecki et al. 2002), or by the introduction of multifractals (Calvet and Fischer 2008).

Another prominent approach to describe this subtle volatility behavior consists of introducing a time change in the stochastic process driving the market prices. Besides stochastic volatility, time changed market models also capture several stylized facts, like non-Normality of returns (the presence of jumps, asymmetry) and negative correlation between the returns and their volatility (see a complete overview in Carr and Wu (2004)). They are motivated by the observation that market participants do not operate uniformly through a trading period, but, on the contrary, the volume, and frequency of transactions greatly vary over time. Following the terminology of Geman (2009), the time process is called the stochastic clock, or business time, while the stochastic process for the underlying market (a Brownian motion, or a more general Lévy process) is said to evolve in operational time.

Historically, the first introduction of a time change in a diffusion process goes back to Bochner (1949); it was first applied to financial modeling in Clark (1973) in the context of the cotton futures market and for a continuous-time change. During the late 1990s and early 2000s, the approach was extended to discontinuous time changes, with the introduction of subordinators (i.e., non-negative Lévy processes, see the theoretical details in Bertoin (1999)). In other words, the business time now admits increasing staircase-like realizations, describing peak periods of activity (following, for instance, earning announcements, central bank reports, or major political events) alternating with less busy periods. Perhaps the best-known subordinators are the Gamma process, like in the Variance Gamma (VG) model by Madan et al. (1998), and the inverse Gaussian process, like in the Normal inverse Gaussian (NIG) model by Barndorff-Nielsen (1997). Let us also mention that subordination has been successfully applied in many other fields of applied science. For instance, Gamma subordination has been employed for modeling the deterioration of production equipment in order to optimize their maintenance (see de Jonge et al. (2017) and references therein), and inverse Gaussian subordination was originally introduced in Barndorff-Nielsen (1977) to model the influence of wind on dunes and beach sands.

Recently, a new type of time subordination, based on fractional calculus, has emerged. Indeed, Lévy processes are closely related to fractional calculus because, for many of them (including stable and tempered stable processes), their probability densities satisfy a space fractional diffusion equation (see details and applications to option pricing in Cartea and del-Castillo-Negrete (2007) and in Luchko et al. (2019)). By also allowing the time derivative to be fractional, as, e.g., in Jizba et al. (2018); Kleinert and Korbel (2016); Korbel and Luchko (2016); Tarasov (2019); Tomovski et al. (2020), it provides a new type of subordinated models: while the order of the space fractional derivative controls the heavy tail behavior of the distribution of returns, the order of the time fractional derivative acts as a temporal subordination parameter whose purpose is to capture time-related phenomena, such as temporal risk redistribution. This model, which we shall refer to as the fractional diffusion (FD) model, is an alternative to time-change models, or to subordinated random walks (Gorenflo et al. 2006).

Regarding the practical implementation and valuation of financial derivatives within subordinated market models, the literature is dominated by numerical techniques. In time changed models notably, tools from Fourier transform (Lewis 2001) or Fast Fourier transform (Carr and Madan 1999), and their many refinements, such as the COS method by Fang and Osterlee (2008) or the PROJ method by Kirkby (2015). These methods are popular, notably because such models’ characteristic functions are known in relatively simple closed-forms. Similarly, Cui et al. (2019) provides a numerical pricing framework for a general time changed Markov processes, and Li and Linetsky (2014) employs eigenfunction expansion techniques. However, recently, closed-form pricing formulas have been derived, for the VG model in Aguilar (2020a) and for the NIG model in Aguilar (2020b). The technique has also been employed in the FD model, for vanilla payoffs in Aguilar et al. (2018) and for more exotic options in Aguilar (2020c).

In this paper, we extend these pricing tools to the calculation of risk sensitivities and to profit-and-loss (P&L) explanation, and we provide comparisons between time changed models (such as the NIG and the VG models) and the FD model. Like for the pricing case, risk sensitivities in the context of time changed market models (and of Lévy market models in general) are traditionally evaluated by means of numerical methods based on Fourier inversion (Eberlein et al. 2010; Takahashi and Yamazaki 2008); in the present paper, we will therefore show that they can be expressed in a tractable way, under the form of fast convergent series whose terms explicitly depend on the model parameters. This will allow for us to construct and compare the performance of option based portfolios, and discuss, both quantitatively and qualitatively, the impact on the parameters on risks and P&L. Related topics, such as volatility modeling, will also be discussed.

1.2. Contributions of the Paper

Our purpose in the present work is to investigate and provide details on the following key points:

(a)

demonstrate that the recent pricing formulas for the VG, NIG and FD models are precise and fast converging, and can be successfully used for other applications (e.g., calculations of volatility curve);

(b)

provides efficient closed-form formulas for first and second-order risk sensitivities (Delta, Gamma) and compare them with numerical techniques; and,

(c)

deduce from these formulas several practical features regarding delta-hedging policies and portfolio performance.

1.3. Structure of the Paper

The paper is organized as follows: in Section 2 we recall some fundamental concepts on Lévy processes and option pricing and, in Section 3, we introduce the class of subordinated market models and their main implications in financial modeling.Subsequently, in Section 4, we mention the various closed pricing formulas that have been obtained for this class of models. Approximating these formulas when options are not far from the money, we establish formulas for computing the market volatility in this configuration, thus generalizing the usual Black-Scholes implied volatility. In Section 5 (resp. Section 6), we derive the expressions for the first (resp. second) order market sensitivities, and for the P&L of a delta-hedged portfolio; the impact of the subordination parameter is discussed, and a comparison with numerical techniques is provided. Last, Section 7 is dedicated to concluding remarks.

2. Exponential Lévy Processes

Let us start by recalling some fundamentals on Lévy processes (see full details in Sato (1999) and in Cont and Tankov (2004) for their applications to financial modeling) and, following the classical setup of Schoutens (2003), how they are implemented for the purpose of option pricing.

2.1. Basics of Lévy Processes

Let $(\Omega ,\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0},\mathbb{P})$ be a probability space that is equipped with its natural filtration. Recall that a Lévy process ${\left\{X_t\right\}}_{t\ge 0}$ is a stochastically continuous process satisfying $X_0=0$ ($\mathbb{P}$-almost surely), and whose increments are independent and stationary. This implies that the characteristic function $\Psi (u,t):={\mathbb{E}}^{\mathbb{P}}\left[e^{iu{X}_t}\right]$ of a Lévy process has a semi-group structure and it admits an infinitesimal generator $\psi \left(u\right)$, called Lévy symbol or characteristic exponent, which satisfies

$$\Psi (u,t)=e^{t\psi \left(u\right)},\psi \left(u\right):=log\Psi (u,1).$$

(1)

The characteristic exponent is entirely determined by the triplet $(a,b,\Pi (\mathrm{d}x\left)\right)$, which corresponds to Lévy–Khintchine representation

$$\psi \left(u\right)=aiu-\frac{1}{2}{b}^2{u}^2+\underset{-\infty }{\overset{+\infty }{\int }}\left(e^{iux}-1-iux{\mathrm{𝟙}}_{\left\{\right|x|<1\}}\right)\Pi \left(\mathrm{d}x\right),$$

(2)

where a is the drift and is b the Brownian (or diffusion) component. The measure $\Pi \left(\mathrm{d}x\right)$, assumed to be concentrated on $\mathbb{R}\backslash \left\{0\right\}$ and satisfy

$$\underset{-\infty }{\overset{+\infty }{\int }}\min (1,x^2)\Pi \left(\mathrm{d}x\right)<\infty ,$$

(3)

is called the Lévy measure of the process, and determines its tail behavior and the distribution of jumps. When $\Pi \left(\mathbb{R}\right)<\infty$, one speaks of a process with finite activity (or intensity); this is the case for jump-diffusion processes, like in the Kou Model (Kou 2002) or the Merton model (Merton 1976), where only a finite number of jumps can occur on each time interval. When $\Pi \left(\mathbb{R}\right)=\infty$, one speaks of a process with infinite activity (or intensity); this class is far richer, because an infinite number of jumps can occur on any finite time interval and, as a consequence, no Brownian component b is even needed to generate a very complex dynamics. A prominent model with infinite activity is the Variance Gamma process, introduced in Madan et al. (1998). When $\Pi \left({\mathbb{R}}_{-}\right)=0$ (i.e., the process has only positive jumps), one speaks of a subordinator.

An important class of Lévy measures, which will be of particular interest to us in this paper, corresponds to the so-called class of tempered stable processes:

$$\Pi \left(\mathrm{d}x\right):=\left[\frac{c_{+}{e}^{-{\lambda }_{+}x}}{x^{1+{\alpha }_{+}}}{\mathrm{𝟙}}_{\{x>0\}}+\frac{c_{-}{e}^{-{\lambda }_{-}\left|x\right|}}{{\left|x\right|}^{1+{\alpha }_{-}}}{\mathrm{𝟙}}_{\{x<0\}}\right]\mathrm{d}x.$$

(4)

This class contains several sub-classes, such as the tempered stable subordinators ($c_{-}=0$) or the stable processes (${\lambda }_{+}={\lambda }_{-}=0$). When $c_{+}=c_{-}:=C$, ${\alpha }_{+}={\alpha }_{-}:=Y$, ${\lambda }_{-}:=G$ and ${\lambda }_{+}:=M$, one speaks of a CGMY process (introduced in Carr et al. (2002)). By requiring the further restriction that $Y=0$, we obtain the Variance Gamma process of Madan et al. (1998); the symmetric case $G=M$ was considered earlier in Madan and Seneta (1990). We also note that the CGMY (and VG) models are members of the KoBoL family, see Boyarchenko and Levendorskiĭ (2000).

2.2. Exponential Lévy Motions

Let $T>0$ and $S:t\in [0,T]\to S_t$ be the market price of some financial asset, seen as the realization of a time dependent random variable ${\left\{S_t\right\}}_{t\in [0,T]}$ on the canonical space $\Omega ={\mathbb{R}}_{+}$. We assume that there exists a risk-neutral measure $\mathbb{Q}$ under which the instantaneous variations of $S_t$ can be written down as:

$$\frac{\mathrm{d}{S}_t}{S_t}=(r-q)\mathrm{d}t+\mathrm{d}{X}_t$$

(5)

where $r\in \mathbb{R}$ is the risk-free interest rate and $q\ge 0$ is the dividend yield (both being assumed to be deterministic and continuously compounded), and where ${\left\{X_t\right\}}_{t\in [0,T]}$ is a Lévy process. Under the dynamics (5), the terminal market price is given by

$$S_T=S_t{e}^{(r-q+\omega )\tau +X_{\tau }},$$

(6)

where $\tau :=T-t$ is the time horizon and $\omega$ is the martingale adjustment (also called convexity adjustment, or compensator), which is determined by the martingale condition

$${\mathbb{E}}^{\mathbb{Q}}\left[S_T\right|{\mathcal{F}}_t]=e^{(r-q)\tau }{S}_t.$$

(7)

Given the form of the exponential process (6), the condition (7) is equivalent to:

$$\omega =-\psi (-i)=-log{\mathbb{E}}^{\mathbb{P}}\left[e^{X_1}\right].$$

(8)

2.3. Option Pricing

Given a path-independent payoff function $\mathcal{P}$, i.e., a positive function depending only on the terminal value $S_T$ of the market price and on some strike parameters $K_1,\dots ,K_N>0$, then the value at time t of a contingent claim delivering a payoff $\mathcal{P}$ at its maturity is equal to the following risk-neutral expectation:

$$\mathcal{C}={\mathbb{E}}^{\mathbb{Q}}\left[e^{-r\tau }\mathcal{P}(S_T,K_1,\dots ,K_n)|{\mathcal{F}}_t\right].$$

(9)

If the Lévy process ${\left\{X_t\right\}}_{t\in [0,T]}$ admits a density $f(x,t)$, then the conditional expectation (9) can be achieved by integrating all possible realizations for the payoff over the Lévy density, thus resulting in:

$$\mathcal{C}=e^{-r\tau }\underset{-\infty }{\overset{+\infty }{\int }}\mathcal{P}(S_t{e}^{(r-q+\omega )\tau +x},K_1,\dots ,K_n)f(x,\tau )\mathrm{d}x.$$

(10)

3. Subordinated Models

In this section, we introduce the class of subordinated market models, which is, models for which time is driven by a particular subordinator. We also provide a review of their main financial applications.

3.1. Exponential VG Model

3.1.1. Model Characteristics

In the exponential VG model by (Madan et al. 1998), one chooses the Lévy process in (5) to be a VG process; this process is defined by

$$X_t^{\left(VG\right)}:=\theta G_t+\sigma W_{G_t}$$

(11)

where $W_t$ is the standard Wiener process, and $\gamma (t,1,\nu )$ is a Gamma process (i.e., a process whose increments $\gamma (t+h,1,\nu )-\gamma (t,1,\nu )$ follow a Gamma distribution with mean $1\times h$ and variance $\nu \times h$). It follows from definition (11) that the VG process is actually distributed according to a so-called Normal variance–mean mixture (see Barndorff-Nielsen et al. (1982)), where the mixing distribution is the Gamma distribution; this distribution materializes the business time, and it is a particular case of a tempered stable subordinator, as it admits the following Lévy measure (see Sato (1999) for instance):

$${\Pi }_G\left(\mathrm{d}x\right)=\frac{1}{\nu }\frac{e^{-\frac{1}{\nu }x}}{x}{\mathrm{𝟙}}_{\{x>0\}}\mathrm{d}x.$$

(12)

It follows from (12) that ${\Pi }_G\left(\mathbb{R}\right)=\infty$, which means that the Gamma process has infinite activity; note also that the Gamma measure (12) is concentrated around 0, which means that most jumps in the business time are small, and become bigger in the high $\nu$ regime. The VG process is actually a tempered stable process itself (and, more precisely, a CGMY process), its Lévy measure admitting the following representation:

$${\Pi }_{VG}\left(\mathrm{d}x\right)=\frac{e^{\frac{\theta x}{{\sigma }^2}}}{\nu \left|x\right|}{e}^{-\frac{\sqrt{\frac{{\theta }^2}{{\sigma }^2}+\frac{2}{\nu }}}{\sigma }\left|x\right|}\mathrm{d}x.$$

(13)

Note that (13) is symmetric around the origin when $\theta =0$ (i.e., positive and negative jumps in asset prices occur with the same probability). The density function of the VG process is obtained by integrating the normal density conditionally to the Gamma time-change, and it reads:

$$f_{VG}(x,t)=\frac{2e^{\theta \frac{x}{{\sigma }^2}}}{{\nu }^{\frac{t}{\nu }}\sqrt{2\pi }\sigma \Gamma \left(\frac{t}{\nu }\right)}{\left(\frac{x^2}{\frac{2{\sigma }^2}{\nu }+{\theta }^2}\right)}^{\frac{t}{2\nu }-\frac{1}{4}}{K}_{\frac{t}{\nu }-\frac{1}{2}}\left(\frac{1}{{\sigma }^2}\sqrt{\frac{2{\sigma }^2}{\nu }+{\theta }^2}\left|x\right|\right)$$

(14)

where $K_a\left(X\right)$ denotes the modified Bessel function of the second kind, sometimes also called MacDonald function (see definition and properties in Abramowitz and Stegun (1972)). The Lévy symbol is known in the exact form:

$${\psi }_{VG}\left(u\right)=-\frac{1}{\nu }log\left(1-i\theta \nu u+\frac{{\sigma }^2\nu }{2}{u}^2\right),$$

(15)

allowing for a simple representation for the VG martingale adjustment:

$${\omega }_{VG}=-{\psi }_{VG}(-i)=\frac{1}{\nu }log\left(1-\theta \nu -\frac{{\sigma }^2\nu }{2}\right).$$

(16)

Remark 1.

Note that, when $\theta =0$ and $\nu \to \infty$, then ${\omega }_{VG}\to -\frac{{\sigma }^2}{2}$, which is the usual Gaussian adjustment, and, in this limit, the exponential VG model degenerates into the Black–Scholes model (Black and Scholes 1973). The limiting regime, VG$(\sigma ,\nu ,0)\stackrel{\nu \to 0}{⟶}$ BS$\left(\sigma \right)$, is illustrated in Figure 1 for decreasing ν. In particular, ν directly controls the excess kurtosis for the VG model.

3.1.2. Financial Applications

As already noted, the presence of the subordination parameter $\nu$ is particularly attractive for modeling time-induced phenomena, as it allows for a non-uniform passage of time. When $\nu$ is small, realizations of the Gamma subordinator are quasi-linear, which corresponds to a situation where the business time is the same as the operational time. On the contrary, for bigger values of $\nu$, realizations of the Gamma process are highly discontinuous and staircase-like (because the process is non decreasing), capturing the alternation of intense and quieter trading periods.

The exponential VG model has been successfully tested on real market data and shown to perform better than Black–Scholes or Jump-Diffusion models in multiple situations, e.g., for European-style options on the HSI index in Lam et al. (2002) or for currency options in Madan and Dual (2005). Several extensions of the model have been subsequently developed, such as the generalization of the subordination to a bivariate or multivariate Brownian motion (Luciano and Schoutens 2006; Semeraro 2008) (with an application to basket options calibration in Linders and Stassen (2015)). Other recent extensions include the possibility of negative jumps in the linear drift rate of the price process in Ivanov (2018). Last, let us also mention that the exponential VG model has also found its way to applications in other fields of quantitative finance, such as credit risk in Fiorani et al. (2007).

3.2. Exponential NIG Model

3.2.1. Model Characteristics

In the exponential NIG model (see Barndorff-Nielsen (1997)), one chooses the Lévy process in (5) to be the NIG process, defined by

$$X_t^{\left(NIG\right)}=\beta {\delta }^2{I}_t+\delta W_{I_t}$$

(17)

where ${\left\{I_t\right\}}_{t\in [0,T]}$ follows an Inverse Gamma distribution of shape $\delta \sqrt{{\alpha }^2-{\beta }^2}$ and mean rate 1. $\alpha >0$ is a tail or steepness parameter controlling the kurtosis of the NIG distribution; the large $\alpha$ regime gives birth to light tails, while small $\alpha$ corresponds to heavier tails. $\beta \in (-\alpha ,\alpha -1)$ is the skewness parameter: $\beta <0$ (resp. $\beta >0$) implies that the distribution is skewed to the left (resp. the right), and $\beta =0$ that the distribution is symmetric. $\delta >0$ is the scale parameter and it plays an analogous role to the variance term ${\sigma }^2$ in the Normal distribution. Let us mention that a location parameter $\mu \in \mathbb{R}$ can also be incorporated, but it has no impact on option prices (see e.g., Aguilar (2020b)), and, therefore, we will assume that it is equal to 0. Let us also note that, again, we are in the presence of a tempered stable subordination, as the Lévy measure of the Inverse Gamma process ${\left\{I_t\right\}}_{t\in [0,T]}$ satisfies

$${\Pi }_{IG}\left(\mathrm{d}x\right)=\frac{e^{-\frac{{\delta }^2({\alpha }^2-{\beta }^2)}{2}x}}{x^{\frac{3}{2}}}{\mathrm{𝟙}}_{\{x>0\}}\mathrm{d}x,$$

(18)

while the Lévy measure of the NIG process itself is given by

$${\Pi }_{NIG}\left(\mathrm{d}x\right):=\frac{\alpha \delta }{\pi }{e}^{\beta x}\frac{K_1\left(\alpha \right|x\left|\right)}{\left|x\right|}\mathrm{d}x.$$

(19)

It follows from definition (17) that, like in the VG case, the NIG process is also distributed according to a Normal variance-mean mixture, where the mixing distribution is now the IG distribution; this mixture is a particular case of the more general class of hyperbolic processes (see discussion and applications to finance in Eberlein and Keller (1995)), the mixing distributions in that case being the Generalized Inverse Gaussian (GIG) distribution. The probability density function for the NIG process is obtained by an integration of the Normal density over the IG distribution and it reads

$$f_{NIG}(x,t):=\frac{\alpha \delta t}{\pi }{e}^{\delta t\sqrt{{\alpha }^2-{\beta }^2}+\beta (x-\mu t)}\frac{K_1\left(\alpha \sqrt{{\left(\delta t\right)}^2+{(x-\mu t)}^2}\right)}{\sqrt{{\left(\delta t\right)}^2+{(x-\mu t)}^2}},$$

(20)

and its Lévy symbol is given by

$${\psi }_{NIG}\left(u\right)=-\delta \left(\sqrt{{\alpha }^2-{(\beta +iu)}^2}-\sqrt{{\alpha }^2-{\beta }^2}\right).$$

(21)

It follows that the NIG convexity adjustment reads

$${\omega }_{NIG}=-{\psi }_{NIG}(-i)=\delta \left(\sqrt{{\alpha }^2-{(\beta +1)}^2}-\sqrt{{\alpha }^2-{\beta }^2}\right).$$

(22)

Remark 2.

When $\alpha \to \infty$ (large steepness regime), then ${\omega }_{NIG}\to -\frac{{\sigma }^2}{2}(1+2\beta )$ where ${\sigma }^2:=\frac{\delta }{\alpha }$; when, furthermore, $\beta =0$ (symmetric process) then one recovers the usual Gaussian adjustment $-\frac{{\sigma }^2}{2}$ and the exponential NIG model degenerates into the Black–Scholes model.

3.2.2. Financial Applications

The exponential NIG model has been proved to provide a distinguished fitting to financial data many times. Let us mention, among others, initial tests for daily returns on Danish and German markets in Rydberg (1997) and, subsequently, on the FTSE All-share index (also known as “Actuaries index”) in Venter and de Jongh (2002). More recently, the impact of high-frequency trading has also been taken into account, and calibrations have been performed on intraday returns, e.g., in Figueroa-López et al. (2012) for different sampling frequencies. Like in the VG case, multivariate extensions have also been considered (see Luciano and Semeraro (2010) and references therein), and applications to credit risk have also been provided (Luciano 2009).

In Figure 2, we display the log-return density for a VG and NIG example, each being recovered from their characteristic functions while using the method of Kirkby (2015). While both models exhibit heavy-tails, the VG model is characterized by a pronounced cusp, especially for shorter maturities. This near singular behavior presents challenges for Fourier pricing methods, and techniques, such as spectral filtering, have been proposed as a remedy Cui et al. (2017); Phelan et al. (2019); Ruijter et al. (2015). In contrast, the closed form pricing formulas presented here exhibit smooth exponential convergence without special handling, as demonstrated in Section 4.

3.3. Fractional Diffusion Model

The FD model by Kleinert and Korbel (2016); Korbel and Luchko (2016) aims at generalizing the Lévy stable model by introducing a time-fractional derivative in the equation governing the probability densities, whose order will be interpreted as a subordination parameter. Before introducing the model and its characteristics, we briefly recall some basics of stable distributions and their link with fractional calculus.

3.3.1. Lévy-Stable Processes and Fractional Derivatives

Taking ${\lambda }_{+}={\lambda }_{-}=0$ and ${\alpha }_{+}={\alpha }_{-}:=\alpha \in (0,2)$ in (4) yields the Lévy measure of the stable (or $\alpha$-stable) process ${\left\{X_t^{\left(stable\right)}\right\}}_{t\in [0,T]}$:

$${\Pi }_{stable}\left(x\right)=\frac{c_{-}}{{\left|x\right|}^{1+\alpha }}{\mathrm{𝟙}}_{\{x<0\}}+\frac{c_{+}}{x^{1+\alpha }}{\mathrm{𝟙}}_{\{x>0\}}.$$

(23)

It is known that, when using Feller’s parametrization,

$$\left\{\begin{array}{c}{\sigma }^{\alpha }:=-(c_{+}+c_{-})\Gamma (-\alpha )cos\frac{\pi \alpha }{2}\hfill \\ \beta :=\frac{c_{+}-c_{-}}{c_{+}+c_{-}}\hfill \end{array}\right.$$

(24)

then the Lévy symbol of the stable process can be written as

$${\psi }_{stable}\left(u\right)={\sigma }^{\alpha }{\left|u\right|}^{\alpha }\left(1-i\beta tan\frac{\alpha \pi }{2}\mathrm{sgn}u\right)+iau$$

(25)

where the drift term a equals the expectation ${\mathbb{E}}^{\mathbb{Q}}\left[X_t^{\left(stable\right)}\right]$ as soon as $\alpha \in (1,2)$, that is, for the class of stable Paretian distributions. This class is the one with the greatest financial meaning and historical importance, with initial calibrations going back to Mandelbrot (1963) with $\alpha \simeq 1.7$ for cotton prices; see a comprehensive overview of these distributions in Zolotarev (1986), and of their financial applications in Mittnik and Rachev (2000). However, we may note that, due to the polynomial decay of (23) on the positive axis, the moment generating function and moments of all order do not exist unless $c_{+}=0$, or, equivalently, in terms of the parametrization (24), $\beta =-1$. This condition, which is known as the maximal negative asymmetry (or skewness) hypothesis, is the key assumption in the Finite Moment Log Stable (FMLS) model by Carr and Wu (2003); it follows from (25) that the FMLS martingale adjustment is

$${\omega }_{FMLS}=\frac{{\left(\frac{\sigma }{\sqrt{2}}\right)}^{\alpha }}{cos\frac{\pi \alpha }{2}}=\frac{1}{\pi }\Gamma \left(\frac{1-\alpha }{2}\right)\Gamma \left(\frac{1+\alpha }{2}\right){\left(\frac{\sigma }{\sqrt{2}}\right)}^{\alpha }$$

(26)

where we have introduced the $\sqrt{2}$ normalization, so as to recover the Gaussian adjustment when $\alpha =2$, and where the second equality is a consequence of the reflection formula for the Gamma function; in the limiting case $\alpha =2$, the stable distribution degenerates into the Normal one, and, therefore, the FMLS model recovers the Black–Scholes model.

Another consequence of the characteristic exponent (25) with $\beta =-1$ is that the FMLS probability density $f_{FMLS}(x,t)$ satisfies the space fractional diffusion equation

$$\frac{\partial f_{FMLS}}{\partial t}+{\omega }_{\alpha }{\mathcal{D}}_x^{\alpha }{f}_{FMLS}(x,t)=0,x\in \mathbb{R},t\in [0,T],$$

(27)

where ${\mathcal{D}}_x^{\alpha }:{=}^{\alpha -2}{\mathcal{D}}_x^{\alpha }$ is a particular case of the Riesz–Feller derivative defined (via its Fourier transform) by

$$\widehat{{}^{\theta }{\mathcal{D}}_x^{\alpha }f}\left(u\right)={\left|k\right|}^{\alpha }{e}^{i\left(\mathrm{sign}u\right)\theta \pi /2}\widehat{f}\left(u\right),\left|\theta \right|\le \min \{\alpha ,\alpha -2\}.$$

(28)

Remark 3.

When $\theta =0$, the Riesz–Feller derivative is simply called Riesz derivative (as the operator inverse to the Riesz potential, see all details and definitions e.g., in the classical monograph Samko et al. (1993)); the choice $\theta =\alpha -2$ in Equation (27) is the fractional calculus analogue to the maximal negative asymmetry hypothesis $\beta =-1$. When $\alpha =2$ then the Riesz-Feller derivative degenerates into the usual second derivative; in that case, (27) becomes the usual heat equation, whose fundamental solution (the heat kernel) is the probability density of the Wiener process.

3.3.2. Model Characteristics

The FD model generalizes the FMLS model, by allowing the time derivative in the Equation (27) for the probability density to be fractional as well:

$$\left({}^{*}{\mathcal{D}}_t^{\gamma }+{\omega }_{FD}{\mathcal{D}}_x^{\alpha }\right)f_{FD}(x,t)=0,x\in \mathbb{R},t\in [0,T],$$

(29)

where $\alpha \in (1,2]$, $\gamma \in (0,\alpha ]$ and ${}^{*}{\mathcal{D}}_{\gamma }^t$ denotes the Caputo fractional derivative (see the definitions and properties in Li et al. (2011) for instance); when $\gamma =1$, it coincides with the usual first-order derivative.

The fundamental solution to (29) has been determined in Mainardi et al. (2001) and admits the following Mellin–Barnes representation

$$f_{FD}(x,t)=\frac{1}{\alpha x}\underset{c-i\infty }{\overset{c+i\infty }{\int }}\left(G_{+}^{*}\left(s\right){\mathrm{𝟙}}_{\{x>0\}}+G_{-}^{*}\left(s\right){\mathrm{𝟙}}_{\{x<0\}}\right){\left(\frac{\left|x\right|}{{(-{\mu }_{\alpha ,\gamma }{t}^{\gamma })}^{\frac{1}{\alpha }}}\right)}^s\frac{\mathrm{d}s}{2i\pi },$$

(30)

where we have defined

$$G_{+}^{*}\left(s\right):=\frac{\Gamma (1-s)}{\Gamma (1-\gamma \frac{s}{\alpha })},G_{-}^{*}\left(s\right):=\frac{\Gamma \left(\frac{s}{\alpha }\right)\Gamma \left(1-\frac{s}{\alpha }\right)\Gamma (1-s)}{\Gamma \left(1-\frac{\gamma }{\alpha }s\right)\Gamma \left(\frac{\alpha -1}{\alpha }s\right)\Gamma \left(1-\frac{\alpha -1}{\alpha }s\right)}.$$

(31)

By analogy with (10), the price of a contingent claim $\mathcal{C}$ delivering a path independent payoff $\mathcal{P}$ at its maturity is defined to be

$$\mathcal{C}=e^{-r\tau }\underset{-\infty }{\overset{+\infty }{\int }}\mathcal{P}(S_t{e}^{(r-q+{\omega }_{FD})\tau +x},K_1,\dots ,K_n)f_{FD}(x,\tau )\mathrm{d}x$$

(32)

and the martingale adjustment is defined in terms of the cumulant generating function by

$${\omega }_{FD}=-log\underset{-\infty }{\overset{+\infty }{\int }}{e}^x{f}_{FD}(x,1)\mathrm{d}x.$$

(33)

It has been shown in Aguilar et al. (2018) that ${\omega }_{FD}$ can be conveniently expressed in the form of a series

$${\omega }_{FD}=-log\sum _{n=0}^{\infty }\frac{{(-1)}^n}{n!}\frac{\Gamma (1+\alpha n)}{\Gamma (1+\gamma \alpha n)}{\omega }_{FMLS}^n$$

(34)

where ${\omega }_{FMLS}$ is the FMLS martingale adjustment that is defined in (26), and under the condition that $\gamma \in (1-\frac{1}{\alpha },\alpha )$; using a first-order Taylor expansion for $log(1+u)$ and the expression (26), we obtain the useful approximation:

$${\omega }_{FD}=\frac{1}{\pi }\frac{\Gamma (1+\alpha )\Gamma \left(\frac{1-\alpha }{2}\right)\Gamma \left(\frac{1+\alpha }{2}\right)}{\Gamma (1+\gamma \alpha )}{\left(\frac{\sigma }{\sqrt{2}}\right)}^{\alpha }+O\left({\sigma }^{2\alpha }\right).$$

(35)

Remark 4.

When $\alpha =2$, we are left with

$${\omega }_{FD}=-\frac{{\sigma }^2}{\Gamma (1+2\gamma )}+O\left({\sigma }^4\right)$$

(36)

which coincides with the Black–Scholes adjustment $-\frac{{\sigma }^2}{2}$ when $\gamma =1$; we call the situation $\alpha =2$ and $\gamma \in (0,2]$ the “subordinated Black-Scholes” (sub-BS) model. This is a slight abuse of terminology, because we are not directly in the presence of a subordinating process (i.e., a non decreasing Lévy process) like in the VG and the NIG cases; subordination is achieved here via the introduction of a fractional time derivative whose order γ acts as a supplementary degree of freedom in the time dynamics. When $\gamma =1$, then the sub-BS model recovers the Black–Scholes model, like the VG model with $\nu \to 0$ and the NIG model with $\alpha \to \infty$; we summarize the situation in

Table 1. Some subordinated market models, and their limiting cases. Time changed models (exponential VG and NIG), FD, and sub-BS models recover the Black–Scholes (BS) model for specific values of their subordination parameters.

Subordinated Model	Limiting Regimes
VG$(\sigma ,\nu ,\theta )$	VG$(\sigma ,\nu ,0)\stackrel{\nu \to 0}{⟶}$ BS$\left(\sigma \right)$
NIG$(\alpha ,\beta ,\delta )$	NIG$(\alpha ,0,\delta )\stackrel{\alpha \to \infty }{⟶}$ BS$\left(\sqrt{\delta /\alpha }\right)$
FD$(\sigma ,\alpha ,\gamma )$	FD$(\sigma ,\alpha ,\gamma )\stackrel{\gamma \to 1}{⟶}$ FMLS$(\sigma ,\alpha )\stackrel{\alpha \to 2}{⟶}$ BS$\left(\sigma \right)$
sub-BS($\sigma$,$\gamma$) $:=$ FD$(\sigma ,2,\gamma )$	sub-BS$(\sigma ,\gamma )\stackrel{\gamma \to 1}{⟶}$ BS$\left(\sigma \right)$

.

3.3.3. Financial Applications

The purpose of the FD model is to allow more flexibility in the risk redistribution of returns. When the tail index $\alpha$ departs from 2, the model shifts the returns towards more significant losses due to the presence of the left fat tail. Similarly, when the order of the fractional derivative $\gamma$ is different from 1, the risk is shifted either towards shorter ($\gamma <1$) or longer ($\gamma >1$) maturities: for instance, when $\gamma <1$, the prices of short term options increase while the prices of long term options slightly decrease (see details, e.g., in Korbel and Luchko (2016)). This behavior is particularly relevant in periods of stressed market conditions: for instance, it has been observed during the 2020 market turmoils (consecutive to the COVID19 pandemics) that short-term implied volatility on the Euro Stoxx 50 index had increased sharply while remaining more stable for long-term options. Several calibrations have been made, notably on market data from turbulent times; in particular, in Kleinert and Korbel (2016), the FD model has been calibrated on data from S&P 500 options traded during November 2008. Such calibrations have shown that $\alpha$ could be quite different from 2: typically, $\alpha \simeq$ 1.6–1.7, that turns out to be quite remote from the log-normal hypothesis ($\alpha =2$), but relatively close to the initial Mandelbrot estimate for the stable law on cotton futures market (Mandelbrot 1963). Moreover, it has been noted that both fractional parameters $\alpha$ and $\gamma$ appeared to vary simultaneously and in the same direction, leaving the diffusion scaling exponent $\gamma /\alpha$ relatively stable.

4. Pricing and Volatility Modelling

In this section, we first recall the pricing formulas that were obtained in recent works for the VG, NIG, and FD models, in the case of a European option C delivering a payoff equal to ${[S_T-K]}^{+}$ at maturity. Then, we discuss some volatility properties when asset prices are not far from the money. In all of the following, we will denote the forward strike price and the log-forward moneyness by

$$F:=K{e}^{-r\tau },k:=log\frac{S_t}{K}+(r-q)\tau .$$

(37)

We will also use the notations

$$k_{VG}:=k+{\omega }_{VG}\tau ,k_{NIG}:=k+{\omega }_{NIG}\tau ,k_{FD}:=k+{\omega }_{FD}\tau$$

(38)

and, in the specific case of the VG model, we will denote

$${\tau }_{\nu }:=\frac{\tau }{\nu }-\frac{1}{2}{\sigma }_{\nu }:=\sigma \sqrt{\frac{\nu }{2}}.$$

(39)

We will assume that the underlying VG (resp. NIG) processes are symmetric, which is, $\theta =0$ (resp. $\beta =0$); this is to simplify the notations, but also because symmetric time-changed models extend the Black–Scholes setup (see Table 1). Therefore, we will be better able to compare the results with usual formulas known in the Black-Scholes model. We will also assume that ${\tau }_{\nu }\notin \mathbb{Q}$; note that this condition is not restrictive, due to the density of $\mathbb{Q}$ in $\mathbb{R}$: if ${\tau }_{\nu }\in \mathbb{Q}$, it is easy to make ${\tau }_{\nu }$ irrational by adding an arbitrary small perturbation, for instance, ${\tau }_{\nu }\to {\tau }_{\nu }+e/{10}^{10}$). Last, whenever the exponential NIG model is concerned, we will always assume that

$$\frac{|k_{NIG}|}{\delta \tau }<1$$

(40)

to ensure the convergence of the series. Please note that this condition is automatically satisfied when options are not far from the money because, in this case, $k_{NIG}$ is small. When $S_t$ is far from K, condition (40) necessitates a restriction on options maturities in order to be satisfied; for a typical set of parameters, maturities shorter than two or three months should be excluded, which remains a reasonable limitation (see details in Aguilar (2020b)).

For convenience, we summarize the pricing formulas for a European call option.

Formula 1

(European call: pricing formulas).

(i) The value at time t of a European call option in the exponential VG model is:

-

(OTM price) If $k_{VG}<0$,

$$\begin{array}{c}{C}_{VG}^{-}(k_{VG},{\sigma }_{\nu })=\frac{F}{2\Gamma \left(\frac{\tau }{\nu }\right)}\sum _{\begin{array}{c}{n}_1=0\\ n_2=1\end{array}}^{\infty }\frac{{(-1)}^{n_1}}{n_1!}\left[\frac{\Gamma (\frac{-n_1+n_2+1}{2}+{\tau }_{\nu })}{\Gamma (\frac{-n_1+n_2}{2}+1)}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^{n_1}{\sigma }_{\nu }^{n_2}\right.\hfill \\ \hfill \left.+2\frac{\Gamma (-2n_1-n_2-1-2{\tau }_{\nu })}{\Gamma (-n_1+\frac{1}{2}-{\tau }_{\nu })}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^{2n_1+1+2{\tau }_{\nu }}{(-k_{VG})}^{n_2}\right].\end{array}$$

(41)

-

(ITM price) If $k_{VG}>0$,

$$C_{VG}^{+}(k_{VG},{\sigma }_{\nu })=S_t{e}^{-q\tau }-F-C_{VG}^{-}(k_{VG},-{\sigma }_{\nu }).$$

(42)

-

(ATM price) If $k_{VG}=0$,

$$C_{VG}^{-}(k_{VG},{\sigma }_{\nu })=C_{VG}^{+}(k_{VG},{\sigma }_{\nu })=\frac{F}{2\Gamma \left(\frac{\tau }{\nu }\right)}\sum _{n=1}^{\infty }\frac{\Gamma (\frac{n+1}{2}+{\tau }_{\nu })}{\Gamma (\frac{n}{2}+1)}{\sigma }_{\nu }^n.$$

(43)

(ii) The value at time t of a European call option in the exponential NIG model is:

$$C_{NIG}=\frac{F\alpha e^{\alpha \delta \tau }}{\sqrt{\pi }}\sum _{\begin{array}{c}{n}_1=0\\ n_2=1\end{array}}^{\infty }\frac{k_{NIG}^{n_1}}{n_1!\Gamma (1+\frac{-n_1+n_2}{2})}{K}_{\frac{n_1-n_2+1}{2}}\left(\alpha \delta \tau \right){\left(\frac{\delta \tau }{2\alpha }\right)}^{\frac{-n_1+n_2+1}{2}}.$$

(44)

(iii) The value at time t of a European call option in the FD model is:

$$C_{FD}=\frac{F}{\alpha }\sum _{\begin{array}{c}{n}_1=0\\ n_2=1\end{array}}^{\infty }\frac{k_{FD}^{n_1}}{n_1!\Gamma (1+\gamma \frac{-n_1+n_2}{\alpha })}{(-{\omega }_{FD}{\tau }^{\gamma })}^{\frac{-n_1+n_2}{\alpha }}.$$

(45)

Proof. 

(i) is proved in Aguilar (2020a), (ii) is proved in Aguilar (2020b) and (iii) in Aguilar et al. (2018). □

The pricing formulas in Formula (1) converge exponentially fast to the true prices, as exhibited in Figure 3 for the VG model (Left). The reference prices are obtained with $N=50$ terms, and they are verified by the method of Kirkby (2015). Recalling Remark 1, VG$(\sigma ,\nu ,0)\stackrel{\nu \to 0}{⟶}$ BS$\left(\sigma \right)$, and fewer terms are required to accurately price the option as $\nu \to 0$. A nearly identical convergence profile is observed for NIG (Right), displayed for several values of $\alpha$. Recalling that NIG$(\alpha ,0,\delta )\stackrel{\alpha \to \infty }{⟶}$ BS$\left(\sqrt{\delta /\alpha }\right)$ (see Table 1), we again see that fewer terms are required in order to accurately price under NIG for larger $\alpha$.

4.1. At-the-Money Forward Approximations

Let us assume that options are at-the-money forward (ATMF), that is, $S_t=F$ (or, equivalently, $k=0$). If we approximate the European call price in the exponential VG model by the first term of the series (41), which is, the term for $n_1=0$, $n_2=1$, we obtain (recall that $\Gamma (3/2)=\sqrt{\pi }/2$):

$$C_{VG}^{-}=\frac{S_t}{\sqrt{2\pi }}\frac{\Gamma (\frac{1}{2}+\frac{\tau }{\nu })}{\Gamma \left(\frac{\tau }{\nu }\right)}\sigma \sqrt{\nu }.$$

(46)

Using the Stirling approximation for the Gamma function, we know that

$$\frac{\Gamma \left(\frac{1}{2}+\frac{\tau }{\nu }\right)}{\Gamma \left(\frac{\tau }{\nu }\right)}\underset{\nu \to 0}{\sim }\sqrt{\frac{\tau }{\nu }},$$

(47)

and, therefore, in the low variance regime, we recover the well-known approximation for the ATMF Black–Scholes price

$$C_{BS}\simeq \frac{S_t}{\sqrt{2\pi }}\sigma \sqrt{\tau }.$$

(48)

The approximation (48) is often known to market practitioners under the form $C\simeq 0.4S\sigma \sqrt{\tau }$, because $1/\sqrt{2\pi }=0.399\dots$, and it was first derived in Brenner and Subrahmanyam (1994).

In a similar way, the ATMF price in the exponential NIG model can be approximated by the first term of the series (44), resulting in

$$C_{NIG}=\frac{S_t\delta \tau e^{\alpha \delta \tau }}{\pi }{K}_0\left(\alpha \delta \tau \right).$$

(49)

Using the large argument behavior of the Bessel function (see Abramowitz and Stegun (1972))

$$K_0\left(\alpha \delta \tau \right)\underset{\alpha \to \infty }{\sim }\sqrt{\frac{\pi }{2\alpha \delta \tau }}{e}^{-\alpha \delta \tau },$$

(50)

it is immediate to see that (49) also recovers the approximation (48) in the large steepness regime, with ${\sigma }^2:=\delta /\alpha$. Likewise, in the FD model, the ATMF price is approximated by the first term of the series (45), resulting in

$$C_{FD}=\frac{S_t}{\alpha }\frac{1}{\Gamma (1+\frac{\gamma }{\alpha })}{(-{\omega }_{FD}{\tau }^{\gamma })}^{\frac{1}{\alpha }}.$$

(51)

Taking $\alpha =2$ and using (36), the ATMF price in the sub-BS model becomes

$$C_{sub-BS}=\frac{S_t}{2}\frac{1}{\Gamma (1+\frac{\gamma }{2})\sqrt{\Gamma (1+2\gamma )}}\sigma {\tau }^{\frac{\gamma }{2}}$$

(52)

which recovers (48) when $\gamma \to 1$.

4.2. Implied Volatility

One key benefit of the subordinated models is their ability to capture the heavy-tails that were observed in financial markets. For the VG model, $\nu$ directly controls the tail-heaviness, as illustrated in Figure 4. In particular, large values of $\nu$ lead to steep implied volatility smiles. The ATMF prices that are obtained in Section 4.1 are helpful to approximate the implied volatility ${\sigma }_I$ of the subordinated models, when $S_t$ is close to F. Denoting by $C_t$ the market price of an ATMF European call option at time t and inverting (46), we immediately see that the VG implied volatility is

$${\sigma }_{VG}=\sqrt{\frac{2\pi }{\nu }}\frac{\Gamma \left(\frac{\tau }{\nu }\right)}{\Gamma (\frac{1}{2}+\frac{\tau }{\nu })}\frac{C_t}{S_t},$$

(53)

and, similarly, inverting Equation (52), the sub-BS implied volatility is

$${\sigma }_{sub-BS}=\frac{2\sqrt{\Gamma (1+2\gamma )}\Gamma (1+\frac{\gamma }{2})}{{\tau }^{\frac{\gamma }{2}}}\frac{C_t}{S_t}.$$

(54)

As expected, VG and sub-BS both implied volatilities recover the BS implied volatility in their limiting regimes ($\nu \to 0$ and $\gamma \to 1$):

$${\sigma }_{BS}=\sqrt{\frac{2\pi }{\tau }}\frac{C_t}{S_t}.$$

(55)

In the NIG case, things are a bit more complicated, because one has to solve

$$\frac{S_t\delta \tau e^{\alpha \delta \tau }}{\pi }{K}_0\left(\alpha \delta \tau \right)=C_t$$

(56)

for which there is no exact solution in an analytical form. Nevertheless, an analytical approximation can be determined by using Hankel’s expression for the Bessel function (see Andrews (1992) or any monograph on special functions), which goes, as follows: define, for $\rho \in \mathbb{R}$,

$$\left\{\begin{array}{cc}\hfill & a_0\left(\rho \right)=1\hfill \\ \hfill & a_k\left(\rho \right)=\frac{(4{\rho }^2-1^2)(4{\rho }^2-3^2)\dots (4{\rho }^2-{(2k-1)}^2)}{k!8^k},k\ge 1,\hfill \end{array}\right.$$

(57)

then, for large z and fixed $\rho$, we have:

$$\begin{array}{cc}\hfill K_{\rho }\left(z\right)& \underset{z\to \infty }{=}\sqrt{\frac{\pi }{2z}}{e}^{-z}\sum _{k=0}^{\infty }\frac{a_k\left(\rho \right)}{z^k}.\hfill \end{array}$$

(58)

In particular, when $4{\rho }^2-1=0$, i.e., when $\rho =\frac{1}{2}$, all the $a_k\left(\rho \right)$ are null in definition (57) when $k\ge 1$, and we are left with:

$$K_{\frac{1}{2}}\left(z\right)=\sqrt{\frac{\pi }{2z}}{e}^{-z}$$

(59)

for all z. Using (58) up to $k=1$ for $z=\alpha \delta \tau$ and inserting into (56), we are left with the quadratic equation

$$X^2-\alpha \sqrt{2\pi }\frac{C_t}{S_t}X-\frac{1}{8}=0(X:=\sqrt{\alpha \delta \tau }),$$

(60)

whose positive solution reads

$$X=\frac{1}{2}\left(\alpha \sqrt{2\pi }\frac{C_t}{S_t}+\sqrt{2\pi {\alpha }^2\frac{C_t^2}{S_t^2}+\frac{1}{2}}\right).$$

(61)

Taylor expanding for large $\alpha$ and turning back to the $\delta$ variable, we obtain

$$\delta =\frac{2\pi \alpha }{\tau }\frac{C_t^2}{S_t^2}+\frac{1}{4\alpha \tau }+O\left(\frac{1}{{\alpha }^3}\right)$$

(62)

which, at first order, recovers (55) for ${\sigma }^2:=\delta /\alpha$. We summarize these results in

Table 2. Volatility modelling for ATMF options in various subordinated models, and their limiting cases.

ATMF Implied Volatility (European Options)
Exponential VG	${\sigma }_{VG}=\sqrt{\frac{2\pi }{\nu }}\frac{\Gamma \left(\frac{\tau }{\nu }\right)}{\Gamma (\frac{1}{2}+\frac{\tau }{\nu })}\frac{C_t}{S_t}$
	Low variance regime ($\nu \to 0$):
	${\sigma }_{VG}\to \sqrt{\frac{2\pi }{\tau }}\frac{C_t}{S_t}$
Exponential NIG	Solve   $\frac{S_t\delta \tau e^{\alpha \delta \tau }}{\pi }{K}_0\left(\alpha \delta \tau \right)=C_t$
	Large steepness regime ($\alpha \to \infty$):
	$\delta =\frac{2\pi \alpha }{\tau }\frac{C_t^2}{S_t^2}+\frac{1}{4\alpha \tau }+O\left(\frac{1}{{\alpha }^2}\right)$
	At order ${\alpha }^0$:
	${\sigma }_{NIG}:=\sqrt{\frac{\delta }{\alpha }}=\sqrt{\frac{2\pi }{\tau }}\frac{C_t}{S_t}$
sub-BS	${\sigma }_{FD}=\frac{2\sqrt{\Gamma (1+2\gamma )}\Gamma (1+\frac{\gamma }{2})}{{\tau }^{\frac{\gamma }{2}}}\frac{C_t}{S_t}$
	Non-fractional regime ($\gamma \to 1$):
	${\sigma }_{FD}\to \sqrt{\frac{2\pi }{\tau }}\frac{C_t}{S_t}$

.

5. First-Order Sensitivities

The sensitivity of a contingent claim $\mathcal{C}$ to the underlying asset, often denoted Delta or $\Delta$, is defined by $\Delta :=\partial \mathcal{C}/\partial S_t$; by deriving Formula (1) with respect to $S_t$ and re-arranging the terms, we obtain the following expressions for European options in subordinated market models:

Formula 2

(European call: Delta).

(i) The Delta at time t of a European call option in the exponential VG model is:

-

(OTM sensitivity) If $k_{VG}<0$,

$$\begin{array}{c}{\Delta }_{VG}^{-}(k_{VG},{\sigma }_{\nu })=\frac{F}{2S_t\Gamma \left(\frac{\tau }{\nu }\right)}\sum _{\begin{array}{c}{n}_1=0\\ n_2=1\end{array}}^{\infty }\frac{{(-1)}^{n_1}}{n_1!}\left[-\frac{n_1\Gamma (\frac{-n_1+n_2+1}{2}+{\tau }_{\nu })}{\Gamma (\frac{-n_1+n_2}{2}+1)}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^{n_1-1}{\sigma }_{\nu }^{n_2-1}\right.\hfill \\ \hfill \left.+2\frac{\Gamma (-2n_1-n_2-2{\tau }_{\nu })}{\Gamma (-n_1+\frac{1}{2}-{\tau }_{\nu })}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^{2n_1+1+2{\tau }_{\nu }}{(-k_{VG})}^{n_2-1}\right].\end{array}$$

(63)

-

(ITM sensitivity) If $k_{VG}>0$,

$${\Delta }_{VG}^{+}(k_{VG},{\sigma }_{\nu })=e^{-q\tau }-{\Delta }_{VG}^{-}(k_{VG},-{\sigma }_{\nu }).$$

(64)

-

(ATM sensitivity) If $k_{VG}=0$,

$${\Delta }_{VG}^{-}(k_{VG},{\sigma }_{\nu })={\Delta }_{VG}^{+}(k_{VG},{\sigma }_{\nu })=\frac{F}{2S_t\Gamma \left(\frac{\tau }{\nu }\right)}\sum _{n=1}^{\infty }\frac{\Gamma (\frac{n}{2}+{\tau }_{\nu })}{\Gamma \left(\frac{n+1}{2}\right)}{\sigma }_{\nu }^{n-1}.$$

(65)

(ii) The Delta at time t of a European call option in the exponential NIG model is:

$${\Delta }_{NIG}=\frac{F\alpha e^{\alpha \delta \tau }}{S_t\sqrt{\pi }}\sum _{\begin{array}{c}{n}_1=0\\ n_2=1\end{array}}^{\infty }\frac{k_{NIG}^{n_1}}{n_1!\Gamma \left(\frac{-n_1+n_2+1}{2}\right)}{K}_{\frac{n_1-n_2}{2}+1}\left(\alpha \delta \tau \right){\left(\frac{\delta \tau }{2\alpha }\right)}^{\frac{-n_1+n_2}{2}}.$$

(66)

(iii) The Delta at time t of a European call option in the FD model is:

$${\Delta }_{FD}=\frac{F}{\alpha S_t}\sum _{\begin{array}{c}{n}_1=0\\ n_2=0\end{array}}^{\infty }\frac{k_{FD}^{n_1}}{n_1!\Gamma (1+\gamma \frac{-n_1+n_2}{\alpha })}{(-{\omega }_{FD}{\tau }^{\gamma })}^{\frac{-n_1+n_2}{\alpha }}.$$

(67)

For illustration, in Figure 5 we compare the Delta of VG$(\sigma ,\nu ,0)$ while using Formula (2) with that of BS$\left(\sigma \right)$. Similarly, we compare the Dollar Gamma using Formula (3). Figure 6 provides a comparison for NIG$(\alpha ,0,\delta )$. For both models, we can see the substantial impact of the heavy-tailed assumption and its implications for hedging. In the next section, we discuss delta hedging in more detail, and provide some simplified approximations for the ATMF case.

5.1. Delta Hedging

The leading term for the Delta of the European option in the exponential VG case is given for $n_1=n_2=1$ in (63) and it reads:

$${\Delta }_{VG}^{-}=\frac{F}{2S_t\Gamma \left(\frac{\tau }{\nu }\right)}\left[\Gamma \left(\frac{1}{2}+{\tau }_{\nu }\right)-2\frac{\Gamma (-3-3{\tau }_{\nu })}{\Gamma (-\frac{1}{2}-{\tau }_{\nu })}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^{-3+2{\tau }_{\nu }}\right];$$

(68)

In the ATMF situation ($S_t=F$, i.e., $k=0$), we can re-write the martingale adjustment and Taylor expand for small $\sigma$

$$k_{VG}={\omega }_{VG}\tau =-\frac{{\sigma }^2}{2}\tau +O\left({\sigma }^4\right)$$

(69)

and, recalling ${\tau }_{\nu }=\frac{\tau }{\nu }-\frac{1}{2}$, we obtain

$${\Delta }_{VG}^{-}=\frac{1}{2}+O\left({\sigma }^{4({\tau }_{\nu }+1)}\right).$$

(70)

It is interesting to note that, in the ATMF situation, ${\Delta }_{VG}^{-}\simeq \frac{1}{2}$, which is, it suffices to be long one unit of the asset S and short two units of an European call written on this asset, to offset the impact of the variations of S on a portfolio; this fact is well-known in the usual Black–Scholes theory, and is therefore preserved in the exponential VG model. The same observation actually also holds for the exponential NIG model: indeed, the leading term in the series (66) (obtained for $n_1=0$, $n_2=1$) reads

$$\frac{F\alpha e^{\alpha \delta \tau }}{S_t\sqrt{\pi }}{K}_{\frac{1}{2}}\left(\alpha \delta \tau \right)\sqrt{\frac{\delta \tau }{2\alpha }}=\frac{F}{2S_t}$$

(71)

where we have used the particular value of the Bessel function of index $\frac{1}{2}$ (59) in order to simplify the expression; when $F=S_t$, we obtain ${\Delta }_{NIG}=\frac{1}{2}$, which, again, turns out to be similar to the usual Black–Scholes behavior. This effect is clearly illustrated in Figure 5 and Figure 6. We can conclude that in both the exponential VG and NIG models, the presence of a time subordination does not modify the delta hedging policy, at least when options are not far from the money. In contrast, the option Gamma is significantly influenced by time subordination, and it is discussed further in Section 6.

In the FD model, things are a bit different; keeping only the leading term ($n_1=n_2=0$) in (67) yields

$${\Delta }_{FD}=\frac{F}{\alpha S_t}\stackrel{S_t\to F}{⟶}\frac{1}{\alpha }$$

(72)

which explicitly depends on the tail parameter $\alpha$ and resumes to $\frac{1}{2}$ in the sub-BS model, for any subordination parameter $\gamma$. Again, the subordination parameter plays no role in the delta-hedging policy of the portfolio, which is entirely governed by the tail parameter $\alpha$; in other words, it suffices to be long one unit of the underlying asset S and short $\alpha$ European calls to offset the effect of the variations of S on the portfolio’s value.

5.2. Comparisons with Numerical Techniques

In this subsection, we show that the series formulas for the first order sensitivity $\Delta$ provided by Formula (2) are a very efficient alternative to Fourier-based computations. Such calculations are typically based on a representation for the price of an European call in terms of Arrau–Debreu securities (see e.g., Lewis (2001))

$${\mathbb{E}}^{\mathbb{Q}}\left[e^{-r\tau }{(S_T-K)}^{+}|S_t\right]=S_t{e}^{-q\tau }{\Pi }_1-K{e}^{-r\tau }{\Pi }_2,$$

(73)

where $e^{-q\tau }{\Pi }_1$ is the option’s Delta. This quantity is known to admit a convenient representation in the Fourier space:

$$\Delta =e^{-q\tau }{\Pi }_1=\frac{1}{2}+\frac{1}{\pi }\underset{0}{\overset{\infty }{\int }}Re\left[\frac{e^{iuk}\tilde{\Psi }(u-i,\tau )}{iu}\right]\mathrm{d}u$$

(74)

where k is the log forward moneyness defined in (37), and the “risk-neutralized” characteristic function is defined by:

$$\tilde{\Psi }(u,t):=e^{iu\omega t}\Psi (u,t)=e^{(iu\omega +\psi (u\left)\right)t}$$

(75)

In the case of the exponential VG and NIG models, the integral in (74) can be easily carried out by inserting the expressions for the Lévy symbol $\psi \left(u\right)$ and the martingale adjustment $\omega$, and by evaluating the integral by any classical recursive algorithm (such as a simple trapezoidal rule, for instance). In

Table 3. First order sensitivity (Delta) of European call options in the exponential VG and NIG models, obtained by truncations of Formula (2), and by a numerical evaluation of (74). Here, $N=n_1=n_2$ is the number of terms in the truncated series. Parameters: $K=4000$, $r=1\%$, $q=0\%$, $\tau =1$.

Exponential VG Model [$\mathit{\sigma }=0.2$, $\mathit{\nu }=0.85$]
	Formula (2)				Lewis (74)
	$\mathit{N}=\mathbf{3}$	$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{10}$	$\mathit{N}=\mathbf{15}$
Deep OTM ($S_t=3000$)	2.1823	0.6347	0.0941	0.0940	0.0940
OTM ($S_t=3500$)	0.4113	0.2567	0.2455	0.2455	0.2455
ATM ($S_t=4040.90$)	0.5703	0.5718	0.5719	0.5719	0.5719
ITM ($S_t=4500$)	0.7569	0.8113	0.8134	0.8134	0.8134
Deep ITM ($S_t=5000$)	0.4729	0.8589	0.9206	0.9206	0.9206
Exponential NIG Model [$\mathbf{\alpha }=\mathbf{9}$, $\mathbf{\delta }=\mathbf{1}.\mathbf{2}$]
	Formula (2)				Lewis (74)
	$\mathit{N}=\mathbf{3}$	$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{10}$	$\mathit{N}=\mathbf{15}$
Deep OTM ($S_t=3000$)	0.2921	0.2722	0.2747	0.2748	0.2748
OTM ($S_t=3500$)	0.4289	0.4309	0.4311	0.4311	0.4311
ATM ($S_t=4234.09$)	0.6336	0.6410	0.6412	0.6412	0.6412
ITM ($S_t=4500$)	0.6936	0.7030	0.7033	0.7033	0.7033
Deep ITM ($S_t=5000$)	0.7827	0.7966	0.7971	0.7971	0.7971

, we compare the results of such numerical evaluations, with several truncations of the series in Formula (2), and for various market configurations. We can observe that the convergence is extremely accurate and fast, notably in the ATM region: this is because, in that case, $k\simeq 0$, which tends to accelerate the overall convergence of the series.

6. Second-Order Sensitivities and Portfolio Performance

6.1. Gamma, Dollar Gamma

The second order derivative of a contingent claim $\mathcal{C}$ with respect to $S_t$ is often denoted by $\Gamma :={\partial }^2\mathcal{C}/\partial S_t^2$. It is closely related to the performance, or Profit and Loss (P&L) of a portfolio: if $t_2>t_1$ are two trading days, then the P&L between $t_1$ and $t_2$ is

$$\mathrm{P}\&\mathrm{L}=\Theta \Delta t+\mathrm{market}\mathrm{P}\&\mathrm{L},\Delta t:=t_2-t_1,$$

(76)

where $\Theta$ is the time sensitivity of the portfolio, and the market P&L is, at order 2:

$$\mathrm{market}\mathrm{P}\&\mathrm{L}:=\Delta (\Delta S_t)+\frac{1}{2}\Gamma {(\Delta S_t)}^2\Delta S_t:=S_{t_2}-S_{t_1}.$$

(77)

Assuming that the portfolio has been delta-hedged, then we are left with

$$\mathrm{market}\mathrm{P}\&\mathrm{L}=\$\Gamma {\left(\frac{\Delta S_t}{S_t}\right)}^2$$

(78)

where the Dollar Gamma has been defined by $\$\Gamma :=\frac{1}{2}{S}_t^2\Gamma$; relation (78) is widely employed in financial engineering, because it allows for expressing the performance of the portfolio as a simple function of the realized variance of the underlying S. In the usual Black–Scholes theory, it is well known that the Dollar Gamma of the European call is

$$\${\Gamma }_{BS}=\frac{S_t}{2\sigma \sqrt{2\pi \tau }}.$$

(79)

Remarkably, as shown by Formula (3), the Gamma of European options admits a simple form in subordinated market models: while the series expansion for the price or the first-order sensitivity is expressed in terms of a double sum, the Gamma can be expressed as a sum over a single index.

Formula 3

(European call: Gamma).

(i) The Gamma at time t of a European call option in the exponential VG model is:

-

(OTM sensitivity) If $k_{VG}<0$,

$$\begin{array}{c}{\Gamma }_{VG}^{-}(k_{VG},{\sigma }_{\nu })=\frac{F}{2S_t^2{\sigma }_{\nu }\Gamma \left(\frac{\tau }{\nu }\right)}\sum _{n=0}^{\infty }\frac{{(-1)}^n}{n!}\left[\frac{\Gamma (-\frac{n}{2}+{\tau }_{\nu })}{\Gamma \left(\frac{-n+1}{2}\right)}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^n\right.\hfill \\ \hfill \left.+2\frac{\Gamma (-2n-2{\tau }_{\nu })}{\Gamma (-n+\frac{1}{2}-{\tau }_{\nu })}{\left(\frac{-k_{VG}}{{\sigma }_{\nu }}\right)}^{2n+2{\tau }_{\nu }}\right].\end{array}$$

(80)

-

(ITM sensitivity) If $k_{VG}>0$,

$${\Gamma }_{VG}^{+}(k_{VG},{\sigma }_{\nu })=-{\Gamma }_{VG}^{-}(k_{VG},-{\sigma }_{\nu }).$$

(81)

-

(ATM sensitivity) If $k_{VG}=0$,

$${\Gamma }_{VG}^{-}(k_{VG},{\sigma }_{\nu })={\Gamma }_{VG}^{+}(k_{VG},{\sigma }_{\nu })=\frac{F}{2\sqrt{\pi }{S}_t^2{\sigma }_{\nu }\Gamma \left(\frac{\tau }{\nu }\right)}\frac{\Gamma ({\tau }_{\nu }-\frac{1}{2})}{\Gamma \left(\frac{\tau }{\nu }\right)}.$$

(82)

(ii) The Gamma at time t of a European call option in the exponential NIG model is:

$${\Gamma }_{NIG}=\frac{F\alpha e^{\alpha \delta \tau }}{S_t^2\sqrt{\pi }}\sum _{n=0}^{\infty }\frac{k_{NIG}^n}{n!\Gamma \left(\frac{-n+1}{2}\right)}{K}_{\frac{n}{2}+1}\left(\alpha \delta \tau \right){\left(\frac{\delta \tau }{2\alpha }\right)}^{-\frac{n}{2}}.$$

(83)

(iii) The Gamma at time t of a European call option in the FD model is:

$${\Gamma }_{FD}=\frac{F}{\alpha S_t^2}\sum _{n=0}^{\infty }\frac{k_{FD}^n}{n!\Gamma (1-\frac{\gamma }{\alpha }(n+1))}{(-{\omega }_{FD}{\tau }^{\gamma })}^{-\frac{n+1}{\alpha }}.$$

(84)

Proof. 

The formulas are all straightforward to obtain, by deriving the series in Formula (2) with respect to $S_t$ and making an appropriate change of variables. For instance, in the NIG case, we have

$$\begin{array}{c}{\Gamma }_{NIG}=\frac{F\alpha e^{\alpha \delta \tau }}{S_t^2\sqrt{\pi }}\left[-\sum _{\begin{array}{c}{n}_1=0\\ n_2=1\end{array}}^{\infty }\frac{k_{NIG}^{n_1}}{n_1!\Gamma \left(\frac{-n_1+n_2+1}{2}\right)}{K}_{\frac{n_1-n_2}{2}+1}\left(\alpha \delta \tau \right){\left(\frac{\delta \tau }{2\alpha }\right)}^{\frac{-n_1+n_2}{2}}\right.\hfill \\ \hfill \left.+\sum _{\begin{array}{c}{n}_1=1\\ n_2=1\end{array}}^{\infty }\frac{k_{NIG}^{n_1-1}}{(n_1-1)!\Gamma \left(\frac{-n_1+n_2+1}{2}\right)}{K}_{\frac{n_1-n_2}{2}+1}\left(\alpha \delta \tau \right){\left(\frac{\delta \tau }{2\alpha }\right)}^{\frac{-n_1+n_2}{2}}\right].\end{array}$$

(85)

Performing the change of variables ${\tilde{n}}_1:=n_1+1$, ${\tilde{n}}_2\to n_2+1$ in the second sum shows that only the terms for ${\tilde{n}}_2=0$ survive; renaming ${\tilde{n}}_1:=n$ yields Formula (83). □

6.2. Properties and Particular Cases

Let us discuss some useful approximations and qualitative properties of Formula (3). First, in the VG case, the leading term ($n=0$) in (80) is

$$\frac{F}{2S_t^2{\sigma }_{\nu }\Gamma \left(\frac{\tau }{\nu }\right)}\left[\frac{\Gamma \left({\tau }_{\nu }\right)}{\sqrt{\pi }}+2\frac{\Gamma (-2{\tau }_{\nu })}{\Gamma (\frac{1}{2}-{\tau }_{\nu })}{\left(-\frac{k_{VG}}{{\sigma }_{\nu }}\right)}^{2{\tau }_{\nu }}\right].$$

(86)

Taylor expanding the VG martingale adjustment for small $\nu$ and assuming that we are not far from the money forward ($S_t\to F$), we have

$$k_{VG}\underset{\nu \to 0}{\sim }k-\frac{{\sigma }^2}{2}\tau \underset{S_t\to F}{\sim }-\frac{{\sigma }^2}{2}\tau ,$$

(87)

therefore, the Gamma writes, at first order:

$${\Gamma }_{VG}^{-}=\frac{1}{2\sqrt{\pi }{S}_t{\sigma }_{\nu }}\frac{\Gamma (\frac{\tau }{\nu }-\frac{1}{2})}{\Gamma \left(\frac{\tau }{\nu }\right)}$$

(88)

and the Dollar Gamma immediately follows:

$$\${\Gamma }_{VG}^{-}=\frac{S_t}{4\sqrt{\pi }{\sigma }_{\nu }}\frac{\Gamma (\frac{\tau }{\nu }-\frac{1}{2})}{\Gamma \left(\frac{\tau }{\nu }\right)}.$$

(89)

While using the functional relation $\Gamma (z+1)=z\Gamma \left(z\right)$ and the Stirling approximation (47), we have:

$$\frac{\Gamma (\frac{\tau }{\nu }-\frac{1}{2})}{\Gamma \left(\frac{\tau }{\nu }\right)}=\frac{1}{\frac{\tau }{\nu }-\frac{1}{2}}\frac{\Gamma (\frac{\tau }{\nu }+\frac{1}{2})}{\Gamma \left(\frac{\tau }{\nu }\right)}\underset{\nu \to 0}{\sim }\sqrt{\frac{\nu }{\tau }}$$

(90)

and, therefore, we obtain the behavior of (89) in the low variance regime:

$$\${\Gamma }_{VG}^{-}\stackrel{\nu \to 0}{⟶}\frac{S_t}{2\sigma \sqrt{2\pi \tau }},$$

(91)

thus recovering the Black–Scholes Dollar Gamma (79) in this limit. It is interesting to note that, contrary to the first order sensitivity Delta (70), which appeared to be independent of $\nu$; this is no longer the case with the second order sensitivity Gamma (89) that explicitly depends on the subordination parameter. In other words, while the subordination parameter does not modify the Delta Hedging policy of the portfolio (when not far from the money), it directly impacts its performance. This observation also holds in the exponential NIG model; indeed, the leading term in (83) for $S_t\to F$ reads

$${\Gamma }_{NIG}=\frac{\alpha e^{\alpha \delta \tau }}{\pi S_t}{K}_1\left(\alpha \delta \tau \right)$$

(92)

and therefore the Dollar Gamma is

$$\${\Gamma }_{NIG}=\frac{\alpha S_t{e}^{\alpha \delta \tau }}{2\pi }{K}_1\left(\alpha \delta \tau \right).$$

(93)

Using the asymptotic behavior for large argument (58) for the Bessel function, we know that

$$K_1\left(\alpha \delta \tau \right)\underset{\alpha \to \infty }{\sim }\sqrt{\frac{\pi }{2\alpha \delta \tau }}{e}^{-\alpha \delta \tau }$$

(94)

and, therefore

$$\${\Gamma }_{NIG}\stackrel{\alpha \to \infty }{⟶}\frac{S_t}{2\sigma \sqrt{2\pi \tau }},{\sigma }^2:=\frac{\delta }{\alpha },$$

(95)

recovering the Black–Scholes Dollar Gamma (79). Last in the FD model, the leading term in the series (84) for $S_t\to F$ is

$${\Gamma }_{FD}=\frac{1}{\alpha S_t}\frac{{(-{\omega }_{FD}{\tau }^{\gamma })}^{-\frac{1}{\alpha }}}{\Gamma (1-\frac{\gamma }{\alpha })},$$

(96)

and, in the sub-BS model ($\alpha =2$), using the approximation (36) for the martingale adjustment,

$${\Gamma }_{sub-BS}=\frac{1}{2S_t}\frac{\sqrt{\Gamma (1+2\gamma )}}{\Gamma (1-\frac{\gamma }{2})\sigma {\tau }^{\frac{\gamma }{2}}}.$$

(97)

Therefore, the Dollar Gamma in the sub-BS model is

$$\${\Gamma }_{sub-BS}=\frac{S_t}{4}\frac{\sqrt{\Gamma (1+2\gamma )}}{\Gamma (1-\frac{\gamma }{2})\sigma {\tau }^{\frac{\gamma }{2}}}.$$

(98)

and, in the non fractional limit ($\gamma \to 1$), we have, again,

$$\${\Gamma }_{sub-BS}\stackrel{\gamma \to 1}{⟶}\frac{S_t}{2\sigma \sqrt{2\pi \tau }}.$$

(99)

In

Table 4. First and second order market sensitivities (ATMF situation) for European call options in various subordinated models, and their limiting cases. Time subordination does not affect the Delta, but it directly impacts the Gamma of options.

	1st Order ($\mathbf{\Delta }$)	2nd Order ($\mathbf{\Gamma }$)
Exponential VG	$\frac{1}{2}$	$\frac{1}{2\sqrt{\pi }{S}_t{\sigma }_{\nu }}\frac{\Gamma (\frac{\tau }{\nu }-\frac{1}{2})}{\Gamma \left(\frac{\tau }{\nu }\right)}$
		Low variance regime ($\nu \to 0$):
		$\frac{1}{S_t\sigma \sqrt{2\pi \tau }}$
Exponential NIG	$\frac{1}{2}$	$\frac{\alpha e^{\alpha \delta \tau }}{\pi S_t}{K}_1\left(\alpha \delta \tau \right)$
		Large steepness regime ($\alpha \to \infty$):
		$\frac{1}{S_t\sigma \sqrt{2\pi \tau }},{\sigma }^2:=\frac{\delta }{\alpha }$
FD	$\frac{1}{\alpha }$	$\frac{1}{\alpha S_t}\frac{{(-{\omega }_{FD}{\tau }^{\gamma })}^{-\frac{1}{\alpha }}}{\Gamma (1-\frac{\gamma }{\alpha })}$
sub-BS	$\frac{1}{2}$	$\frac{1}{2S_t}\frac{\sqrt{\Gamma (1+2\gamma )}}{\Gamma (1-\frac{\gamma }{2})\sigma {\tau }^{\frac{\gamma }{2}}}$
		Non fractional regime ($\gamma \to 1$):
		$\frac{1}{S_t\sigma \sqrt{2\pi \tau }}$

, we summarize these observations, as well as the properties that are discussed for the first-order sensitivity in Section 5.

7. Concluding Remarks

In this article, we have provided a review of several subordinated market models and recalled their main properties. We have also recalled recent formulas while used for European option pricing in this context. Our main conclusions are the following:

(a)

The pricing formulas are smooth and fast converging, and provide excellent agreement with efficient numerical techniques (such as the PROJ method). Moreover, these formulas can provide useful approximations for at-the-money options, and allow for the construction of volatility curves.

(b)

We have derived several analytical formulas for risk sensitivities and shown that they also provide excellent agreement with standard numerical (Fourier) evaluations.

(c)

Thanks to these formulas, we were able to show that the presence of a time subordination in the VG, NIG, and FD models has a minimal impact on the delta hedging policy of an at-the-money option, but, on the contrary, has a direct impact on the P&L of the delta hedged portfolio.

Future work should include extending the pricing and sensitivities formulas to path-dependent instruments or to options written on several assets. It would also be interesting to determine whether these analytical results could be extended if the risk-neutral hypothesis is replaced, for instance, by approaches based on optimal quadratic hedging or utility functions.