Explicit option valuation in the exponential NIG model
aa r X i v : . [ q -f i n . P R ] J u l Explicit option valuation in the exponential NIG model
Jean-Philippe Aguilar ∗ July 08, 2020
Abstract
We provide closed-form pricing formulas for a wide variety of path-independent options, in theexponential Lévy model driven by the Normal inverse Gaussian process. The results are obtained inboth the symmetric and asymmetric model, and take the form of simple and quickly convergent series,under some condition involving the log-forward moneyness and the maturity of instruments. Proofs arebased on a factorized representation in the Mellin space for the price of an arbitrary path-independentpayoff, and on tools from complex analysis. The validity of the results is assessed thanks to severalcomparisons with standard numerical methods (Fourier-related inversion, Monte-Carlo simulations) forrealistic sets of parameters. Precise bounds for the convergence speed and the truncation error are alsoprovided.
Keywords:
Lévy process; Normal inverse Gaussian process; Stochastic volatility; Option pricing.
AMS subject classifications (MSC 2020):
JEL Classifications:
C00, C02, G10, G13. ∗ BRED Banque Populaire, 18 Quai de la Râpée, FR-75012 Paris, Email: [email protected]
Introduction
Whether the dramatic COVID-19 events and the subsequent turmoils in global markets were unpredictable"black swan" events in the sense of Taleb (2010) or, on the contrary, could have been forecasted (or atleast, following the terminology of Giannone & al. (2008), "nowcasted") will undoubtedly be the matterof intense debates. But what is already certain is that they demonstrate, yet again, that the kurtosis inthe distribution of asset returns far exceeds the tails of the Normal one, and that market volatility is notconstant over time; it should therefore be a minimal requirement for any reliable market model that theyinclude (at least) these two stylized facts.It has now long been known that exponential - sometimes also called geometrical - Lévy models fulfilthese conditions. Such models have been introduced in quantitative finance during the late 1990s / early2000s in several influential works, and make the assumption that asset log returns are driven by somedrifted Lévy process: a Normal inverse Gaussian (NIG) process in Barndorff-Nielsen (1995, 1997), aVariance Gamma (VG) process in Madan et al. (1998), a hyperbolic or a generalized hyperbolic processin Eberlein and Keller (1995); Eberlein (2001), a CGMY process in Carr & al. (2002) or a stable (or α -stable) process in Mittnik & Rachev (2000); Carr and Wu (2003). Stable distributions, in particular, maybe noted for their historical importance, having been considered a credible candidate for the modelling ofasset prices as early as in the 1960s by Mandelbrot (1963) in the context of the cotton market, thus pavingthe way to the more generic setup of exponential Lévy models. Readers who may be less familiar withthe broad family of Lévy processes and their applications to finance are invited to refer to the classicalreferences Bertoin (1996); Schoutens (20003); Cont & Tankov (2004); Rachev et al. (2011).In the present work, we will be particularly interested in the class of exponential Lévy models whoseLévy process is distributed according to a NIG distribution, namely, the class of exponential NIG models.NIG distributions were originally introduced for physical purpose, more precisely to model the complexbehavior of dunes and beach sands, in the seminal article Barndorff-Nielsen (1977); as noted above, theyhave subsequently been introduced for financial purpose approximately two decades later, because theyfeature several degrees of freedom that have a direct empirical interpretation in terms of financial timeseries. First, they possess fat tails, allowing for the presence of extreme variations of prices (positive ornegative jumps); when the tail parameter goes to infinity, then the NIG distribution degenerates into theNormal distribution and the exponential NIG model recovers the Black-Scholes model (Black & Scholes(1973)). Second, NIG distributions can be skewed, allowing to capture the asymmetry that can be observed1n the distribution of jumps (price drops occurring more often than raises). Last, but not least, a NIGprocess can be interpreted as a drifted Brownian motion whose time follows an inverse Gamma process- this is a consequence of the fact that the NIG distribution is actually a particular case of a so-calledNormal variance-mean mixture, the mixing distribution being the inverse Gaussian (IG) distribution; theNIG process is therefore a time changed Lévy process, which allows for stochastic volatility modelling andrelated phenomena, such as clustering or negative correlation between the returns and their volatility (seedetails in Carr and Wu (2004)). Let us also mention that, as observed by Mechkov (2015) the NIG processis also a limiting case of the Heston model (Heston (1993)) in the fast reversion limit, which makes it anideal candidate for calibrating the short time part of a Heston volatility surface.Of course, since it was introduced, the exponential NIG model has been proved to provide a verygood fitting to financial data many times. Let us mention, among others, initial tests for daily returnson Danish and German markets in Barndorff-Nielsen (1995); Rydberg (1997) and subsequently on theFTSE All-share index (also known as "Actuaries index") in Venter & de Jongh (2002). More recently, theimpact of high frequency trading has also been taken into account, and calibrations have been performedon intraday returns e.g. in Figueroa-López et al. (2012) for different sampling frequencies. Let us alsomention that multivariate extensions of the exponential NIG model, i.e., featuring a different time changefor different assets, have also been considered (see Luciano & Semeraro (2010) and references therein).As one could expect, pricing contingent claims turns out to be a tougher task in the exponentialNIG model than it is in the usual Black-Scholes framework. Numerical methods are largely favored,including Monte-Carlo valuation methods (Ribeiro and Webber (2003)), numerical evaluation of Fourier(Lewis (2001)) and Fast Fourier (Carr & Madan (1999)) transforms. The success of Fourier transformmethods is strongly linked to the relative simplicity of the characteristic function of most exponentialLévy models, and has opened the way to a wide range of other transform based approaches: they include,among others, the COS method by Fang & Osterlee (2008), the Hilbert transform method (see notablya recent application to time-changed Lévy processes in Zeng and Kwok (2014)) or the local basis FramePROJection (PROJ) method by Kirkby (2015). Efforts have also been made towards analytic evaluation orapproximations: in Ivanov (2013), a closed-form formula (in terms of Appel functions) for the Europeancall is derived in the particular case where the NIG distribution has a tail parameter of / , and inAlbrecher & Predota (2004) approximations and bounds are provided for Asian options.In this paper, we would like to show that it is actually possible to obtain tractable closed-form pricingformulas in the exponential NIG model, for a broad range of path independent instruments. This is2ade possible by a remarkable property allowing to express the Mellin transform of an arbitrary pathindependent option as the product of the Mellin transforms of its payoff and of the NIG probabilitydensity. Inverting it by means of residue summation yields the option price, computed under the form ofquick convergent residue series whose terms are directly expressed in terms of the model’s parameters. ThisMellin residue summation method has been used very recently within the framework of other exponentialLévy models, namely in the Finite Moment Log Stable (FMLS) model in Aguilar & Korbel (2019) and inthe exponential VG model in Aguilar (2020); in the present paper, we will therefore demonstrate that thetechnique is also well-suited to the exponential NIG model. Moreover, we will establish pricing formulasfor both the symmetric and the asymmetric NIG processes, while the formulas in the VG case in Aguilar(2020) were mainly obtained for the symmetric VG process. Due to the nature of the residues series,however, we will need to introduce a restriction on the model parameters to ensure the convergence to theprice. We will show that this condition is compliant with most of the implied parameters calibrated in theliterature; moreover, when options are not far from the money, it is automatically satisfied.The paper is organized as follows: in section 2, we start by recalling fundamental concepts on the NIGprocess and its implementation via exponential Lévy models. In section 3, we focus on the symmetricNIG process: after establishing the pricing formula in the Mellin space for an arbitrary path independentinstrument, we evaluate, analytically, the price of the European and digital options, as well as payoffsfeaturing more exotic attributes (power options, log contracts, . . . ). In section 4 we extend the pricingformula to the more general case of the asymmetric model, and provide analytic formulas for the digital andEuropean prices. In section 5, precise bounds for the convergence speed and the truncation errors of theseries are obtained, and the validity of the results is assessed by comparing them with classic numericalmethods (Fourier inversion, Monte Carlo simulations). For the reader’s convenience, the paper is alsoequipped with two appendices: in appendix A we provide a short overview of the Mellin transform, andin appendix B we recall some important special function identities that are used throughout the paper. In this section we recall important concepts on NIG distributions and processes; more details can befound in the initial articles by Barndorff-Nielsen or in subsequent review articles like Hanssen & Øigård(2001); Papantolen (2008). We also introduce the exponential NIG model, following the classical setup ofexponential Lévy models such as defined e.g. in Schoutens (20003); Tankov (2010).3 .1 The Normal inverse Gaussian process
The Normal inverse Gaussian (NIG) process can be defined by in several different ways: classically, it isdefined either as a process whose increment follow a NIG distribution, in terms of its Lévy measure, or asa time-changed Lévy process. Let us also mention that, as remarked in Mechkov (2015), the NIG processcan also be seen as the limit of a Fast Reverting Heston (FRH) process.
NIG density
The NIG distribution, denoted by
NIG( α, β, δ, µ ) , is a four-parameter distribution whosedensity function is: f ( x ) := αδπ e δ √ α − β + β ( x − µ ) K (cid:16) α p δ + ( x − µ ) (cid:17)p δ + ( x − µ ) . (1)The function z → K ( z ) is the modified Bessel function of the second kind and of index 1 (sometimesalso called Macdonald function, see definitions and properties in appendix B). α > is a tail or steepness parameter controlling the kurtosis of the distribution; the large α regime gives birth to light tails, whilesmall α corresponds to heavier tails. β ∈ ( − α, α − is the skewness parameter: β < (resp. β > )implies that the distribution is skewed to the left (resp. the right), and β = 0 that the distribution issymmetric around the location parameter µ ∈ R . δ > is the scale parameter and plays an analogue roleto the variance term σ in the Normal distribution; when β = 0 , the Normal distribution is itself recoveredin the large steepness regime: NIG( α, , δ, µ ) −→ α →∞ N ( µ, σ ) , σ := δα . (2)We say that a stochastic process { X t } t ≥ is a NIG process if it has NIG distributed increments, thatis if X t + h − X t ∼ NIG( α, β, δh, µh ) for all h ≥ ; it follows from (1) that the density of the processconditionally to X = 0 is (with a slight abuse of notations): f ( x, t ) := αδtπ e δt √ α − β + β ( x − µt ) K (cid:16) α p ( δt ) + ( x − µt ) (cid:17)p ( δt ) + ( x − µt ) . (3)It is also possible to define the NIG process as a time-changed drifted Brownian motion: if { I t } t ≥ is aprocess distributed according to an Inverse Gamma density of shape δ p α − β and mean rate 1 and if { W t } t ≥ is a standard Wiener process, then the process X t = βδ I t + δW I t (4)4s a centered NIG process ( µ = 0 ). The process { I t } t ≥ is a tempered stable subordinator; it has positivejumps, and therefore is interpreted as a business time that can differ from the operational time , theoccurence of jumps corresponding to periods of intense business activity. A similar interpretation holdsfor instance in the case of the Variance Gamma process, which features another example of temperedstable subordination (via a Gamma process). Lévy symbol
The NIG process is a (pure jump) Lévy process whose characteristic function Ψ( u, t ) := E [ e iuX t ] can be written down as Ψ( u, t ) = e tψ ( u ) , where the characteristic exponent, or Lévy symbol, isknown in exact form: ψ ( u ) := log Ψ( u,
1) = iµu − δ (cid:16)p α − ( β + iu ) − p α − β (cid:17) . (5)The process admits the Lévy-Khintchine triplet ( a, , ν (d x )) , where the drift a and the Lévy measure ν are defined by a := µ + 2 αδπ Z sinh( βx ) K ( αx ) d xν (d x ) := αδπ e βx K ( α | x | ) | x | d x, (6)allowing to write down the characteristic exponent (5) in terms of its Lévy-Khintchine representation: ψ ( u ) = iau + Z R ( e iux − − iux {| x | < } ) ν (d x ) . (7)Let us observe that it follows from the definition of the Lévy measure ν that the NIG process has infinitevariation and infinite intensity (i.e. ν ( R ) = ∞ ), and therefore possesses a very rich dynamics with infinitenumber of jumps on any time interval - this is why no Brownian component is even needed in the Lévy-Khintchine triplet. We should also note that the NIG process has all its moments finite, which is not thecase with (double-sided) α -stable processes for instance: this is because the Bessel function admits theasymptotic behavior (see (144)) K ( | x | ) ∼ | x |→∞ r π | x | e −| x | , (8)and therefore the tails of the NIG measure ν are less heavy than the tails of the α -stable measure (whichhas polynomial decrease in / | x | α ) . In other words, the jumps in the NIG process are not as big as for α -stable processes, but allow finiteness of moments and therefore of option prices; in the α -stable case,5his would be achieved only for spectrally negative processes (i.e., having negative jumps only). Fast reverting Heston limit
Following Mechkov (2015), we choose the following formulation for theHeston dynamics: d x t = − σ z t d t + σ √ z t d W (1) t d z t = a (cid:16) (1 − z t )d t + γ √ z t d W (2) t (cid:17) (9)where the two Brownian motions have correlation ρ , and we consider the fast reversion limit a → ∞ inthe CIR process driving the stochastic multiplier z t . Then, in this limit, the process { x t } t ≥ is distributedaccording to a NIG distribution: x t ∼ NIG p − ργσ + γ σ γσ (1 − ρ ) , − γσ − ρ γσ (1 − ρ ) , σ p − ρ γ t, − σργ t ! . (10) Model specification
Let
T > and S : t ∈ [0 , T ] → S t be the market price of some financial asset,seen as the realization of a time dependent random variable { S t } t ∈ [0 ,T ] on the canonical space Ω = R + equipped with its natural filtration. We assume that there exists a risk-neutral measure Q under whichthe instantaneous variations of S t can be written down as: d S t S t = ( r − q ) d t + d X t (11)where r ≥ is the risk-free interest rate and q ≥ is the dividend yield (both assumed to be deterministicand continuously compounded), and where { X t } t ∈ [0 ,T ] is the NIG process. The solution to the stochasticdifferential equation (11) it the exponential process S T = S t e ( r − q + ω ) τ + X τ , ω := − ψ ( − i ) , (12)where τ := T − t is the time horizon and ω is the martingale adjustment (also called convexity adjust-ment, or compensator) determined by the martingale condition E Q [ S T | S t ] = e ( r − q ) τ S t ; it follows from thedefinition of the Lévy symbol (5) that this adjustment is equal to: ω = − µ + δ (cid:16)p α − ( β + 1) − p α − β (cid:17) . (13)6t is interesting to note that, in the large steepness regime, (13) has the following asymptotic behavior: ω ∼ α →∞ − µ − σ β ) , σ := δα . (14)Taking µ = 0 (centered process) and β = 0 (symmetric process), (14) recovers the the Gaussian martingaleadjustment − σ / , and the exponential NIG model (11) degenerates into the Black-Scholes model. Contingent claim valuation
Given a path-independent payoff function P , i.e., a positive functiondepending only on the terminal value S T of the market price and on some strike parameters K , . . . , K N > , then the value at time t of a contingent claim delivering a payoff P at maturity is equal to the followingrisk-neutral expectation: C = E Q (cid:2) e − rτ P ( S T , K , . . . , K n ) | S t (cid:3) . (15)The conditional expectation (15) can be achieved by integrating all possible realizations for the payoff overthe probability density of the NIG process, thus resulting in: C = e − rτ + ∞ Z −∞ P ( S t e ( r − q + ω ) τ + x , K , . . . , K n ) f ( x, τ ) d x. (16) In this section, we assume that β = 0 , i.e., that the process { X t } t ∈ [0 ,T ] in (11) is distributed according tothe symmetric distribution NIG( α, , δt, µt ) . First, we establish a general pricing formula for an arbitrarypath independent instrument; then, we apply this formula to the analytic evaluation of several optionsand contracts. Let us start by establishing a representation for the symmetric NIG density f ( x, t ) under the form of aMellin-Barnes integal. Lemma 3.1.
For any c ∈ R + , the following holds true: f ( x, t ) = α π e αδt c + i ∞ Z c − i ∞ Γ (cid:16) s (cid:17) K − s ( αδt ) (cid:18) δtα (cid:19) s | x − µt | − s d s iπ . (17)7 roof. Taking β = 0 in (3) yields: f ( x, t ) = αδtπ e αδt K (cid:16) α p ( δt ) + ( x − µt ) (cid:17)p ( δt ) + ( x − µt ) . (18)Using the Mellin transform for the Bessel function (see table 7 in appendix A with ν = 1 ) and the Mellininversion formula (127), we can write: K (cid:16) α p ( δt ) + ( x − µt ) (cid:17)p ( δt ) + ( x − µt ) = 12 δτ c + i ∞ Z c − i ∞ Γ (cid:16) s (cid:17) K − s ( αδτ ) (cid:18) δtα (cid:19) s | x − µt | − s d s iπ (19)for any c > . Inserting into (18) yields the representation (17).Let us now introduce the double-sided Mellin transform of the payoff function: P ∗ ( s ) = ∞ Z −∞ P (cid:16) S t e ( r − q + ω ) τ + x , K , . . . , K n (cid:17) | x − µτ | − s d x (20)and assume that it exists for Re ( s ) ∈ ( c − , c + ) for some real numbers c − < c + . Then, as a consequenceof the risk-neutral pricing formula (16) and of lemma 3.1, we immediately obtain: Proposition 3.2 (Factorization in the Mellin space) . Let c ∈ (˜ c − , ˜ c + ) where (˜ c − , ˜ c + ) := ( c − , c + ) ∩ R + is assumed to be nonempty. Then the value at time t of a contingent claim delivering a payoff P ( S T , K , . . . , K n ) at its maturity t = T is equal to: C = α π e ( αδ − r ) τ c + i ∞ Z c − i ∞ Γ (cid:16) s (cid:17) P ∗ ( s ) K − s ( αδτ ) (cid:18) δτα (cid:19) s d s iπ . (21)Throughout the paper, our purpose will be to express the complex integral (21) as a sum of residuesassociated to the singularities of the integrand. Schematically, we will therefore be able to express theprice of a contingent claim under the form of a series: α π e ( αδ − r ) τ × X (cid:20) residues of Γ (cid:16) s (cid:17) P ∗ ( s ) × particular values of K − s ( αδτ ) × powers of δτα (cid:21) . (22)As we will see, the residues turn out to have to be computed in the multidimensional sense, because,depending on the payoff’s complexity, the evaluation of P ∗ ( s ) can call for the introduction of a second8ellin variable s (in the asymmetric case, we will see that one even needs a third Mellin variable s ).However, as only Gamma functions are involved, these residues are straightforward to compute, even inthe C n sense.Before proceeding to pricing itself, let us introduce the notation for the forward strike F and the logforward moneyness k : F := Ke − ( r − q ) τ , k := log S t F + ωτ = log S t K + ( r − q + ω ) τ. (23)It will also be useful to introduce k := k + µτ ; taking β = 0 in the definition of the martingale adjustment(13), we have: k = log S t K + (cid:16) r − q + δ (cid:16)p α − − α (cid:17)(cid:17) τ. (24)Note that k is independent of the location µ (in both the symmetric and asymmetric cases). Last, weneed to introduce a restriction on the parameters, that will be fundamental for the series to converge: Assumption 1.
In all of the following, and unless otherwise stated, we will assume that the model’s inputsare such that | k | δτ < . (25) We start our applications of proposition 3.2 with the determination of the price of the digital (also calledbinary) options, and of the vanilla European option.
Digital option (asset-or-nothing)
The asset-or-nothing call option consists in receiving a unit of theunderlying asset S T , on the condition that it exceeds a predetermined strike price K . The payoff cantherefore be written down as: P a/n ( S T , K ) := S T { S T >K } . (26) Formula 1 (Asset-or-nothing call) . The value at time t of an asset-or-nothing call option is: C a/n = Kαe ( αδ − r ) τ √ π ∞ X n =0 n =0 k n n !Γ(1 + − n + n ) K n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n . (27)9 roof. Step 1: Let us first assume that k < . We remark that, using notations (23), we can write P a/n (cid:16) S t e ( r − q + ω ) τ + x , K (cid:17) = K e k + x { x> − k } . (28)Using a Mellin-Barnes representation for the exponential term (see table 7 in appendix A): e k + x = c + i ∞ Z c − i ∞ ( − − s Γ( s )( k + x ) − s d s iπ ( c > (29)and inserting into (20), we get: P ∗ ( s ) = K c + i ∞ Z c − i ∞ ( − − s Γ( s ) ∞ Z − k ( k + x ) − s ( x − µτ ) − s d x d s iπ (30) = K c + i ∞ Z c − i ∞ ( − − s Γ( s )Γ(1 − s )Γ( s + s − s ) ( − k ) − s − s +1 d s iπ (31)where the x -integral exists because − ( k + µτ ) = − k > by hypothesis. Using proposition 3.2 and theLegendre duplication formula (136), we obtain the price of the asset-or-nothing call: C a/n = Kαe ( αδ − r ) τ √ π c + i ∞ Z c − i ∞ c + i ∞ Z c − i ∞ ( − − s Γ( s )Γ(1 − s )Γ( s + s − s +12 ) ( − k ) − s − s +1 K − s ( αδτ ) × (cid:18) δτ α (cid:19) s d s iπ d s iπ (32)which converges in the subset { ( s , s ) ∈ C , < Re ( s ) < , Re ( s + s ) > } and can be analyticallycontinued outside this polyhedron, except when the Gamma functions in the numerator are singular, thatis, when their arguments equal a negative integer. If we consider the singularities induced by Γ( s ) at s = − n , n ∈ N and by Γ( s + s − at s + s − − n , n ∈ N , then, the associated residuesare straightforward to compute via the change of variables u := s + s − , v := s , and via the singularbehavior (126) for the Gamma functions; they read: Kαe ( αδ − r ) τ √ π ( − n ( − n n ! ( − n n ! Γ(1 + n )Γ(1 + − n + n ) ( − k ) n K n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n . (33)Simplifying and summing all residues (33) yields the announced series (27).10tep 2: Let us now assume that k > : in that case, the x -integral on the interval ( − k , ∞ ) in (30) doesnot converge. But, as { X t } t ∈ [0 ,T ] is a Q -martingale, we can write: E Q [ S T { S T >K } | S t ] = S t e ( r − q ) τ − E Q [ S T { S T We remark that, using notations (23), we can write: P eur ( Se ( r − q + ω ) τ + x , K ) = K ( e k + x − { x> − k } . (41)Then, we use the Mellin-Barnes representation (see table 7 in appendix A): e k + x − c + i ∞ Z c − i ∞ ( − − s Γ( s )( k + x ) − s d s iπ ( − < c < (42)and we proceed exactly the same way than for proving Formula 1; note that the n -summation in (40)now starts in n = 1 instead of n = 0 , because the strip of convergence of (42) is reduced to < − , > instead of < , ∞ > in (29).Let us examine the series (40) in the large steepness regime ( α → ∞ ). It follows from the asymptoticbehavior of the Bessel function for large arguments (144) that: K n − n ( αδτ ) ∼ α →∞ √ π √ αδτ e − αδτ , (43)and from (14) that: k ∼ α →∞ log S t K + (cid:18) r − q − δ α (cid:19) τ. (44)Therefore, denoting σ := δα , we obtain C ( α →∞ ) eur = Ke − rτ ∞ X n =0 n =1 n !Γ(1 + − n + n ) (cid:18) log S t K + (cid:18) r − q − σ (cid:19) τ (cid:19) n (cid:18) σ τ (cid:19) − n n (45)12hich is the series expansion of the Black-Scholes formula for the European call that was derived in Aguilar(2019). Digital option (cash-or-nothing) The payoff of the cash-or-nothing call option is P c/n ( S T , K ) = { S T >K } (46)and therefore the option price itself is: C c/n = 1 K (cid:0) C a/n − C eur ) (cid:1) . (47)Using formulas 1 and 2, it is immediate to see that: Formula 3 (Cash-or-nothing call) . The value at time t of a cash-or-nothing call option is: C c/n = αe ( αδ − r ) τ √ π ∞ X n =0 k n n !Γ(1 − n ) K n +12 ( αδτ ) (cid:18) δτ α (cid:19) − n +12 . (48)In (48), only terms for n = 0 and n = 2 p + 1 , p ∈ N actually survive (because of the divergence ofthe Gamma function in the denominator when n = 2 p , p ≥ ). Therefore, using the particular values ofthe Gamma function at negative half-integers (134) and of the Bessel function for ν = (145), we canre-write formula 3 as: C c/n = e − rτ 12 + απ e αδτ ∞ X p =0 ( − p k p +10 p !(2 p + 1) K p +1 ( αδτ ) (cid:18) δτα (cid:19) − p . (49)The representation (49) is less compact than formula 3, however it allows for a direct computation of theput option: indeed, using E Q [ { S T >K } | S t ] = 1 − E Q [ { S T Gap option A gap (sometimes called pay-later) call has the following payoff: P gap ( S T , K , K ) = ( S T − K ) { S T >K } (52)and degenerates into the European call when trigger and strike prices coincide ( K = K = K ). From thedefinition (52), it is immediate to see that the value at time t of the Gap call is: C gap = C a/n − K C c/n (53)where the value of the asset-or-nothing and cash-or-nothing calls are given by formulas 1 and 3 for K = K . Power options Power options deliver a non linear payoff and are an easy way to increase the leverageratio of trading strategies; the payoffs of the digital power calls are P pow.c/n ( S T , K ) = { S aT >K } P pow.a/n ( S T , K ) = S aT { S aT >K } (54)for some a > , and the power European call is: P pow.eur ( S T , K ) := [ S aT − K ] + . (55)Introducing the notation k a := log S t K a + ( r − q + ω ) τ , k ,a := k a + µτ (56)then we can remark that: P pow.a/n ( S t e ( r − q + ω ) τ + x , K ) = Ke a ( k a + x ) { x> − k a } . (57)14herefore, using the representations (see table 7 in appendix A) e a ( k a + x ) = c + i ∞ Z c − i ∞ ( − − s a − s Γ( s )( k a + x ) − s d s iπ ( c > (58)and e a ( k a + x ) − c + i ∞ Z c − i ∞ ( − − s a − s Γ( s )( k a + x ) − s d s iπ ( − < c < (59)and proceeding exactly the same way than for proving formulas 1, 2 and 3, we obtain: Formula 4 (Power options) . The values at time t of the power options are:- Asset-or-nothing power call: C pow.a/n = Kαe ( αδ − r ) τ √ π ∞ X n =0 n =0 a n k n ,a n !Γ(1 + − n + n ) K n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n ; (60) - European power call: C pow.eur = Kαe ( αδ − r ) τ √ π ∞ X n =0 n =1 a n k n ,a n !Γ(1 + − n + n ) K n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n ; (61) - Cash-or-nothing power call: C pow.c/n = αe ( αδ − r ) τ √ π ∞ X n =0 k n ,a n !Γ(1 − n ) K n +12 ( αδτ ) (cid:18) δτ α (cid:19) − n +12 . (62)It is clear that, for the series (60), (61) and (62) to converge, assumption 1 has to be satisfied by k ,a and no longer by k , that is: (cid:12)(cid:12)(cid:12)(cid:12) k ,a δτ (cid:12)(cid:12)(cid:12)(cid:12) < . (63) Log options, log contract Log options are, basically, options on the rate of return of the underlying(Wilmott (2006)). The payoff of a log call and of a log put are: P log call ( S T , K ) := [log S T − log K ] + , P log put ( S T , K ) := [log K − log S T ] + . (64)15he log contract, introduced by Neuberger (1994), is a forward contract that is obtained by being long ofa log call and short of a log put, resulting in P log contract ( S T , K ) = log S T K . (65)Note that a delta-hedged log contract with K = 1 is actually a synthetic variance swap: indeed, bydenoting the quadratic variation of S by < S > and using Itô’s lemma, it is well known that, in theBlack-Scholes model, E Q [ < S > T − < S > t | S t ] = 2 E Q (cid:20) − log S T S t + S T S t − | S t (cid:21) . (66)In the more general framework of exponential Lévy models, the overall multipliers in the r.h.s. of (66)are different from 2 and have been determined in Carr & Wu (2012); for instance in the symmetric NIGmodels, it is equal to α ( α −√ α − which, as expected, tends to when α → ∞ . Let us therefore show howto derive pricing formulas for the log options and the log contract in this model: remarking that, usingnotations (23), P log call ( S t e ( r − q + ω ) τ + x , K ) = [ k + x ] + , (67)it follows that the Mellin transform for the payoff function (20) reads, for the log call: P ∗ ( s ) = ∞ Z − k ( k + x ) ( x − µτ ) − s d x = ( − k ) − s ( s − s − (68)and, using proposition 3.2, that the log call price itself writes: C log = αe ( αδ − r ) τ π c + i ∞ Z c − i ∞ Γ( s )( s − s − 1) ( − k ) − s K − s ( αδτ ) (cid:18) δτα (cid:19) s d s iπ (69)where c > . Similarly, the log put writes: P log = αe ( αδ − r ) τ π c + i ∞ Z c − i ∞ Γ( s )( s − s − k − s K − s ( αδτ ) (cid:18) δτα (cid:19) s d s iπ (70)where c > . Summing all residues arising at s = 2 , s = 1 and s = − n , n ∈ N , grouping the termsand simplifying yields: 16 ormula 5 (Log options, log contract) . The value at time t of a log option is:- Log call: C log = e − rτ " k αe αδτ π ∞ X n =0 ( − n − k n n !(2 n − K n ( αδτ ) (cid:18) δτα (cid:19) − n +1 ; (71) - Log put: P log = e − rτ " − k αe αδτ π ∞ X n =0 ( − n − k n n !(2 n − K n ( αδτ ) (cid:18) δτα (cid:19) − n +1 ; (72) - Log contract: C log − P log = e − rτ k . (73)Recall that, when α → ∞ , ω ∼ − µ − σ , where σ := δα and therefore the log contract (73) becomes ( C log − P log ) ( α →∞ ) = e − rτ (cid:18) log S t K + ( r − q − σ τ (cid:19) (74)which, taking K = 1 , is the formula originally obtained by Neuberger (1994) for the price of a log contractin the Black-Scholes model. Capped payoffs Suppose that we wish introduce a cap to limit the exercise range of a digital optionfor example; in this case, the payoff of the cash-or-nothing call would read: P capped c/n ( S T , K ) := { K − For any c ∈ R , the following holds true: f ( x, t ) = α π e γδt × Z c + i R ( − − s β − s Γ (cid:16) s (cid:17) Γ( s ) K − s ( αδt ) (cid:18) δtα (cid:19) s | x − µt | − s ( x − µt ) − s d s d s (2 iπ ) . (81) Proof. Like in the proof of lemma 3.1, we introduce the Mellin representation (19) for the Bessel functionthat holds for c ∈ R , and we introduce a supplementary representation for the exponential term (see18able 7 in appendix A): e β ( x − µt ) = c + i ∞ Z c − i ∞ ( − − s β s Γ( s ) ( x − µt ) − s d s iπ (82)that holds for c ∈ R + . Inserting (19) and (82) into the density (3) yields the reprensentation (81).Let us now introduce the asymmetric analogue to the P ∗ ( s ) function (20): P ∗ ( s , s ) = ∞ Z −∞ P (cid:16) S t e ( r − q + ω ) τ + x , K , . . . , K n (cid:17) | x − µτ | − s ( x − µτ ) − s d x (83)and assume that it exists for ( Re ( s ) , Re ( s )) ∈ P for a certain subset P ⊂ R . Then, as a consequenceof the risk-neutral pricing formula (16) and of lemma 4.1, we immediately obtain: Proposition 4.2 (Factorization in the Mellin space) . Let c ∈ ˜ P where ˜ P := P ∩ R is assumed to benonempty. Then the value at time t of a contingent claim delivering a payoff P ( S T , K , . . . , K n ) at itsmaturity t = T is equal to: C = α π e ( γδ − r ) τ Z c + i R ( − − s β − s Γ (cid:16) s (cid:17) Γ( s ) P ∗ ( s , s ) K − s ( αδτ ) (cid:18) δτα (cid:19) s d s d s (2 iπ ) . (84) To illustrate some applications of proposition 4.2, we compute the price of the digital and Europeanoptions, whose payoffs were defined in subsection 3.2. We also recall the notation for the Pochhammersymbol ( a ) n := Γ( a + n )Γ( a ) . Formula 6 (Asset-or-nothing call) . The value at time t of an asset-or-nothing call option is: C a/n = Kαe ( γδ − r ) τ √ π ∞ X n ,n ,n =0 ( − n + n + 1) n k n β n n ! n !Γ(1 + − n + n + n ) K n − n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n n . (85) Proof. Step 1: The proof starts like the proof of formula 1, by assuming k < , by remarking that P a/n (cid:16) S t e ( r − q + ω ) τ + x , K (cid:17) = K e k + x { x> − k } . (86)19nd by introducing a Mellin-Barnes representation for the exponential term in the option’s payoff: e k + x = c + i ∞ Z c − i ∞ ( − − s Γ( s )( k + x ) − s d s iπ ( c > . (87)Therefore, the P ∗ ( s , s ) function (83) reads: P ∗ ( s , s ) = K c + i ∞ Z c − i ∞ ( − − s Γ( s ) ∞ Z − k ( k + x ) − s ( x − µτ ) − s − s d x d s iπ (88) = K c + i ∞ Z c − i ∞ ( − − s Γ( s )Γ(1 − s )Γ( s + s + s − s + s ) ( − k ) − s − s − s +1 d s iπ (89)where the x -integral exists because k < . Using proposition 4.2, we obtain the price of the asset-or-nothing call: C = Kα π e ( γδ − r ) τ Z c + i R ( − − s − s β − s Γ( s )Γ( s )Γ( s )Γ(1 − s )Γ( s + s + s − s + s ) ( − k ) − s − s − s +1 × K − s ( αδτ ) (cid:18) δτα (cid:19) s d s d s d s (2 iπ ) (90)which converges in the subset { ( s , s , s ) ∈ C , Re ( s ) > , Re ( s ) > , < Res ( s ) < , Re ( s + s + s ) > } and can be analytically continued outside this polyhedron, except when the Gamma functions in thenumerator are singular. If we consider the singularities induced by Γ( s ) at s = − n , n ∈ N , by Γ( s ) at s = − n , n ∈ N and by Γ( s + s + s − at s + s + s − − n , n ∈ N , then, the associatedresidues are straightforward to compute via the change of variables u := s + s + s − − , v := s , w = s and via the singular behavior (126) for the Gamma functions; they read: Kαe ( γδ − r ) τ π ( − n + n β n ( − n n ! ( − n n ! ( − n n ! Γ(1 + n )Γ( − n + n − n +12 )Γ( − n + n + 1) ( − k ) n × K − − n n n ( αδτ ) (cid:18) δτ α (cid:19) − n n n . (91)Using the Legendre duplication formula (136) and the definition of the Pochhammer symbol (137), wewrite: Γ( − n + n + n +12 )Γ( − n + n + 1) = √ π − n + n + n ( − n + n + 1) n Γ(1 + − n + n + n ) . (92)20nserting into (91), simplifying and summing all residues for n , n , n ∈ N yields the series (85).Step 2: Like in the proof of formula 1, extension to the case k > is performed thanks to the parity E Q [ S T { S T >K } | S t ] = S t e ( r − q ) τ − E Q [ S T { S T Like in the proof of formula 2, we remark that we can write: P eur ( Se ( r − q + ω ) τ + x , K ) = K ( e k + x − { x> − k } . (95)Then, we use the Mellin-Barnes representation (see table 7 in appendix A): e k + x − c + i ∞ Z c − i ∞ ( − − s Γ( s )( k + x ) − s d s iπ ( − < c < (96)and we proceed exactly the same way than for proving Formula 6; the n -summation in (94) starts in n = 1 instead of n = 0 , because the strip of convergence of (96) is reduced to < − , > instead of < , ∞ > in (87).By difference of (85) and (94), we immediately obtain the formula for the cash-or-nothing call: Formula 8 (Cash-or-nothing call) . The value at time t of a cash-or-nothing call option is: C c/n = αe ( γδ − r ) τ √ π ∞ X n ,n =0 ( − n + 1) n k n β n n ! n !Γ(1 + − n + n ) K n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n . (97) In this section, we start by determining what restriction is induced by assumption 1 in terms of accessibleoption maturities, and we provide some precise estimates for the convergence speed and the truncation21rrors of the series. Then, we compare the various pricing formulas established in the above with severalnumerical tools, and demonstrate the reliability and efficiency of the results. We start by remarking that the at-the-money (ATM) situation ( S t = K ) is a favorable situation forsatisfying assumption 1. Indeed, in that case, we have: | k | δτ = (cid:12)(cid:12)(cid:12)(cid:12) r − qδ + p α − ( β + 1) − p α − β (cid:12)(cid:12)(cid:12)(cid:12) . (98)In the symmetric model in particular, it is clear that − r − qδ < r − qδ + p α − − α < r − qδ (99)and therefore assumption 1 is satisfied as soon as r − q < δ ; according to the implied parameters in table1, the smallest calibrated value for δ is . , therefore assumption 1 is satisfied (independently of α andof other market parameters) as soon as the risk-free interest rate is smaller than , which is of coursethe case for most financial applications.In the more general non at-the-money and non symmetric case, satisfying assumption 1 necessitatessome restriction on the option’s maturities, depending on the moneyness situation. Assuming that µ = 0 (as option prices are not sensitive to µ ) and, introducing ρ ± := log S t K ± δ − r + q − ω , (100)then it is not hard to see that:- If S t > K (in-the-money (ITM) situation), then assumption 1 is satisfied if τ > ρ + or τ < ρ − ;- If S t < K (out-of-the-money (OTM) situation), then assumption 1 is satisfied if τ > ρ − or τ < ρ + .In table 1, we illustrate this rule on several implied NIG parameters, calibrated in the literature on variousoption markets: OBX options in Saebø (2009), S&P 500 options in Matsuda (2006); Albrecher & Schoutens(2005) or Euro Stoxx 50 (SX5E) options in Schoutens & al. (2004).22able 1: Maturities allowing that assumption 1 is satisfied, for some sets of implied NIG parameters.Other parameters: K = 4000 , r = 1% , q = 0% and S t = 3500 (OTM) or S t = 4500 (ITM).NIG parameters Accessible maturities α β δ OTM ITMSaebø (2009) 8.9932 -4.5176 1.1528 τ > . τ > . Matsuda (2006) 20.7408 -11.7308 0.2483 τ > . τ > . Schoutens & al. (2004) 16.1975 -3.1804 1.0867 τ > . τ > . Albrecher & Schoutens (2005) 18.4815 -4.8412 0.4685 τ > . τ > . In this subsection we estimate the rest of some series arising in our pricing formulas, in order to determinewhat truncation has to be applied to obtain a desired level of precision in option prices. For simplicityof notations, we perform the analysis in the symmetric model, but extension to the asymmetric case isstraightforward. Cash-or-nothing Let us observe that the general term of the cash-or-nothing series (49) is the samethan the R p +1 term introduced in (37) in the proof of formula 1: R p +1 := 1 √ π p + 1 ( − p p p ! k p +10 K p +1 ( αδτ ) (cid:18) δτ α (cid:19) − p . (101)Using the bound (38), we therefore know that, for ǫ > , there exists a rank p ǫ such that the general termof the series in the cash-or-nothing formula (49) is bounded by | R p ǫ +1 | ∼ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p πp ǫ (2 p ǫ + 2) k eαδτ (cid:18) k ( δτ ) (cid:19) p ǫ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < ǫ. (102)As a consequence of assumption 1, | k eαδτ | < and therefore, denoting by ⌈ X ⌉ the least integer greater orequal to a real number X , it suffices to choose p ǫ = & log αǫ | k δτ | ' (103)to be sure that all terms of order p ≥ p ǫ are O ( ǫ ) in the series (49). Turning back to the n -variable (i.e. n = 2 p + 1 ), it follows from (103) that, definying n ǫ := 2 p ǫ + 1 , (104)23hen all terms of order n ≥ n ǫ are O ( ǫ ) in the series of formula 3, and that the error in the option priceitself is bounded by αe ( αδ − r ) τ √ π ǫ (105)after the computation of n ǫ + 1 terms. Asset-or-nothing Recall the notations introduced in the proof of formula 1 for the general term of theseries: R n ,n := k n n !Γ(1 + − n + n ) K n − n ( αδτ ) (cid:18) δτ α (cid:19) − n n (106)and for the terms on the line n = 0 : R n := k n n !Γ(1 − n ) K n ( αδτ ) (cid:18) δτ α (cid:19) − n . (107)Let us fix n ∈ N and consider (cid:12)(cid:12)(cid:12)(cid:12) R n ,n +1 R n ,n (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Γ(1 + − n + n )Γ(1 + − n + n +12 ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K n − n ( αδτ ) K n − n ( αδτ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r δτ α . (108)From the particular values of the Gamma functions (134), the ratio of Gamma functions in (108) is smalleror equal to √ π , and the ratio of Bessel functions is smaller than 1, as a consequence of the symmetry andmonotonicity relations (140) and (141). Hence, (cid:12)(cid:12)(cid:12)(cid:12) R n ,n +1 R n ,n (cid:12)(cid:12)(cid:12)(cid:12) < r πδτ α (109)and, consequently, | R n ,n | < | R n , | = | R n | for any n in N as soon as τ < απδ . (110)Under this condition, all R n ,n terms are therefore O ( ǫ ) as soon as n , n ≥ n ǫ where n ǫ is the onedetermined in (104), and, consequently, the error in the option price given formula 1 is bounded by Kαe ( αδ − r ) τ √ π ǫ (111)24fter the computation of ( n ǫ + 1) terms. Note that if (110) is not satisfied, the series still converges butthe maximum is not attained on the line n = 0 , which complicates the estimation of the number of termsto compute. We may nevertheless observe that (110) is a very reasonable condition: for instance, usingthe implied parameters given in table 1 for SX5E options, we find τ < . , which is very close to themaximal expiry (10 years) quoted for options written on this underlying. European Exactly the same analysis can be performed on the European option, resulting in an errorfor the option price given by formula 2 bounded by Kαe ( αδ − r ) τ √ π ǫ (112)after the computation of n ǫ ( n ǫ + 1) terms (because the n summation starts at n = 1 ). To illustratethese observations, we summarize in table 2 the minimal rank, number of terms and price errors obtainedfor the digital and European options for some realistic market parameters.Table 2: Rank n ǫ beyond which the series terms in formulas 1, 2 and 3 are ǫ ) , and correspondingtruncation error on the option prices. Parameters: S t =3800, K = 4000 , r = 1% , q = 0% , τ = 1 , α = 8 . , δ = 1 . . Asset-or-nothing (Formula 1) ǫ Minimal rank n ǫ Number of terms ( n ǫ + 1) Price error − − − . × − − 10 121 . × − European (Formula 2) ǫ Minimal rank n ǫ Number of terms n ǫ ( n ǫ + 1) Price error − − − . × − − 10 110 . × − Cash-or-nothing (Formula 3) ǫ Minimal rank n ǫ Number of terms n ǫ + 1 Price error − − − . × − − 10 11 . × − .3 Comparisons with Fourier techniques Lewis formula We recall that, following Lewis (2001), digital option prices admit convenient representa-tion involving the risk-neutral characteristic function and the log-forward moneyness; the asset-or-nothingcall can be written as C a/n = S t e − qτ 12 + 1 π ∞ Z Re (cid:20) e iuk Ψ L ( u − i, τ ) iu (cid:21) d u , (113)and the cash-or-nothing call as C c/n = e − rτ 12 + 1 π ∞ Z Re (cid:20) e iuk Ψ L ( u, τ ) iu (cid:21) d u , (114)where, here, k := log S t K + ( r − q ) τ , and where the characteristic function Ψ( u, t ) = e tψ ( u ) has beennormalized by the martingale adjustment: Ψ L ( u, t ) := e iuωt Ψ( u, t ) = e iuωt + iµut − δt (cid:16) √ α − ( β + iu ) − √ α − β (cid:17) , (115)so that the martingale condition Ψ L ( − i, t ) = 1 holds true. In table 3, we compare the asset-or-nothingprices obtained by an application of formula 1 (truncated at n = n = max ) and of formula 6 (truncatedat n = n = n = max ), with a numerical evaluation of the Lewis formula (113) performed via a classicalrecursive algorithm on [0 , ] . Same comparison is made in table 4 for the cash-or-nothing prices. Weobserve the excellent agreement between our analytical result and numerical ones, as well as the fastconvergence of the series. The convergence is particularly accelerated in the ATM situation (for instancein the symmetric model, only 3 terms are needed to obtain a precision of − in the cash-or-nothingprice). It is slightly slower for deep OTM options: this is because k ∼ log S when S → , and thereforethe positive powers of k tend to slow down the overall convergence speed. Note also that the convergenceis more rapid in the symmetric than in the asymmetric model, because we choose an implied parameter | β | > complying with the calibrations in table 1; if we had chosen | β | < , then the positive powers of β would have accelerated the convergence of the asymmetric series.26able 3: Prices of asset-or-nothing call options, obtained by truncations of formulas 1 and 6, and by anumerical evaluation of (113). Parameters: K = 4000 , r = 1% , q = 0% , τ = 1 , α = 8 . , δ = 1 . . Symmetric model [ β = 0 ]Formula 1 Lewis (113) max = 3 max = 5 max = 10 max = 15 Deep OTM ( S t = 3000 ) 861.9096 796.515 804.8118 804.9099 804.9097OTM ( S t = 3500 ) 1495.76986 1493.3213 1493.5276 1493.5278 1493.5278ATM 2309.8330 2313.6169 2313.7110 2313.7110 2313.7110ITM ( S t = 4500 ) 3163.3516 3170.7414 3170.9431 3170.9431 3170.9431Deep ITM ( S t = 5000 ) 3986.4269 3999.5086 3999.8854 3999.8852 3999.8852 Asymmetric model [ β = − . ]Formula 6 Lewis (113) max = 10 max = 20 max = 30 max = 50 Deep OTM ( S t = 3000 ) 1084.9112 991.4964 990.8328 990.8302 990.8302OTM ( S t = 3500 ) 1814.0381 1705.6678 1704.8935 1704.8905 1704.8905ATM 2593.7092 2480.0154 2479.11828 2479.1149 2479.1149ITM ( S t = 4500 ) 3310.5927 3252.0495 3250.4093 3250.4089 3250.4089Deep ITM ( S t = 5000 ) 3777.9899 4003.6194 3989.4277 3989.7291 3989.7293 Carr-Madan formula Regarding European options, we recall the representation given in Carr & Madan(1999) based on the introduction of a dampling factor a to avoid the divergence in u = 0 ; namely, let Ψ CM ( u, t ) := e iu [log S t +( r − q + ω ) t ] Ψ( u, t ) , (116)then the European call price admits the representation: C eur = e − a log K − rτ π ∞ Z e − iu log K Re (cid:20) Ψ CM ( u − ( a + 1) i, τ ) a + a − u + i (2 a + 1) u (cid:21) d u, (117)where a < < a max , and a max is determined by the square integrability condition Ψ CM ( − ( a +1) i, τ ) < ∞ .In table 5 we compare the European prices obtained by formula 2 (truncated at n = n = max ) andformula 7 (truncated at n = n = n = max ), with a numerical evaluation of the Carr-Madan formula(117) on the interval [0 , ] . We also observe the excellent agreement between our analytical results andthe numerical ones, as well as the accelerated convergence for very short term options. For instance, when τ = (1 + 5) iterations are enough to obtain a precision of − in the option price in the symmetricmodel; this is because, when S t is close to K , then k ∼ ( r − q + ω ) τ and therefore when τ → thepositive powers of k arising in formulas 2 and 7 accelerate the convergence of the series. Note that, on thecontrary, the short maturity case is not a favorable situation for a numerical evaluation of the Carr-Madanformula, because of the presence of oscillations of the integrand that considerably slow down the numerical27able 4: Prices of cash-or-nothing call options, obtained by truncations of formulas 3 and 8, and by anumerical evaluation of (113). Parameters: K = 4000 , r = 1% , q = 0% , τ = 2 , α = 8 . , δ = 1 . . Symmetric model [ β = 0 ]Formula 3 Lewis (114) max = 3 max = 5 max = 10 max = 15 Deep OTM ( S t = 3000 ) 0.2127 0.2092 0.2095 0.2095 0.2095OTM ( S t = 3500 ) 0.3076 0.3073 0.3073 0.3073 0.3073ATM 0.4054 0.4054 0.4054 0.4054 0.4054ITM ( S t = 4500 ) 0.4973 0.4973 0.4973 0.4973 0.4973Deep ITM ( S t = 5000 ) 0.5793 0.5793 0.5793 0.5793 0.5793 Asymmetric model [ β = − . ]Formula 8 Lewis (114) max = 10 max = 20 max = 30 max = 50 Deep OTM ( S t = 3000 ) 0.2579 0.2360 0.2357 0.2357 0.2357OTM ( S t = 3500 ) 0.3523 0.3244 0.3240 0.3240 0.3240ATM 0.4544 0.4077 0.4074 0.4074 0.4074ITM ( S t = 4500 ) 0.5740 0.4823 0.4827 0.4827 0.4827Deep ITM ( S t = 5000 ) 0.7634 0.7277 0.5733 0.5452 0.5489Fourier inversion process.Table 5: Prices of European call options of various maturities, obtained by truncations of formulas 2 and7, and by a numerical evaluation of (117). Parameters: S t = K = 4000 , r = 1% , q = 0% , α = 8 . , δ = 1 . . Symmetric model [ β = 0 ]Formula 2 Carr-Madan (117)Maturity max = 3 max = 5 max = 10 max = 15 Asymmetric model [ β = − . ]Formula 7 Carr-Madan (117)Maturity max = 10 max = 20 max = 30 max = 50 Let n ∈ N \{ } and define the family of independent and identically distributed random variables Z ( i ) , i = 1 . . . n , all distributed according to the symmetric NIG distribution Z ( i ) ∼ NIG( α, , δ, µ ) , and define C ( i ) log := e − rτ (cid:20) log (cid:18) S t K e ( r − q + ω ) τ + Z ( i ) (cid:19)(cid:21) + = e − rτ [ k + Z ( i ) ] + (118)28s well as C ( n ) log := 1 n n X i =1 C ( i ) log . (119)We know from the strong law of large numbers that C ( n ) log converges to the price of the log call option,more precisely that C ( n ) log −→ E Q " e − rτ (cid:20) log S T K (cid:21) + | S t (120)almost surely when n → ∞ . Similarly, regarding power options, we define (in the European case): C ( i ) pow := e − rτ h S a e a (( r − q + ω ) τ + Z ( i ) ) − K i + , C ( n ) pow := 1 n n X i =1 C ( i ) pow (121)and, for the capped digital option, C ( i ) capped c/n := e − rτ {− k , − Monte Carlo (119) Formula 5 n = 100 n = 500 n = 1000 n max = 1 n max = 3 n max = 5 OTM ( S t = 3500) S t = 4000) S t = 4500) Power option (European) Monte Carlo (121) Formula 4 n = 100 n = 1000 n = 5000 max = 20 max = 40 max = 60 OTM ( S t = 3500) S t = 4000) S t = 4500) Capped option (digital) Monte Carlo (122) Series (77) n = 100 n = 1000 n = 5000 n max = 1 n max = 5 n max = 10 OTM ( S t = 3500) S t = 4000) S t = 4500) In this paper, we have proved two general formulas for pricing arbitrary path independent instrumentsin the exponential NIG model, in the symmetric and asymmetric cases. These formulas allow to expressthe Mellin transform of the instrument’s price as the product of the Mellin transform of the instrument’spayoff and of the NIG probability density. Inverting the formulas by means of residue theory in C and C n has allowed us to derive practical closed-form pricing formulas for various path independent optionsand contracts, under the form of quickly convergent series. The convergence of the series is guaranteed assoon as a simple condition of the log forward moneyness and on the option’s maturity is fulfilled. We havetested our results by comparing them with classical numerical methods, and provided precise estimate forthe convergence speed; notable feature is that a very reasonable number of terms is required to obtain anexcellent level of precision, and that the convergence is particularly fast for short term and at-the-moneyoptions.Future work should include, among others, an extension of the Mellin residue summation method topath independent instruments on several assets, and to path dependent instruments. Asian options withcontinuous geometric payoffs, in particular, should be investigated, because the characteristic function forthe geometric average is known exactly in the exponential NIG model (see Fusai & Meucci (2008)), for30oth fixed and floating strikes.Extension of the technique to Generalized Hyperbolic (GH) Lévy motions (see Prause (1999); Eberlein(2001)) should also be considered; indeed, the probability density of the GH distribution has a very similarform to the NIG density (1), which, at first sight, allows for the same convenient representation in termsof Mellin-Barnes integrals for the Bessel kernel. However, GH distributions are not convolution-closed,that is, the Lévy processes they generate are not necessarily distributed according to a GH distribution forincrements of length t = 1 (exceptions being the NIG process, which, as we know, is distributed accordingto a NIG distribution NIG( α, β, δt, µt ) for all t , as well as the generalized Laplace or Variance Gammadistribution). As a consequence, a Mellin-Barnes integral representation for the density of the GH processis not straightforward to derive, but could nevertheless be obtained from the moment generating function,after suitable transformations from the Laplace space to the Mellin space. Acknowledgments The author thanks Ryan McCrickerd for insightful comments and discussions. References Abramowitz, M. and Stegun, I., Handbook of Mathematical Functions, Dover Publications, Mineola, NY(1972)Aguilar, J. Ph., On expansions for the Black-Scholes prices and hedge parameters, Journal of MathematicalAnalysis and Applications , 973-989 (2019)Aguilar, J. Ph. and Korbel, J., Simple Formulas for Pricing and Hedging European Options in the FiniteMoment Log-Stable Model, Risks , 36 (2019)Aguilar, J. Ph., Some pricing tools for the Variance Gamma model, International Journal of Theoreticaland Applied Finance (2020)Albrecher, H. and Predota, M., On Asian option pricing for NIG Lévy processes, Journal of Computationaland Applied Mathematics , 153-168 (2004)Albrecher, H. and Schoutens, W., Static hedging of Asian options under stochastic volatility models usingFast Fourier Transform. In: A. Kyprianou et al. (Eds), Exotic Options and Advanced Lévy models pp.129-148, John Wiley & Sons, Hoboken, NJ (2005)Andrews, L.C., Special Functions of Mathematics for Engineers, McGraw-Hill Book Company, New York(1992)Barndorff-Nielsen, O., Exponentially decreasing distributions for the logarithm of particle size, Proceedingsof the Royal Society of London , 401-419 (1977)Barndorff-Nielsen, O., Normal inverse Gaussian distributions and the modeling of stock returns, Researchreport no 300, Department of Theoretical Statistics, Aarhus University (1995)31arndorff-Nielsen, O., Normal inverse Gaussian distributions and stochastic volatility models, Scandina-vian Journal of Statistics , 1-133 (1997)Bateman, H., Tables of Integral Transforms (vol. I and II), McGraw-Hill Book Company, New York (1954)Bertoin, J., Lévy Processes, Cambridge University Press, Cambridge, New York, Melbourne (1996)Black, F. and Scholes, M., The Pricing of Options and Corporate Liabilities, Journal of Political Economy , 637-654 (1973)Carr, P. and Madan, D., Option valuation using the Fast Fourier Transform, Journal of ComputationalFinance , 61-73 (1999)Carr, P., Geman, H., Madan, D., Yor, M., The Fine Structure of Asset Returns: An Empirical Investiga-tion, Journal of Business , 305-332 (2002)Carr, P. and Wu, L., The Finite Moment Log Stable Process and Option Pricing, The Journal of Finance , 753-777 (2003)Carr, P. and Wu, L., Time-changed Lévy processes and option pricing, Journal of Financial Economics , 113-141 (2004)Carr, P., Lee R. and Wu, L., Variance swaps on time-changed Lévy processes, Finance and Stochastics , 335-355 (2012)Cont, R. and Tankov, P., Financial Modelling with Jump Processes, Chapman & Hall, New York (2004)Eberlein, E. and Keller,U., Hyperbolic distributions in finance, Bernoulli , 281-299 (1995)Eberlein, E., Application of Generalized Hyperbolic Lévy Motions to Finance. In: Lévy Processes,Barndorff-Nielsen O.E., Resnick S.I., Mikosch T. (eds), Birkhauser, Boston, MA (2001)Fang, F. and Oosterlee, C.W., A novel pricing method for European options based on Fourier cosine seriesexpansions, SIAM Journal on Scientific Computing , 826-848 (2008)Figueroa-López, J.E., Lancette, S.R, , Lee, K. and Mi, Y., Estimation of NIG and VG models for highfrequency financial data. In: Handbook of Modeling High-Frequency Data in Finance, F. Viens, M.C.Mariani, I. Florescu (eds.), John Wiley & Sons, Hoboken, NJ (2012)Flajolet, P., Gourdon, X. and Dumas, P., Mellin transforms and asymptotics: Harmonic sums, TheoreticalComputer Science , 3-58 (1995)Fusai, G. and Meucci, A., Pricing discretely monitored Asian options under Lévy processes, Journal ofBanking & Finance , 2076–2088 (2008)Giannone, D., Reichlin, L. and Small, D., Nowcasting: The real-time informational content of macroeco-nomic data, Journal of Monetary Economics , 665-676 (2008)Glasserman, P., Monte Carlo methods in financial engineering, Springer Science & Business Media Vol.53,New York (2004)Hanssen, A. and Øigård, T.A., The Normal inverse Gaussian distribution: a versatile model for heavy-tailed stochastic processes, Proceedings - ICASSP, IEEE International Conference on Acoustics, Speechand Signal Processing , 3986-3988 (2001)Heston, S., A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond andCurrency Options, The Review of Financial Studies , 327-343 (1993)32vanov, R.V., Closed Form Pricing of European Options for a Family of Normal Inverse Gaussian Processes,Journal of Stochastic Models , 435-450 (2013)Kirkby, J. L., Efficient Option Pricing by Frame Duality with the Fast Fourier Transform, SIAM Journalon Financial Mathematics , 713-747 (2015)Lewis, A.L., A simple option formula for general jump-diffusion and other exponential Lévy processes,Available at SSRN: https://ssrn.com/abstract=282110 (2001)Luciano, E. and Semeraro, P., Multivariate time changes for Lévy asset models: Characterization andcalibration, Journal of Computational and Applied Mathematics , 1937-1953 (2010)Madan, D., Carr, P. and Chang, E., The Variance Gamma Process and Option Pricing, European FinanceReview , 79-105 (1998)Mandelbrot, B., The Variation of Certain Speculative Prices, The Journal of Business , 384-419(1963)Matsuda, K., Calibration of Lévy Option Pricing Models: Applications to S& P 500 Futures option, PhDThesis City University of New York (2006)Mechkov, S., Fast-Reversion Limit of the Heston Model, Available at SSRN:https://ssrn.com/abstract=2418631 (2015)Mittnik, S. and Rachev, S., Stable Paretian models in finance, John Wiley & Sons, Hoboken, NJ (2000)Neuberger, A., The log contract, Journal of Portfolio Management , 74-80 (1994)Papantoleon, A., An introduction to Lévy Processes with applications in finance, arXiv:0804.0482 (2008)Prause, K, The generalized hyperbolic model: estimation, financial derivatives and risk measures, PhDthesis, Institut für Mathematische Statistik, Albert-Ludwigs-Universität Freiburg (1999)Rachev, S., Kim, Y., Bianchi, M., Fabozzi, F., Financial models with Lévy processes and volatility clus-tering, John Wiley & Sons, Hoboken, NJ (2011)Ribeiro, C. and Webber, N., A Monte Carlo Method for the Normal Inverse Gaussian Option ValuationModel using an Inverse Gaussian Bridge, City University preprint (2003)Rydberg, T., The Normal inverse Gaussian Lévy process: simulation and approximation, Communicationsin Statistics. Stochastic Models , 887-910 (1997)Saebø, K., Pricing Exotic Options with the Normal Inverse Gaussian Market Model using Numerical PathIntegration, Master’s Thesis Norwegian University of Science and Technology (2009)Schoutens, W., Lévy processes in finance: pricing financial derivatives, John Wiley & Sons, Hoboken, NJ(2003)Schoutens, W., Simons, E., and Tistaert, J., A perfect calibration! Now what?, Wilmott magazine (2004)Su, Y., and Fu., M.C., Importance sampling in derivative securities pricing, 2000 Winter SimulationConference Proceedings Vol. 1. (2000)Taleb, N.N., The Black Swan: The Impact of the Highly Improbable, Random House Publishing Group,New York (2010) 33ankov, P., Pricing and Hedging in Exponential Lévy Models: Review of Recent Results. In: Paris-Princeton Lectures on Mathematical Finance. Lecture Notes in Mathematics, vol 2003, Springer, Berlin,Heidelberg (2010)Venter, J. and de Jongh, P., Risk estimation using the Normal inverse Gaussian distribution, The Journalof Risks , 1-25 (2002)Wilmott, P., Paul Wilmott on Quantitative Finance, Wiley & Sons, Hoboken, NJ, 2006Zeng, P. and Kwok, Y.K., Pricing barrier and Bermudan style options under time-changed Lévy processes:fast Hilbert transform approach, SIAM Journal on Scientific Computing , B450-B485 (2014) A Brief review of the Mellin transform We present an overview of the one-dimensional Mellin transform; this theory is explained in full detail inFlajolet et al. (1995), and table of Mellin transforms can be found in any monograph on integral transforms(see e.g. Bateman (1954)).1. The Mellin transform of a locally continuous function f defined on R + is the function f ∗ defined by f ∗ ( s ) := ∞ Z f ( x ) x s − d x. (124)The region of convergence { α < Re ( s ) < β } into which the integral (124) converges is often called thefundamental strip of the transform, and sometimes denoted < α, β > .2. The Mellin transform of the exponential function is, by definition, the Euler Gamma function: Γ( s ) = ∞ Z e − x x s − d x (125)with strip of convergence { Re ( s ) > } . Outside of this strip, it can be analytically continued, except atevery negative s = − n integer where it admits the singular behavior Γ( s ) ∼ s →− n ( − n n ! 1 s + n , n ∈ N . (126)In table 7 we summarize the main Mellin transforms used in this paper, as well as their convergence strips.3. The inversion of the Mellin transform is performed via an integral along any vertical line in the strip34able 7: Mellin pairs used throughout the paper. f ( x ) f ∗ ( s ) Convergence strip e − ax a − s Γ( s ) < , ∞ >e − ax − a − s Γ( s ) < − , >K ν ( ax ) a − s s − Γ (cid:0) s − ν (cid:1) Γ (cid:0) s + ν (cid:1) < | Re ( ν ) | , ∞ > K ν ( a √ x + b )( x + b ) ν a s s − b s − ν Γ( s ) K ν − s ( ab ) < , ∞ > of convergence: f ( x ) = c + i ∞ Z c − i ∞ f ∗ ( s ) x − s d s iπ c ∈ ( α, β ) (127)and notably for the exponential function one gets the so-called Cahen-Mellin integral : e − x = c + i ∞ Z c − i ∞ Γ( s ) x − s d s iπ , c > . (128)4. When f ∗ ( s ) is a ratio of products of Gamma functions of linear arguments: f ∗ ( s ) = Γ( a s + b ) . . . Γ( a m s + b m )Γ( c s + d ) . . . Γ( c l s + d l ) (129)then one speaks of a Mellin-Barnes integral , whose characteristic quantity is defined to be ∆ = m X k =1 a k − l X j =1 c j . (130) ∆ governs the behavior of f ∗ ( s ) when | s | → ∞ and thus the possibility of computing (127) by summingthe residues of the analytic continuation of f ∗ ( s ) right or left of the convergence strip: ∆ < f ( x ) = − X Re ( s ) >β Res (cid:2) f ∗ ( s ) x − s (cid:3) , ∆ > f ( x ) = X Re ( s ) <α Res (cid:2) f ∗ ( s ) x − s (cid:3) . (131)For instance, in the case of the Cahen-Mellin integral one has ∆ = 1 and therefore: e − x = X Re ( s ) < Res (cid:2) Γ( s ) x − s (cid:3) = ∞ X n =0 ( − n n ! x n (132)as expected from the usual Taylor series of the exponential function.35 Some useful special functions identities We list some properties of special functions that are used throughout the paper; more details can be founde.g. in Abramowitz & Stegun (1972); Andrews (1992). B.1 Gamma function Particular values The Gamma function Γ( s ) has been defined in (125) for Re ( s ) > ; integrating byparts shows that it satisfies the functional relation Γ( s + 1) = s Γ( s ) ; as Γ(1) = 1 , it follows that Γ( n + 1) = n ! , n ∈ N (133)and that the analytic continuation of Γ( s ) to the negative half-plane is singular at every negative integer − n with residue ( − n n ! . Other useful identities include Γ( ) = √ π and, more generally, Γ (cid:18) − n (cid:19) = ( − n n n !(2 n )! √ π Γ (cid:18) 12 + n (cid:19) = (2 n )!4 n n ! √ π. (134)for n ∈ N . Stirling approximation We recall the well-known Stirling approximation for the factorial: n ! ∼ n →∞ √ πn n n e − n . (135) Legendre duplication formula For any s ∈ C , we have: Γ (cid:0) s (cid:1) Γ( s ) = √ π s − (cid:0) s +12 (cid:1) . (136) Pochhammer symbol The Pochhamer symbol ( a ) n , sometimes denoted by the Appel symbol ( a, n ) ,and also called rising factorial, is defined by ( a ) n := Γ( a + n )Γ( a ) , a / ∈ Z − . (137)36he definition (137) extends continuously to negative integers thans to the functional relation Γ( s + 1) = s Γ( s ) , thanks to the relation: ( − k ) n = ( − n k !( k − n )! 0 ≤ n < k n > k. (138)where k ∈ N . B.2 Bessel functions The modified Bessel function of the second kind, also called MacDonald function, can be defined by theMellin integral K ν ( z ) := 12 (cid:16) z (cid:17) ν ∞ Z e − t − z t t − ν − d t (139)for | arg z | < π . It follows that K ν ( z ) has the symmetry property: K ν ( z ) = K − ν ( z ) (140)and has monotonous absolute values: ≤ ν < ν = ⇒ | K ν ( z ) | < | K ν ( z ) | . (141) Large index When ν → ∞ , one has the following behavior: K ν ( z ) ∼ ν →∞ r π ν (cid:16) ez ν (cid:17) − ν . (142) Large argument (Hankel’s expansion) Define the following sequence: a ( ν ) = 1 a k ( ν ) = (4 ν − )(4 ν − ) . . . (4 ν − (2 k − ) k !8 k , k ≥ . (143)Then, for large z and fixed ν , we have: K ν ( z ) = z →∞ r π z e − z ∞ X k =0 a k ( ν ) z k . (144)37n particular, when ν − , i.e. when ν = , all the a k ( ν ) are null in definition (143) when k ≥ ,and we are left with: K ( z ) = r π z e − zzK + } | S t (cid:3) = 1 − E Q (cid:2) { K −