Pricing, risk and volatility in subordinated market models

Aguilar, Jean-Philippe; Kirkby, Justin Lars; Korbel, Jan Article Pricing, risk and volatility in subordinated market models Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Aguilar, Jean-Philippe; Kirkby, Justin Lars; Korbel, Jan (2020) : Pricing, risk and volatility in subordinated market models, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 8, Iss. 4, pp. 1-27, https://doi.org/10.3390/risks8040124 This Version is available at: https://hdl.handle.net/10419/258077 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ risks Article Pricing, Risk and Volatility in Subordinated Market Models Jean-Philippe Aguilar 1,* , Justin Lars Kirkby 2and Jan Korbel 3,4,5,6 1Covéa Finance, Quantitative Research Team, 8-12 rue Boissy d’Anglas, FR75008 Paris, France 2School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA 30318, USA; [email protected] 3Section for the Science of Complex Systems, Center for Medical Statistics, Informatics, and Intelligent Systems (CeMSIIS), Medical University of Vienna, Spitalgasse 23, 1090 Vienna, Austria; [email protected] 4Complexity Science Hub Vienna, Josefstädterstrasse 39, 1080 Vienna, Austria 5 Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University, 11519 Prague, Czech Republic 6The Czech Academy of Sciences, Institute of Information Theory and Automation, Pod Vodárenskou Vˇeží 4, 182 00 Prague 8, Czech Republic *Correspondence: [email protected] Received: 26 October 2020; Accepted: 13 November 2020; Published: 17 November 2020   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. Keywords: Lévy process; subordination; option pricing; risk sensitivity; stochastic volatility; Greeks; time-change 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 Risks 2020,8, 124; doi:10.3390/risks8040124 www.mdpi.com/journal/risks Risks 2020,8, 124 2 of 27 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). Risks 2020,8, 124 3 of 27 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 2we 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 7is 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 (Ω , F , {Ft}t≥0 , P) be a probability space that is equipped with its natural filtration. Recall that a Lévy process {Xt}t≥0 is a stochastically continuous process satisfying X0= 0 ( P -almost surely), and whose increments are independent and stationary. This implies that the characteristic function Ψ(u , t):=EP[eiuXt] of a Lévy process has a semi-group structure and it admits an infinitesimal generator ψ(u), called Lévy symbol or characteristic exponent, which satisfies Ψ(u,t) = etψ(u),ψ(u):=log Ψ(u,1). (1) Risks 2020,8, 124 4 of 27 The characteristic exponent is entirely determined by the triplet (a , b , Π( d x)) , which corresponds to Lévy–Khintchine representation ψ(u) = a iu −1 2b2u2+ +∞ Z −∞eiux −1−iux 1 {|x|<1}Π(dx), (2) where a is the drift and is b the Brownian (or diffusion) component. The measure Π( d x) , assumed to be concentrated on R\{0}and satisfy +∞ Z −∞ min(1, x2)Π(dx)<∞, (3) is called the Lévy measure of the process, and determines its tail behavior and the distribution of jumps. When Π(R)<∞ , 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 Π(R) = ∞ , 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 Π(R−) = 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: Π(dx):="c+e−λ+x x1+α+ 1 {x>0}+c−e−λ−|x| |x|1+α− 1 {x<0}#dx. (4) This class contains several sub-classes, such as the tempered stable subordinators ( c−= 0) or the stable processes ( λ+=λ−= 0). When c+=c−:=C , α+=α−:=Y , λ−:=G and λ+:=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˘ i(2000). 2.2. Exponential Lévy Motions Let T> 0 and S:t∈[ 0, T]→St be the market price of some financial asset, seen as the realization of a time dependent random variable {St}t∈[0,T] on the canonical space Ω=R+ . We assume that there exists a risk-neutral measure Qunder which the instantaneous variations of Stcan be written down as: dSt St = (r−q)dt+dXt(5) where r∈R is the risk-free interest rate and q≥ 0 is the dividend yield (both being assumed to be deterministic and continuously compounded), and where {Xt}t∈[0,T] is a Lévy process. Under the dynamics (5), the terminal market price is given by ST=Ste(r−q+ω)τ+Xτ, (6) Risks 2020,8, 124 5 of 27 where τ:=T−t is the time horizon and ω is the martingale adjustment (also called convexity adjustment, or compensator), which is determined by the martingale condition EQ[ST|Ft] = e(r−q)τSt. (7) Given the form of the exponential process (6), the condition (7) is equivalent to: ω=−ψ(−i) = −logEPheX1i. (8) 2.3. Option Pricing Given a path-independent payoff function P , i.e., a positive function depending only on the terminal value ST of the market price and on some strike parameters K1 , . . . , KN> 0, then the value at time t of a contingent claim delivering a payoff Pat its maturity is equal to the following risk-neutral expectation: C=EQe−rτP(ST,K1, . . . , Kn)|Ft. (9) If the Lévy process {Xt}t∈[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: C=e−rτ +∞ Z −∞P(Ste(r−q+ω)τ+x,K1, . . . , Kn)f(x,τ)dx. (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(VG) t:=θGt+σWGt(11) where Wt is the standard Wiener process, and γ(t ,1, ν) is a Gamma process (i.e., a process whose increments γ(t+h ,1, ν)−γ(t ,1, ν) follow a Gamma distribution with mean 1 ×h and variance ν×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): ΠG(dx) = 1 ν e−1 νx x 1 {x>0}dx. (12) It follows from (12) that ΠG(R) = ∞ , 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 ν regime. The VG process is actually a Risks 2020,8, 124 6 of 27 tempered stable process itself (and, more precisely, a CGMY process), its Lévy measure admitting the following representation: ΠVG(dx) = eθx σ2 ν|x|e−rθ2 σ2+2 ν σ|x|dx. (13) Note that (13) is symmetric around the origin when θ= 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: fVG(x,t) = 2eθx σ2 νt ν√2πσΓ(t ν) x2 2σ2 ν+θ2!t 2ν−1 4 Kt ν−1 2 1 σ2r2σ2 ν+θ2|x|!(14) where Ka(X) 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: ψVG(u) = −1 νlog 1−iθνu+σ2ν 2u2, (15) allowing for a simple representation for the VG martingale adjustment: ωVG =−ψVG(−i) = 1 νlog 1−θν −σ2ν 2. (16) Remark 1. Note that, when θ= 0and ν→∞ , then ωVG → −σ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 (σ , ν ,0 )ν→0 −→ BS (σ) , is illustrated in Figure 1for decreasing ν . In particular, ν directly controls the excess kurtosis for the VG model. -1.5 -1 -0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 Figure 1. Black–Scholes (circles) as the limit of VG(σ,ν,0)ν→0 −→ BS(σ). Risks 2020,8, 124 7 of 27 3.1.2. Financial Applications As already noted, the presence of the subordination parameter ν is particularly attractive for modeling time-induced phenomena, as it allows for a non-uniform passage of time. When ν 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 ν , 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(NIG) t=βδ2It+δWIt(17) where {It}t∈[0,T] follows an Inverse Gamma distribution of shape δpα2−β2 and mean rate 1. α> 0 is a tail or steepness parameter controlling the kurtosis of the NIG distribution; the large α regime gives birth to light tails, while small α corresponds to heavier tails. β∈(−α , α− 1 ) is the skewness parameter: β< 0 (resp. β> 0) implies that the distribution is skewed to the left (resp. the right), and β= 0 that the distribution is symmetric. δ> 0 is the scale parameter and it plays an analogous role to the variance term σ2 in the Normal distribution. Let us mention that a location parameter µ∈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 {It}t∈[0,T]satisfies ΠIG(dx) = e−δ2(α2−β2) 2x x3 2 1 {x>0}dx, (18) while the Lévy measure of the NIG process itself is given by ΠNIG(dx):=αδ πeβxK1(α|x|) |x|dx. (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 Risks 2020,8, 124 8 of 27 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 fNIG(x,t):=αδt πeδt√α2−β2+β(x−µt)K1αp(δt)2+ (x−µt)2 p(δt)2+ (x−µt)2, (20) and its Lévy symbol is given by ψNIG(u) = −δqα2−(β+iu)2−qα2−β2. (21) It follows that the NIG convexity adjustment reads ωNIG =−ψNIG(−i) = δqα2−(β+1)2−qα2−β2. (22) Remark 2. When α→∞ (large steepness regime), then ωNIG → −σ2 2( 1 + 2 β) where σ2:=δ α ; when, furthermore, β= 0(symmetric process) then one recovers the usual Gaussian adjustment −σ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. Risks 2020,8, 124 15 of 27 it is immediate to see that (49) also recovers the approximation (48) in the large steepness regime, with σ2:=δ/α. Likewise, in the FD model, the ATMF price is approximated by the first term of the series (45), resulting in CFD =St α 1 Γ(1+γ α)(−ωFDτγ)1 α. (51) Taking α=2 and using (36), the ATMF price in the sub-BS model becomes Csub−BS =St 2 1 Γ(1+γ 2)pΓ(1+2γ)στ γ 2(52) which recovers (48) when γ→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, ν directly controls the tail-heaviness, as illustrated in Figure 4. In particular, large values of ν lead to steep implied volatility smiles. The ATMF prices that are obtained in Section 4.1 are helpful to approximate the implied volatility σI of the subordinated models, when St is close to F . Denoting by Ct the market price of an ATMF European call option at time tand inverting (46), we immediately see that the VG implied volatility is σVG =r2π ν Γ(τ ν) Γ(1 2+τ ν) Ct St, (53) and, similarly, inverting Equation (52), the sub-BS implied volatility is σsub−BS =2pΓ(1+2γ)Γ(1+γ 2) τγ 2 Ct St. (54) As expected, VG and sub-BS both implied volatilities recover the BS implied volatility in their limiting regimes (ν→0 and γ→1): σBS =r2π τ Ct St. (55) In the NIG case, things are a bit more complicated, because one has to solve Stδτeαδτ πK0(αδτ) = Ct(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 ρ∈R,      a0(ρ) = 1 ak(ρ) = (4ρ2−12)(4ρ2−32). . . (4ρ2−(2k−1)2) k!8k,k≥1, (57) Risks 2020,8, 124 16 of 27 then, for large zand fixed ρ, we have: Kρ(z) = z→∞rπ 2ze−z ∞ ∑ k=0 ak(ρ) zk. (58) 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 Figure 4. Implied volatility smiles of VG (σ , ν ,0 ) obtained by Formula (1) . Params: σ= 0.3, τ= 0.2, r=1%, q=0%, St=4000. The moneyness is determined by F:=Stexp((r−q)τ). In particular, when 4 ρ2− 1 = 0, i.e., when ρ=1 2 , all the ak(ρ) are null in definition (57) when k≥ 1, and we are left with: K1 2(z) = rπ 2ze−z(59) for all z . Using (58) up to k= 1 for z=αδτ and inserting into (56) , we are left with the quadratic equation X2−α√2πCt St X−1 8=0(X:=√αδτ), (60) whose positive solution reads X=1 2 α√2πCt St +s2πα2C2 t S2 t +1 2!. (61) Taylor expanding for large αand turning back to the δvariable, we obtain δ=2πα τ C2 t S2 t +1 4ατ +O1 α3(62) which, at first order, recovers (55) for σ2:=δ/α. We summarize these results in Table 2. Risks 2020,8, 124 17 of 27 Table 2. Volatility modelling for ATMF options in various subordinated models, and their limiting cases. ATMF Implied Volatility (European Options) Exponential VG σVG =q2π ν Γ(τ ν) Γ(1 2+τ ν) Ct St Low variance regime (ν→0): σVG →q2π τCt St Exponential NIG Solve Stδτeαδτ πK0(αδτ) = Ct Large steepness regime (α→∞): δ=2πα τ C2 t S2 t +1 4ατ +O1 α2 At order α0: σNIG :=qδ α=q2π τCt St sub-BS σFD =2√Γ(1+2γ)Γ(1+γ 2) τγ 2 Ct St Non-fractional regime (γ→1): σFD →q2π τCt St 5. First-Order Sensitivities The sensitivity of a contingent claim C to the underlying asset, often denoted Delta or ∆ , is defined by ∆:=∂C/∂St ; by deriving Formula (1) with respect to St 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 kVG <0, ∆− VG(kVG,σν) = F 2StΓ(τ ν) ∞ ∑ n1=0 n2=1 (−1)n1 n1!"−n1Γ(−n1+n2+1 2+τν) Γ(−n1+n2 2+1)−kVG σνn1−1 σn2−1 ν +2Γ(−2n1−n2−2τν) Γ(−n1+1 2−τν)−kVG σν2n1+1+2τν (−kVG)n2−1#. (63) - (ITM sensitivity) If kVG >0, ∆+ VG(kVG,σν) = e−qτ−∆− VG(kVG,−σν). (64) - (ATM sensitivity) If kVG =0, ∆− VG(kVG,σν) = ∆+ VG(kVG,σν) = F 2StΓ(τ ν) ∞ ∑ n=1 Γ(n 2+τν) Γ(n+1 2)σn−1 ν. (65) (ii) The Delta at time t of a European call option in the exponential NIG model is: ∆NIG =Fαeαδτ St√π ∞ ∑ n1=0 n2=1 kn1 NIG n1!Γ(−n1+n2+1 2)Kn1−n2 2+1(αδτ)δτ 2α−n1+n2 2. (66) Risks 2020,8, 124 18 of 27 (iii) The Delta at time t of a European call option in the FD model is: ∆FD =F αSt ∞ ∑ n1=0 n2=0 kn1 FD n1!Γ(1+γ−n1+n2 α)(−ωFDτγ)−n1+n2 α. (67) For illustration, in Figure 5we compare the Delta of VG (σ , ν ,0 ) while using Formula (2) with that of BS (σ) . Similarly, we compare the Dollar Gamma using Formula (3) . Figure 6provides a comparison for NIG (α ,0, δ) . 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. 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2000 3000 4000 5000 6000 7000 8000 0 1000 2000 3000 4000 5000 6000 7000 8000 Figure 5. Delta ( Left ) and Dollar Gamma ( Right ) of a call option under VG (σ , ν ,0 ) using Formula (2) and Formula (3) . Greeks of BS (σ) are provided for reference (dash lines). Params: σ= 0.3, τ= 1, r=1%, q=0%, K=4000. 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 0 500 1000 1500 2000 2500 Figure 6. Delta ( Left ) and Dollar Gamma ( Right ) of a call option under NIG (α ,0, δ) using Formulas (2) and (3) . Greeks of BS (σ= 0.3 ) are provided for reference (dash lines). Params: δ= 1.2, τ= 1, r= 1%, q=0%, K=4000. Risks 2020,8, 124 19 of 27 5.1. Delta Hedging The leading term for the Delta of the European option in the exponential VG case is given for n1=n2=1 in (63) and it reads: ∆− VG =F 2StΓ(τ ν)"Γ1 2+τν−2Γ(−3−3τν) Γ(−1 2−τν)−kVG σν−3+2τν#; (68) In the ATMF situation ( St=F , i.e., k= 0), we can re-write the martingale adjustment and Taylor expand for small σ kVG =ωVGτ=−σ2 2τ+Oσ4(69) and, recalling τν=τ ν−1 2, we obtain ∆− VG =1 2+Oσ4(τν+1). (70) It is interesting to note that, in the ATMF situation, ∆− VG ≃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 n1=0, n2=1) reads Fαeαδτ St√πK1 2(αδτ)rδτ 2α=F 2St(71) where we have used the particular value of the Bessel function of index 1 2(59) in order to simplify the expression; when F=St , we obtain ∆NIG =1 2 , which, again, turns out to be similar to the usual Black–Scholes behavior. This effect is clearly illustrated in Figures 5and 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 ( n1=n2= 0) in (67) yields ∆FD =F αSt St→F −→ 1 α(72) which explicitly depends on the tail parameter α and resumes to 1 2 in the sub-BS model, for any subordination parameter γ . Again, the subordination parameter plays no role in the delta-hedging policy of the portfolio, which is entirely governed by the tail parameter α ; in other words, it suffices to be long one unit of the underlying asset S and short α 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 ∆ 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)) EQe−rτ(ST−K)+|St=Ste−qτΠ1−Ke−rτΠ2, (73) Risks 2020,8, 124 20 of 27 where e−qτΠ1 is the option’s Delta. This quantity is known to admit a convenient representation in the Fourier space: ∆=e−qτΠ1=1 2+1 π ∞ Z0 Re "eiuk ˜ Ψ(u−i,τ) iu #du(74) where k is the log forward moneyness defined in (37) , and the “risk-neutralized” characteristic function is defined by: ˜ Ψ(u,t):=eiuωtΨ(u,t) = e(iuω+ψ(u))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 ψ(u) and the martingale adjustment ω , and by evaluating the integral by any classical recursive algorithm (such as a simple trapezoidal rule, for instance). In Table 3, 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≃ 0, which tends to accelerate the overall convergence of the series. 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=n1=n2 is the number of terms in the truncated series. Parameters: K=4000, r=1%, q=0%, τ=1. Exponential VG Model [σ=0.2,ν=0.85] Formula (2)Lewis (74) N=3N=5N=10 N=15 Deep OTM (St=3000) 2.1823 0.6347 0.0941 0.0940 0.0940 OTM (St=3500) 0.4113 0.2567 0.2455 0.2455 0.2455 ATM (St=4040.90) 0.5703 0.5718 0.5719 0.5719 0.5719 ITM (St=4500) 0.7569 0.8113 0.8134 0.8134 0.8134 Deep ITM (St=5000) 0.4729 0.8589 0.9206 0.9206 0.9206 Exponential NIG Model [α=9,δ=1.2] Formula (2)Lewis (74) N=3N=5N=10 N=15 Deep OTM (St=3000) 0.2921 0.2722 0.2747 0.2748 0.2748 OTM (St=3500) 0.4289 0.4309 0.4311 0.4311 0.4311 ATM (St=4234.09) 0.6336 0.6410 0.6412 0.6412 0.6412 ITM (St=4500) 0.6936 0.7030 0.7033 0.7033 0.7033 Deep ITM (St=5000) 0.7827 0.7966 0.7971 0.7971 0.7971 6. Second-Order Sensitivities and Portfolio Performance 6.1. Gamma, Dollar Gamma The second order derivative of a contingent claim C with respect to St is often denoted by Γ:=∂2C/∂S2 t. It is closely related to the performance, or Profit and Loss (P&L) of a portfolio: if t2>t1 are two trading days, then the P&L between t1and t2is P&L =Θ ∆t+marketP&L, ∆t:=t2−t1, (76) Risks 2020,8, 124 21 of 27 where Θis the time sensitivity of the portfolio, and the market P&L is, at order 2: marketP&L :=∆(∆St) + 1 2Γ(∆St)2∆St:=St2−St1. (77) Assuming that the portfolio has been delta-hedged, then we are left with marketP&L =$Γ∆St St2 (78) where the Dollar Gamma has been defined by $Γ:=1 2S2 tΓ ; 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 $ΓBS =St 2σ√2πτ . (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 kVG <0, Γ− VG(kVG,σν) = F 2S2 tσνΓ(τ ν) ∞ ∑ n=0 (−1)n n!"Γ(−n 2+τν) Γ(−n+1 2)−kVG σνn +2Γ(−2n−2τν) Γ(−n+1 2−τν)−kVG σν2n+2τν#. (80) - (ITM sensitivity) If kVG >0, Γ+ VG(kVG,σν) = −Γ− VG(kVG,−σν). (81) - (ATM sensitivity) If kVG =0, Γ− VG(kVG,σν) = Γ+ VG(kVG,σν) = F 2√πS2 tσνΓ(τ ν) Γ(τν−1 2) Γ(τ ν). (82) (ii) The Gamma at time t of a European call option in the exponential NIG model is: ΓNIG =Fαeαδτ S2 t√π ∞ ∑ n=0 kn NIG n!Γ(−n+1 2)Kn 2+1(αδτ)δτ 2α−n 2. (83) (iii) The Gamma at time t of a European call option in the FD model is: ΓFD =F αS2 t ∞ ∑ n=0 kn FD n!Γ(1−γ α(n+1)) (−ωFDτγ)−n+1 α. (84) Risks 2020,8, 124 22 of 27 Proof. The formulas are all straightforward to obtain, by deriving the series in Formula (2) with respect to Stand making an appropriate change of variables. For instance, in the NIG case, we have ΓNIG =Fαeαδτ S2 t√π   − ∞ ∑ n1=0 n2=1 kn1 NIG n1!Γ(−n1+n2+1 2)Kn1−n2 2+1(αδτ)δτ 2α−n1+n2 2 + ∞ ∑ n1=1 n2=1 kn1−1 NIG (n1−1)!Γ(−n1+n2+1 2)Kn1−n2 2+1(αδτ)δτ 2α−n1+n2 2   . (85) Performing the change of variables ˜ n1:=n1+ 1, ˜ n2→n2+ 1 in the second sum shows that only the terms for ˜ n2=0 survive; renaming ˜ n1:=nyields 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 F 2S2 tσνΓ(τ ν)"Γ(τν) √π+2Γ(−2τν) Γ(1 2−τν)−kVG σν2τν#. (86) Taylor expanding the VG martingale adjustment for small ν and assuming that we are not far from the money forward (St→F), we have kVG ∼ ν→0k−σ2 2τ∼ St→F−σ2 2τ, (87) therefore, the Gamma writes, at first order: Γ− VG =1 2√πStσν Γ(τ ν−1 2) Γ(τ ν)(88) and the Dollar Gamma immediately follows: $Γ− VG =St 4√πσν Γ(τ ν−1 2) Γ(τ ν). (89) While using the functional relation Γ(z+1) = zΓ(z)and the Stirling approximation (47), we have: Γ(τ ν−1 2) Γ(τ ν)=1 τ ν−1 2 Γ(τ ν+1 2) Γ(τ ν)∼ ν→0rν τ(90) and, therefore, we obtain the behavior of (89) in the low variance regime: $Γ− VG ν→0 −→ St 2σ√2πτ , (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 ν ; this is no longer the case with Risks 2020,8, 124 23 of 27 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 St→Freads ΓNIG =αeαδτ πSt K1(αδτ)(92) and therefore the Dollar Gamma is $ΓNIG =αSteαδτ 2πK1(αδτ). (93) Using the asymptotic behavior for large argument (58) for the Bessel function, we know that K1(αδτ)∼ α→∞rπ 2αδτe−αδτ (94) and, therefore $ΓNIG α→∞ −→ St 2σ√2πτ ,σ2:=δ α, (95) recovering the Black–Scholes Dollar Gamma (79) . Last in the FD model, the leading term in the series (84) for St→Fis ΓFD =1 αSt (−ωFDτγ)−1 α Γ(1−γ α), (96) and, in the sub-BS model (α=2), using the approximation (36) for the martingale adjustment, Γsub−BS =1 2StpΓ(1+2γ) Γ(1−γ 2)στ γ 2. (97) Therefore, the Dollar Gamma in the sub-BS model is $Γsub−BS =St 4pΓ(1+2γ) Γ(1−γ 2)στ γ 2. (98) and, in the non fractional limit (γ→1), we have, again, $Γsub−BS γ→1 −→ St 2σ√2πτ . (99) In Table 4, we summarize these observations, as well as the properties that are discussed for the first-order sensitivity in Section 5. Risks 2020,8, 124 24 of 27 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 (∆) 2nd Order (Γ) Exponential VG 1 21 2√πStσν Γ(τ ν−1 2) Γ(τ ν) Low variance regime (ν→0): 1 Stσ√2πτ Exponential NIG 1 2αeαδτ πStK1(αδτ) Large steepness regime (α→∞): 1 Stσ√2πτ ,σ2:=δ α FD 1 α1 αSt (−ωFD τγ)−1 α Γ(1−γ α) sub-BS 1 21 2St √Γ(1+2γ) Γ(1−γ 2)στ γ 2 Non fractional regime (γ→1): 1 Stσ√2πτ 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. Author Contributions: J.-P.A. and J.L.K. performed conceptualization and calculations and prepared the original manuscript draft. J.-P.A., J.L.K. and J.K. reviewed and finalized the manuscript. All authors have read and agreed to the published version of the manuscript. Funding: J. K. acknowledges support from Grant Agency of the Czech Republic, grant No. 19-16066S, and grant No. 20-17295S. Conflicts of Interest: The authors declare no conflict of interest. References Abramowitz, Milton, and Irene Stegun. 1972. Handbook of Mathematical Functions. Mineola: Dover Publications. Aguilar, Jean-Philippe. 2020a. Some pricing tools for the Variance Gamma model. International Journal of Theoretical and Applied Finance 23: 2050025.