1 Introduction↩︎

The absence of normality in stock returns has been observed since at least the pioneering essay by Bachelier [1] in the early twentieth century and has motivated the introduction of stochastic processes that go beyond the Brownian motion (notably in terms of higher-order moments) to model the logarithm of stock prices. Among these processes, Lévy processes (see the classical monographs [2], [3] among many other references) play a leading role as they capture many stylized facts such asymmetry in the distributions and fat tails, while preserving the independence and stationarity of financial returns. Among the most popular Lévy processes in quantitative finance, one finds the Variance Gamma (VG) process of [4], the bilateral Gamma (BG) process of [5], the CGMY process of [6] or the KoBoL process of [7]; more recent extensions also include additive processes (see [8], [9]) that relax the stationarity hypothesis, or convolutions of generalized inverse Gaussian subordinators in [10] for instance.

VG, BG, CGMY and KoBoL processes actually belong to a broader class of processes, known as tempered stable (TS) processes (see [11]), that are also sometimes called generalized tempered stable (GTS) processes and generalize the so-called truncated Lévy flights introduced in Physics in [12] (see also [13]) with the purpose to capture two different behaviors: a central part that is approximately Lévy stable, and tails that are approximately exponential in order to ensure finiteness of the second moment. Since then, TS processes have also proved popular in finance notably in the framework of exponential Lévy pricing models as described in [14] (see also more information on the change of measure and Esscher transform for TS processes in [15]), thanks to their tractability and the good empirical results they provide. Among the many financial applications of TS processes, one can mention energy market modeling and pricing of associated derivatives in [16], [17], incorporation of stochastic volatility in [18], [19] or multivariate extensions and basket option pricing in [20], [21]. Equity indices calibration has also been conducted in [22], and simulation algorithms for TS processes have been provided in [23]–[26].

One of the reasons for the popularity of TS processes (and of Lévy processes in general) in option pricing applications is the availability of their characteristic function in closed form, allowing for the usage of Fourier pricing methods, that are easy to implement and benefit from many refinements that have been achieved to increase their performance. In particular, the introduction of the Fast Fourier Transform (FFT) in Finance in [27] paved the way to the current state-of-the-art Fourier methods such as the COS and PROJ methods (see [28] and [29] respectively), and extensive research has been conducted to extend the methods to exotic options (see [30]–[32]). However, Fourier related methods are known to be subject to numerical accuracies stemming from multiple valued function such as complex square roots or logarithms (see [33], [34]); moreover, at least one prefixed parameter has to be selected to perform the computations. For the FFT, the option price has to be modified (dampened) to ensure square integrability, and the convergence of the method is highly sensitive to the choice of the damping parameter, in particular in the multivariate case (see [35], where optimization of damping factors is discussed to accelerate the convergence of numerical quadrature methods in the context of Fourier pricing). For the COS method, no damping parameters are required, however a truncation range has to be selected, into which the log returns density is approximated by a (finite) cosine expansion; this truncation is generally chosen through a “rule of thumb" based on the first four cumulants of the distribution, that can actually lead to serious mispricings, see [36]; in some cases, the numerical inaccuracies of the FFT and COS methods (and similar Fourier techniques like the Lewis-Lipton method [37]) have also been observed to compensate the model error itself, giving birth to what has been called in [38], [39] a”ghost calibration" and to spurious volatility smiles and skews. Similarly, the PROJ method also does not require the selection of a control parameter and displays great efficiency and robustness (see for instance [30], [32]), however, like in the COS case, it necessitates the introduction of a truncation for the density support, that is also determined through a cumulant-based rule and needs to be adapted when it comes to very precise option pricing; additionally, and like in the COS case, PROJ pricing also includes a series truncation, corresponding to the discretization of the integration interval and controlling the precision of the characteristic function’s projection over the Riesz basis.

The purpose of this paper is therefore to develop an option pricing method under TS processes that minimizes the numerical disadvantages described above, and allows for an arbitrary degree of precision (up to machine capabilities) within a competitive computational time. The solution we propose is based on series expansions generated by multidimensional residues associated to Mellin-Barnes integrals. We note that the Mellin transform has already been used for option pricing purposes (notably via transformations of partial differential equations) for instance in [40]–[43], but the combination with residue calculus is more recent and dates back to [44] in the context of path-independent contracts under exponential Lévy models. Refining this technique, we are able to derive new and fast converging series representations for TS distributions and also for digital and European option prices, and we recover and improve the formulas previously obtained in [44] and [45] in the VG, BG and one-sided cases. Moreover, our method has the significant advantage of being hyperparameter-free, the precision of the computed prices depending only on the chosen summation order.

The paper is organized as follows: we first recall in section 2 the definition of TS distributions, their associated Lévy processes as well as their particular cases. Then, in section 3, we provide series expansions of the TS distributions by means of multidimensional Mellin transform inversion; these formulas are new to the literature (they only existed in the single-sided case so far). In section 4, we derive series expansions for digital and European option under TS processes and, in section, 5, we focus on the particular cases covered by the TS model, and derive simpler notable for the BG and VG models. In section 6, we adapt our methodology to the one-sided case and provide particularly simple pricing formulas, taking the form of series expansions; finally, we conduct in section 7 detailed numerical analysis and provide comparison with state-of-the-art pricing methods. To further aid in research and reproducibility, the implementation code is provided in the publicly available Github repository TS-Pricing¹ (written in Python). Section 8 is dedicated to conclusive remarks.

2 TS distributions↩︎

2.1 Notations↩︎

In all of the following, we will consider a probability space \((\Omega,\mathcal{F},\mathbb{P})\) and, unless otherwise mentioned, expectations will be taken under \(\mathbb{P}\), that is \(\mathbb{E} [.] = \mathbb{E} ^{\mathbb{P}} [.]\). Given a random variable \(X:\Omega \rightarrow \mathbb{R}\), we define its characteristic function (Fourier transform) as \[\forall u\in{\mathbb{R}},\quad \varphi_X (u) = \mathbb{E} \left[ e^{\mathrm{i}u X} \right] = \int_{\mathbb{R}} e^{\mathrm{i}u x} f_X(x) \mathrm{d}x ,\] where the second equality holds if we assume that \(X\) possesses a \(\mathbb{P}\)-density \(f_X\). In the following, we will also denote by \(\varphi_X\) the analytic continuation of the characteristic function whenever it exists.

2.2 One-sided and double-sided tempered stable distributions↩︎

We start by recalling some essential facts about \(\mathrm{TS}\) distributions. The reader can find further details in [11] (as well as multiple other references therein).

2.2.1 One-sided TS distributions↩︎

One-sided \(\mathrm{TS}\) distributions are probability distributions supported by \(\mathbb{R}_+\) and whose characteristic function \(\varphi_{\mathrm{TS}_+}\) is of the form: \[\forall u\in{\mathbb{R}},\quad\varphi_{\mathrm{TS}_+}(u) =\exp\left( \int_{\mathbb{R}}\left(e^{\mathrm{i}u x}-1\right)\nu_{\mathrm{TS}_+}(x)\mathrm{d}x \right)\] where the Lévy measure \(\nu_{\mathrm{TS}_+}\) is given by: \[\label{Levy95measure95TS43} \forall x\in{\mathbb{R}}, \quad \nu_{\mathrm{TS}_+}(x) = \alpha\frac{e^{-\lambda x}}{x^{1+\beta}}\mathbb{1}_{(0,\infty)}(x) ,\tag{1}\] for \(\alpha,\lambda\in(0,\infty)^2\) and \(\beta\in [0,1)\). We will denote such distributions by \(\mathrm{TS}_+(\alpha,\lambda,\beta)\), and their density by \(f_{\mathrm{TS}_+}(x;\alpha,\beta,\lambda)\); when there is no ambiguity on the parameters and to simplify the notations, we will simply write \(f_{\mathrm{TS}_+}(x)\). If furthermore \(\beta\neq 0\), then the characteristic function of a random variable \(X\sim \mathrm{TS}_+(\alpha,\beta,\lambda)\) is known under closed form, and can be written explicitly as: \[\label{char95TS43} \forall u\in{\mathbb{R}}, \quad\varphi_{\mathrm{TS}_+}(u) = \exp\left({\alpha\Gamma(-\beta)\left((\lambda-\mathrm{i}u)^{\beta}-\lambda^{\beta}\right)}\right).\tag{2}\] As \(\mathrm{TS}_+(\alpha,\beta,\lambda)\) distributions are infinitely divisible, they generate a Lévy process \(X^+\) whose characteristic function is: \[\forall (u,t)\in{\mathbb{R}}\times(0,\infty),\quad\varphi_{+}(u,t) = \exp\left(t\alpha\Gamma(-\beta)\left((\lambda-\mathrm{i}u)^{\beta}-\lambda^{\beta}\right)\right).\] The process \(X^+\) is an almost surely non decreasing process (i.e., it is a subordinator), whose increments are independent and distributed according to: \[\forall t\geq s\geq 0, \quad X_t^+-X_s^+\sim \mathrm{TS}_+(\alpha (t-s),\lambda,\beta).\]

2.2.2 Double-sided TS distributions↩︎

We now define the double-sided TS distribution or, simply, the TS distribution, by the following convolution: \[\label{TS95def} \mathrm{TS}(\alpha_+,\beta_+,\lambda_+,\alpha_-,\beta_-,\lambda_-) = \mathrm{TS}_+(\alpha_+,\beta_+,\lambda_+) \star \mathrm{TS}_-(\alpha_-,\beta_-,\lambda_-)\tag{3}\] where \(\alpha_\pm, \lambda_\pm \in {\mathbb{R}}\) and \(\beta_\pm \in [0,1)\). We note that some authors also speak of generalized tempered stable distributions (GTS), see for instance [3]. As for one sided TS distributions, TS distributions are infinitely divisible and, as such, they generate a Lévy process; it follows from definition 3 and from 2 and ?? that its characteristic function is, for \(\beta_\pm \in (0,1)\), \[\begin{gather} \label{eq:char95TS} \forall (u,t)\in{\mathbb{R}}\times(0,\infty),\quad\varphi_{\mathrm{TS}}(u,t) = \exp\left(t\alpha_+\Gamma(-\beta_+)\left((\lambda_+-\mathrm{i}u)^{\beta_+}-\lambda_+^{\beta_+}\right)\right.\\ +\left.t\alpha_-\Gamma(-\beta_-)\left((\lambda_-+\mathrm{i}u)^{\beta_-}-\lambda_-^{\beta_-}\right) \right). \end{gather}\tag{4}\]

Extension to the case \(1< \beta_\pm < 2\) has been considered (see discussion in [11]), however we will stick to the initial hypothesis \(\beta_\pm\in[0,1)\) because it covers most cases of financial significance (see section 2.4 thereafter). Last, it follows from definition 3 that the density of a \(\mathrm{TS}\) distribution admits the representation \[\label{TSdensity95conv} f_\mathrm{TS}(x;\alpha_+,\beta_+,\lambda_+,\alpha_-,\beta_-,\lambda_-) = \int_0^\infty f_{\mathrm{TS}_+}(x+y;\alpha_+,\beta_+,\lambda_+)f_{\mathrm{TS}_+}(y;\alpha_-,\beta_-,\lambda_-)\mathrm{d}y.\tag{5}\] for all \(x\in{\mathbb{R}}\).

2.3 A continuity property↩︎

Let us derive an interesting continuity property that will be useful in the option pricing context

2.4 Particular cases↩︎

TS distributions recover many well-known distributions that are popular within the option pricing and financial engineering space; below are some examples.

• \(\beta_+=\beta_-\) in 4 yields the characteristic function of the KoBoL distribution (see [7], [47]).

• If \(\alpha_+=\alpha_-=\alpha\) and \(\beta_+=\beta_-=\beta\) then the TS distribution becomes a CGMY (Carr-Geman-Madan-Yor) distribution, where \(C:=\alpha\), \(G:=\lambda_-\), \(M=\lambda_+\), \(Y=\beta\). We note, moreover, that the case \(Y\in(0,1)\) is the case where the associated Lévy process is of infinite activity and finite variation, a situation which is typical for risk-neutral processes of most equity indices for instance (see [6]).

• If \(\beta_\pm=0\), one speaks of a bilateral Gamma (BG) distribution (see [5]) and in that case the characteristic function becomes: \[\label{eq:char95BG} \forall u\in{\mathbb{R}},\quad\varphi_{\mathrm{BG}}(u) = \left( \frac{\lambda_+}{\lambda_+ - \mathrm{i}u} \right)^{\alpha_+} \left( \frac{\lambda_-}{\lambda_+ + \mathrm{i}u} \right)^{\alpha_-} .\tag{7}\]

• When \(\alpha_+=\alpha_-=\alpha\) and \(\beta_+=\beta_-=0\), one recovers the Variance Gamma (VG) distribution of [4]; when furthermore \(\lambda_+=\lambda_-\), the distribution is symmetric, and was introduced in [48].

3 Series expansion of TS distributions↩︎

The purpose of this section is to derive series representations for single and double-sided TS distributions. Our approach is the following: first, we derive a Mellin-Barnes (MB) integral representation for the corresponding density function, and, in a second time, we apply residue calculus to express these integrals under the form of residues series expansions. Of course, series (and asymptotic) representations for Lévy densities have already been discussed in the literature but they often focus on Gram-Charlier expansions (see for instance [49], [50] and references therein) or on expansions for the cumulative distribution function (see [51]). Mellin residue calculus for TS densities has been utilized in [52] to get a density expansion in the one-side case, that we will recover in proposition 4; the general double-sided series expansion we will obtain in proposition 5 is the first of its kind (to the best of our knowledge).

3.1 Notations↩︎

As a preamble, we first introduce notations that will be used extensively in the article, especially for writing multiple MB integrals (that is, contour integrals involving ratios of products of gamma functions of linear arguments, see [46] or any monograph on special functions) in a lighter way.

To this aim, let us first introduce some quantities. Let \(n\in\mathbb{N}\) be the dimension of the MB integral (\(n=2\) for a double integral, \(n=3\) for a triple integral etc.), and let \((m_u,m_d)\in\mathbb{N}^2\) denote the number of gamma functions in the numerator and in the denominator respectively. We also define two other integers \((\ell_u,\ell_d)\in\left\llbracket 1, n \right\rrbracket^2\) denoting the dimension of the linear span generated by the complex variables involved in the numerator and denominator respectively. Let us finally consider the two matrices \((U,D)\in\mathcal{M}_{m_u,\ell_u}({\mathbb{R}})\times\mathcal{M}_{m_d,\ell_d}({\mathbb{R}})\) and the two vectors \((\boldsymbol{u},\boldsymbol{d})\in{\mathbb{R}}^{m_u}\times{\mathbb{R}}^{m_d}\) containing the coefficients of the arguments of the gamma functions, and let us define the integration subset \[\boldsymbol{c}+\mathrm{i}{\mathbb{R}}^n := \{ \boldsymbol{z}\in{\mathbb{C}}^n, z_1 = c_1 + \mathrm{i}y_1, \dots, z_n = c_n + \mathrm{i}y_n \} , \quad (y_1,\dots,y_n) \in {\mathbb{R}}^n.\] With these notations, an \(n\)-dimensional MB integral can be written as: \[\int\limits_{\boldsymbol{c}+\mathrm{i}{\mathbb{R}}^n} \boldsymbol{\Gamma}^{\boldsymbol{\pi}}\left[U\boldsymbol{s}_{\leqslant\ell_u}+\boldsymbol{u}|~D\boldsymbol{s}_{\leqslant\ell_d}+\boldsymbol{d}\right] \prod_{i=1}^n q_i^{s_i}\frac{\mathrm{d}\boldsymbol{s}}{(\mathrm{i}2\pi)^n}\] where \(\boldsymbol{q}\in{\mathbb{R}}^n\), \(\boldsymbol{s}_{\leqslant\ell_u}:=(s_1,...,s_{\ell_u})\), \(\boldsymbol{s}_{\leqslant\ell_d}:=(s_1,...,s_{\ell_d})\), \(\mathrm{d}\boldsymbol{s} := \mathrm{d}s_1 \wedge ... \wedge \mathrm{d}s_n\), and where \(\boldsymbol{\Gamma}^{\boldsymbol{\pi}}\left[\cdot|~\cdot\right]\) is defined as: \[\forall (\boldsymbol{v},\boldsymbol{w})\in {\mathbb{R}}^{m_u}\times{\mathbb{R}}^{m_d}, \quad \boldsymbol{\Gamma}^{\boldsymbol{\pi}}\left[\boldsymbol{v}|~\boldsymbol{w}\right] := \frac{\prod_{i=1}^{m_u} \Gamma\left(v_i\right)}{\prod_{i=1}^{m_d} \Gamma\left(w_i\right)}.\] Finally, for a matrix \(M\) with \(n_M\) lines, and for a set of \(k\) indices \(\{i_m\}_{m\in\left\llbracket 1, k \right\rrbracket}\in\left\llbracket 1, n_M \right\rrbracket^{k}\), we will denote \(M_{\{i_1,...,i_k\}}\) the sub-matrix composed of the lines \(\{i_1,...,i_k\}\). In order to clarify this formalism, we provide an explanatory illustration in example 1.

3.2 Mellin-Barnes representation of TS densities↩︎

We give here the MB representations for the densities of the single and double-sided TS distributions. We start with the single-sided case, for which as MB representation was already obtained in [52] with a slightly different approach.

From now on, we associate to a \(\mathrm{TS}(\alpha_+,\beta_+,\lambda_+,\alpha_-,\beta_-,\lambda_-)\) distribution the quantities \(a_\pm\), \(\underaccent{\bar}{\lambda}\) and \(\gamma\) defined as: \[\label{eq:constants-ts} \begin{cases} a_+ := -\alpha_+ \Gamma(-\beta_+),\\ a_- := -\alpha_- \Gamma(-\beta_-),\\ \underaccent{\bar}{\lambda} := \lambda_+ + \lambda_-,\\ \gamma := a_+\lambda_+^{\beta_+}+a_-\lambda_-^{\beta_-}. \end{cases}\tag{8}\]

3.3 Series expansion of TS densities↩︎

Thanks to the MB representations that we have established in propositions 2 and 3, we will now be able to derive series representations for the single and double-sided TS distributions using residue calculus. We start with the one-sided case, for which we are able to recover a more compact version of the expression already obtained in [52].

Let us now extend proposition 4 to the more general double-sided TS case. To that extent, and whenever double-sided TS distributions are considered throughout the rest of the paper, we will impose the additional condition \((\beta_+,\beta_-)\in ((0,1)\cap({\mathbb{R}}\backslash{\mathbb{Q}}))^2\). The irrationality of the \(\beta_\pm\) parameters is not mandatory in theory, but it allows to avoid multiple poles in the arising MB integrals, and thus allows for simpler expressions for the residue series. Furthermore, we note that, with the pointwise convergence established in proposition 1, this condition appears to be absolutely not restrictive thanks to the density of \({\mathbb{Q}}\) in \({\mathbb{R}}\).

In figure 1, we compare the series expansion from 5 (truncated at \(\boldsymbol{n}=(60,60,60)\) with the numerical Fourier inversion of the TS characteristic function 4 with \(t=1\); we can observe that both approaches display excellent agreement.

4 Option pricing under TS processes↩︎

TS distributions have already been widely used in financial modeling context (see for example [11], [21], [58], [59]). The availability of the characteristic function in closed form allows to use Fourier pricing techniques such as Lewis [37] or Carr-Madan [27], as well as their many refinements such as Hilbert [60], COS [28] and PROJ [29], the latter being known to be particularly fast and efficient. All of these methods feature one (or several) hyperparameters to be optimally selected before any pricing procedure, a drawback that is of course undesirable in practice. In this section, we will therefore focus on providing hyperparameter free TS pricing for digital and European options based on residue series expansion; to that end, we will first derive MB representations for the option prices and then compute them use residue calculus.

4.1 Option pricing framework↩︎

Let us consider a filtered probability space \((\Omega, \mathcal{F}, \left(\mathcal{F}_t\right)_{t\geqslant 0}, \mathbb{P})\) where \(\left(\mathcal{F}_t\right)_{t\geqslant 0}\) denotes the natural filtration generated by the TS process \(\left(X_t\right)_{t\geqslant 0}\). We assume that the price at time \(t\geqslant 0\) of some financial asset can be written as: \[\label{St95dynamics} S_t := S_0e^{(r-q+\zeta)t + X_t} ,\tag{12}\] where \(S_0>0\), \(r\) and \(q\) are respectively the risk-free rate and the dividend yield (both assumed to be continuously compounded) and \(\zeta\) is the martingale adjustment (or compensator) chosen in such a way that the discounted process \((e^{-(r-q)t}S_t)_{t\geqslant 0}\) is a martingale. As we know that, for any real valued Lévy process \(Y\), \(e^Y\) is a martingale if and only if \(\mathbb{E} \left[e^{Y_1}\right]=1\) (see [3]), it is easy to see from 4 that the discounted price process is a martingale if and only if \(\lambda_+ > 1\) (a condition that we will require throughout the rest of the paper) and \[\label{eq:def-zeta} \zeta = -\ln \varphi_\mathrm{TS}(-\mathrm{i}) = a_+ ( (\lambda_+-1)^{\beta_+} - \lambda_+^{\beta_+} ) + a_- ( (\lambda_-+1)^{\beta_-} - \lambda_-^{\beta_-} ) .\tag{13}\] Last, we recall that at time \(t\geqslant 0\), the price \(C_t\) of a contingent claim delivering a (integrable) payoff \(\mathcal{P}(S_T)\) at maturity \(T\geqslant t\) is given by the conditional expectation \[\label{pricing} C_t = \mathbb{E} \left[ e^{-r(T-t)} \mathcal{P} (S_T)|\mathcal{F}_t \right].\tag{14}\]

In the following, we will assume without loss of generality that we conduct pricing at \(t=0\), and we will focus on the two most important path-independent claims, namely digital options (that can be either cash-or-nothing (CN) or asset-or-nothing (AN) options) and European (EUR) options. The payoffs of the digital options of maturity \(T\) and strike price \(K\) are given by \[\begin{cases} \mathcal{P}_{{\mathrm{AN}}}(S_T) = S_T \mathbb{1}_{\{S_T>K\}}, \\ \mathcal{P}_{{\mathrm{CN}}}(S_T) = \mathbb{1}_{\{S_T>K\}} , \end{cases}\] and the payoff of a European option of maturity \(T\) and strike \(K\) is: \[\label{eq:payoff95Eureopean} \mathcal{P}_{{\mathrm{EUR}}}(S_T) = \mathcal{P}_{{\mathrm{AN}}}(S_T) - K \mathcal{P}_{{\mathrm{CN}}}(S_T) = (S_T-K)\mathbb{1}_{\{S_T>K\}}.\tag{15}\]

For the two digital options, it is easy to see that the prices 14 can be carried out via the following integrals: \[\label{eq:digital-option-def} \begin{cases} \displaystyle C_{\mathrm{CN}}= e^{-rT}\int_{-k}^\infty f_{X_T}(x) \mathrm{d}x,\\ \displaystyle C_{\mathrm{AN}}= Ke^{k-rT}\int_{-k}^\infty e^xf_{X_T}(x) \mathrm{d}x, \end{cases}\tag{16}\] where \(f_{X_T}\) is the density function of \(X_T\) and \(k:= \ln(S_0/K) + (r-q+\zeta)T\) is the log-forward moneyness. Our strategy is the following: we first derive MB integral representations for the two digital calls with \(k<0\), then we express these integral as infinite series of residues. We extend the formulas for \(k>0\) thanks to the symmetry relations of the TS density ?? , and, last, we get the European call price by difference using 15 .

4.2 Digital options↩︎

As mentioned before, we give the MB representations of digital options and derive their associated series expansions.

4.2.1 Mellin-Barnes integral representations↩︎

We start by deriving the MB representation for the CN option in proposition 6; extension to the AN case will be handled in the subsequent proposition 7. The matrices and vectors \(U_{\mathrm{TS}}\), \(u_{\mathrm{TS}}\), \(D_{\mathrm{TS}}\), and \(d_{\mathrm{TS}}\) are the ones defined earlier in proposition. 3

4.2.2 Series expansions↩︎

Let us know express the MB integrals of propositions 6 and 7 as residue series. To that extent, given \(x_1<0\), \(x_2>0\) and an integer vector \(\boldsymbol{n} := (n_1,n_2,n_3)\in\mathbb{N}^3\), let us define the coefficients: \[\begin{align} \label{eq:coefficients-ts} \begin{cases} \begin{aligned} c_{\boldsymbol{n}}^{(1)}(x_1;x_2) := &\;(-\beta_-n_3)_{n_1} \frac{\Gamma(1-n_1+\beta_+n_2+\beta_-n_3)}{\Gamma(1+\beta_+n_2)\Gamma(-\beta_+n_2)} \\ &\times \underaccent{\bar}{\lambda}^{n_1} x_2^{-n_1+\beta_+n_2+\beta_-n_3} \gamma_{\text{inc}}(n_1-\beta_+n_2-\beta_-n_3,-x_1 x_2), \\ c_{\boldsymbol{n}}^{(2)}(x_1;x_2) := &\;(1+\beta_+n_2)_{n_1} \frac{\Gamma(-1-n_1-\beta_+n_2-\beta_-n_3)}{\Gamma(-\beta_+n_2)\Gamma(-\beta_-n_3)} \\ &\times \underaccent{\bar}{\lambda}^{1+n_1+\beta_+n_2+\beta_-n_3} x_2^{-1-n_1} \gamma_{\text{inc}}(1+n_1,-x_1x_2 ), \\ c_{\boldsymbol{n}}^{(3)}(x_2) := & \frac{\Gamma(n_1-\beta_+n_2-\beta_-n_3)}{\Gamma(1+n_1-\beta_+n_2)\Gamma(-\beta_-n_3)} \times \underaccent{\bar}{\lambda}^{-n_1+\beta_+ n_2+\beta_-n_3} x_2^{n_1}, \end{aligned} \end{cases} \end{align}\tag{22}\] where \((z)_n\) denotes the Pochhammer symbol (descending factorial) and \(\gamma_{\text{inc}}\) is the lower incomplete gamma function (see [56]). Furthermore, the coefficients 22 admit the special limiting cases: \[\label{eq:cn-coeff-boundary} \begin{cases} \displaystyle c_{ (n_1,0,0)}^{(1)}(x_1;x_2) := \frac{\beta_+}{\beta_-+\beta_+}\mathbb{1}_{\{n_1 = 0\}},\\ \displaystyle c_{(n_1,0,0)}^{(2)}(x_1;x_2) := 0, \\ \displaystyle c_{(0,0,0)}^{(3)}(x_2) := \frac{\beta_-}{\beta_-+\beta_+}. \end{cases}\tag{23}\] We finally define the overall coefficient \(c_{\boldsymbol{n}}\) as: \[\label{eq:coefficients-ts-sum} c_{\boldsymbol{n}}(x_1,x) := \frac{(-1)^{n_1}}{n_1!} \left( c_{\boldsymbol{n}}^{(1)}(x_1;x_2) +c_{\boldsymbol{n}}^{(2)}(x_1;x_2) \right) +c_{\boldsymbol{n}}^{(3)}(x_2).\tag{24}\]

In the following remark, we show that it is straightforward to extend proposition 8 to the case of a positive moneyness.

The AN option can be treated almost exactly the same way as the CN one, the difference being that we now express the MB integral ?? (featuring a \((\lambda_+-1)^{-s_4}\) term) instead of the MB integral ?? (featuring a \(\lambda_+^{-s_4}\) term). We obtain:

As for CN options, the case of positive moneyness easily handled, as shown in remark 5.

4.3 European options↩︎

As mentioned before, the European option price is simply obtained by differences of the AN (proposition 9) and CN (propsition 8) prices. We immediately obtain:

Last, the case of positive moneyness \(k>0\) is again straightforward to handle, thanks to remarks 5 and 4.

5 Particular cases↩︎

The pricing series ?? is valid under the most general TS processes; however, as we will see, it can be simplified in some particular cases that correspond to models that are popular in the literature. These simplifications are particularly remarkable in the cases of bilateral Gamma and Variance Gamma pricing.

5.1 KoBoL↩︎

The KoBoL process (see [7], [47]) corresponds to a TS process with the restriction \(\beta_+=\beta_-\); to obtain the corresponding European call price, it therefore suffices to take \(\beta_+=\beta_-\) in the coefficients 22 and the overall coefficient 24 , the latter being denoted by \(c_{\boldsymbol{n}}^{\text{KoBoL}}(k,x)\) in this case. Applying the pricing formula ?? (proposition 10) thus yields the KoBoL call price: \[\label{eur-series-KoBoL} C_{\mathrm{EUR}}^{\text{KoBoL}} = Ke^{-rT}\left[e^{k-\zeta T}-1 + e^{\gamma T} \sum_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{1}+n_2+n_3}}{n_{1}!n_2!n_3!} (a_+T)^{n_2}(a_-T)^{n_3} c_{{\mathrm{EUR}},\boldsymbol{n}}^{\text{KoBoL}}(k;\lambda_+) \right]\tag{26}\] where \[c_{{\mathrm{EUR}},\boldsymbol{n}}^{\text{KoBoL}}(k,\lambda_+) := e^kc_{\boldsymbol{n}}^{\text{KoBoL}}(k,\lambda_+-1) -c_{\boldsymbol{n}}^{\text{KoBoL}}(k,\lambda_+).\]

5.2 CGMY↩︎

The CGMY process of [6] is a KoBoL process with furthermore \(\alpha_+ = \alpha_- := C\) (and in this context the traditional usage is to write \(\beta_+=\beta_-:=Y)\). Interestingly, one can still use the KoBoL coefficients like in 26 ; the only difference actually lies in the \(a_\pm\) terms that are now equal and can therefore be grouped into the same power, yielding the CGMY call price: \[\label{eur-series-CGMY} C_{\mathrm{EUR}}^{\text{CGMY}} = Ke^{-rT}\left[e^{k-\zeta T}-1 + e^{\gamma T} \sum_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{1}}}{n_{1}!n_2!n_3!} (C \Gamma(-Y)T)^{n_2+n_3} c_{{\mathrm{EUR}},\boldsymbol{n}}^{\text{KoBoL}}(k;M) \right]\tag{27}\] where we have used the traditional notation \(\lambda_+=M\). As before, extension to in-the-money situation (\(k>0\) is straightforward for both KoBoL and CGMY processes, thanks to remarks 5 and 4.

5.3 Bilateral Gamma↩︎

The BG process of [5] corresponds to a TS process with the restriction \(\beta_+=\beta_-=0\). As such, it is therefore a particular case of a KoBoL process, but not of a CGMY process. A series pricing formula has already been derived in [5], however it is valid only for at-the-money options (\(k=0\)); a formula valid for all \(k\) has been derived in [45] however it involves three summation indices. In this section, we will obtain a much simpler formula, valid for all moneyness and involving only one summation index. For a \(\mathrm{BG}(\alpha_+,\lambda_+,\alpha_-,\lambda_-)\) process, let us define the following five quantities: \[\begin{cases} \displaystyle m_T := \frac{\lambda_+^{\alpha_+ T}\lambda_-^{\alpha_- T}}{\underline{\lambda}^{\underline{\alpha}T}\Gamma(\alpha_+ T)\Gamma(\alpha_- T)\Gamma(1-\alpha_+ T)},\\ \displaystyle \underline{\lambda} := \lambda_++\lambda_-,\\ \displaystyle \overline{\lambda} := \frac{1}{2}(\lambda_+-\lambda_-),\\ \displaystyle \underline{\alpha} := \frac{1}{2}(\alpha_++\alpha_-),\\ \displaystyle \overline{\alpha} := \frac{1}{2}(\alpha_+-\alpha_-). \end{cases}\] Moreover, the moneyness is now defined by \(k_{\mathrm{BG}}:=\log (S_0/K)+(r-q+\zeta_\mathrm{BG})T\) where \[\zeta_\mathrm{BG}:= -\ln \varphi_\mathrm{BG}(-\mathrm{i}) =\left(\frac{\lambda_+}{\lambda_+-1}\right)^{\alpha_+} \left(\frac{\lambda_-}{\lambda_-+1}\right)^{\alpha_-}.\] In the rest of this section, we will still denote \(k_\mathrm{BG}\) by \(k\), on order to preserve the simplicity of the notations.

In [5], the authors show that the density function \(f_\mathrm{BG}\) of the BG distribution can be written as: \[\forall x\in(0,+\infty),\quad f_\mathrm{BG}(x) =m_1\Gamma(1-\alpha_+)\Gamma(\alpha_-) x^{\underline{\alpha}-1} e^{-(\lambda_+-\lambda_-)x/2} W_{\overline{\alpha},\underline{\alpha}-1/2}(\underline{\lambda}x).\] Introducing the MB representation [56] for the Whittaker function yields \[\label{eq:MB-density-BG} f_\mathrm{BG}(x)= m_1 e^{-\lambda_+x} \int_{c+\mathrm{i}{\mathbb{R}}} \Gamma\left(\underline{\alpha}+s\right)\Gamma(1-\underline{\alpha}+s)\Gamma(-\overline{\alpha}-s)\underline{\lambda}^{-s}x^{-s-1+\underline{\alpha}} \frac{\mathrm{d}s}{\mathrm{i}2\pi}\tag{28}\] where \(c\in(-\underline{\alpha}\vee(\underline{\alpha}-1),-\overline{\alpha})\).

Before providing the series expansion for the option prices, let us first define for \(x_1<0\), \(x_2>0\), for an integer \(n\in\mathbb{N}\), and for a maturity \(T>0\), the coefficients: \[\label{eq:bg-series-coefficients} \begin{cases} \displaystyle c_n^{\mathrm{BG}, (1)}(x_1;x_2) := \Gamma(1 - 2 \underline{\alpha}T - n) \Gamma(\alpha_-T + n) \underline{\lambda}^{n + \underline{\alpha}T} x_2^{-n - 2 \underline{\alpha}} \gamma_\mathrm{inc} (n + 2 \underline{\alpha}T, -x_1 x_2), \\ \displaystyle c_n^{\mathrm{BG}, (2)}(x_1;x_2) := \Gamma(- \alpha_+T + 1 + n) \Gamma(2 \underline{\alpha}T - 1 - n) \underline{\lambda}^{- \underline{\alpha} T+ 1 + n} x_2^{-1 - n} \gamma_\mathrm{inc} (n+ 1, -x_1 x_2),\\ c_n^{\mathrm{BG}}(x_1;x_2) :=c_n^{\mathrm{BG}, (1)}(x_1;x_2)+c_n^{\mathrm{BG}, (2)}(x_1;x_2), \end{cases}\tag{29}\] and the function \(h\) defined by: \[\label{eq:h95function} h:x\in (-\underline{\lambda},\underline{\lambda})\mapsto \frac{\lambda_+^{\alpha_+T}\lambda_-^{\alpha_-T} \Gamma(2 \underline{\alpha}T)}{\underline{\lambda}^{2\underline{\alpha}T}\Gamma(1+\alpha_+T)\Gamma(\alpha_-T)} ~{_2}F_1 \left( 2\underline{\alpha}T,1 ; 1+\alpha_+T ; \frac{x}{\underline{\lambda}}\right),\tag{30}\] where \(_2F_1\) is the Gauss hypergeometric function (see [46] or any other classical monograph on special functions).

As before in the general TS case (see 15 ), the price of the European claim \(C_{\mathrm{EUR}}^{\mathrm{BG}}\) can be written as \(C_\text{Eur}^{\mathrm{BG}} = C_\text{AN}^{\mathrm{BG}} - KC_\text{CN}^{\mathrm{BG}}\).

5.4 Variance Gamma↩︎

The VG process of [4] is a BG process with the restriction \(\alpha_+=\alpha_-=\alpha\). The pricing formulas from proposition 12 remain valid, with \(\underline{\alpha}=\alpha\) and \(\overline{\alpha}=0\), and with the traditional VG parametrization: \[\label{VG95notations} \left\{ \begin{align} & \sigma^2 = \frac{2\alpha}{\lambda_+\lambda_-}, \\ & \theta = \frac{\alpha}{\lambda_+} - \frac{\alpha}{\lambda_-}, \\ & \nu = \frac{1}{\alpha}. \end{align} \right.\tag{33}\] When furthermore \(\lambda_+=\lambda_-=\lambda\), then \(\theta=0\) and the VG model is symmetric (see [48]); replacing in 30 and using Legendre’s duplication formula yields: \[h(x) = \frac{ \Gamma(\frac{1}{2} + \frac{T}{\nu})}{2\sqrt{\pi}\Gamma(1+\frac{T}{\nu})} ~{_2}F_1 \left( \frac{2T}{\nu},1 ; 1+\frac{T}{\nu} ; \frac{x}{\lambda}\right)\] and the ATM price ?? is: \[\begin{gather} \label{eq:atm-symvg} C_\mathrm{Eur}^\mathrm{sVG} = Ke^{-rT} \times \left( \left(1-\frac{\sigma^2\nu}{2}\right)^{-\frac{T}{\nu}}-1\right.\\ \left. + \frac{ \Gamma(1/2 + \frac{T}{\nu})}{2\sqrt{\pi}\Gamma(1+\frac{T}{\nu})} \left( {_2}F_1 \left( \frac{2T}{\nu},1 ; 1+\frac{T}{\nu} ; \frac{1}{2}\right) -{_2}F_1 \left( \frac{2T}{\nu},1 ; 1+\frac{T}{\nu} ; 1-\frac{\sigma\sqrt{\nu}}{2^{3/2}}\right)\right) \right). \end{gather}\tag{34}\]

6 Option pricing under one-sided tempered stable processes↩︎

Let us now focus on another important class of TS processes, namely one-sided TS processes, whose Lévy measure is either supported by the positive or negative real axis only. This class is actually important, and contains several models used for option pricing or credit risk purpose for instance in [61], [62]. We will see that for one-sided TS distributions, option prices can be expressed as series over a single index, which constitutes an improvement when compared to the series obtained in [44] where 3 summation indices are needed, and will allow for an accelerated pricing convergence.

Throughout this section, the underlying log-price process \((X_t)_{t\geqslant 0}\) is therefore assumed to be a \(\mathrm{TS}(\alpha,\beta,\lambda)\) process. We still keep the definition \(a:=-\alpha\Gamma(-\beta)\) and we adapt the definition of \(\gamma\) by setting now \(\gamma:=a\lambda^\beta\). The moneyness is now defined by \(k_{\mathrm{TS}_+}:=\log (S_0/K)+(r-q+\zeta_{\mathrm{TS}_+})T\) where \[\zeta_{\mathrm{TS}_+} := -\ln \varphi_{\mathrm{TS}_+} (-\mathrm{i}) =-a((\lambda-1)^\beta - \lambda^\beta)\] In the rest of this section, we will still denote \(k_{\mathrm{TS}_+}\) by \(k\), on order to preserve the simplicity of the notations.

6.1 Mellin-Barnes representation↩︎

As done for double-sided TS distributions, we first introduce Mellin-Barnes representations for digital option prices.

6.2 Series expansion↩︎

Let us start by defining the coefficient \(c_n^{+}(x_1;x_2)\) as follows: for a tuple \((x_1,x_2)\in (-\infty,0)\times (0,\infty)\) and \(n\in\mathbb{N}\), we let: \[c_n^+(x_1;x_2) := x_2^{\beta n}\Gamma(-n\beta,-x_1x_2).\]

As usual, it follows from 15 ) that the price of the European claim \(C_{\mathrm{EUR}}^{+}\) can be written as \(C_\text{Eur}{^+} = C_\text{AN}^{+} - KC_\text{CN}^{+}\) and we have: \[C_{{\mathrm{EUR}}}^+ = Ke^{(\gamma-r)T} \sum_{n\in\mathbb{N}} \frac{(-1)^{n}(aT)^{n}}{n!\Gamma(-n\beta)}c_{{\mathrm{EUR}},n}^+(k;\lambda)\] where \(c_{{\mathrm{EUR}},n}^+(k;\lambda) := e^kc_n^+(k;\lambda-1)-c_n^+(k;\lambda)\).

7 Numerical results↩︎

In this section, we provide numerical checks as well as a detailed performance study of the Mellin pricing series we obtained for the general double-side TS model, and also for specific cases that are of practical importance as they recover popular models (bilateral gamma, Variance Gamma, spectrally one-sided processes...); implementation is publicly available in the repository TS-Pricing (written in Python). As previously mentioned, one of the major advantages of TS distributions is the availability of an explicit characteristic function, allowing for the use of Fourier based methods for option pricing, e.g. Lewis-Lipton [37], Carr-Madan [27], Hilbert [60], COS [28], PROJ [29] and making these methods ideal challengers to assess the accuracy and efficiency of our Mellin pricing series.

Throughout this section, we will use the following set of parameters: \[\label{eq:numerical-params} \begin{cases} \lambda_+ =0.44,\\ \beta_+ = 0.1 + e/10,\\ \alpha_+ = 1.4,\\ \lambda_- = 0.35 ,\\ \beta_- = 0.5-\pi/100,\\ \alpha_- = 0.4.\\ \end{cases}\tag{35}\] The risk-free rate \(r\) and dividend yield \(q\) will be set to the values \(r=0.02\) and \(q=0.05\). All the experiences have been conducted on a personal laptop with CPU Intel Core i5-1021u-1.60GHz and operating system Ubuntu 24.04.

7.1 Numerical checks↩︎

The first set of tests that we conduct is aimed at checking the validity of the pricing formulas we have obtained. We consider the most general case (i.e., the double-sided TS model), and we compare the European call prices obtained via the Mellin series in proposition 10 to the prices obtained via the PROJ method; of course the choice of the PROJ method is arbitrary and any other Fourier related pricing method could be chosen as a reference, however we deliberately choose to compare to PROJ prices because the accuracy and efficiency of this method over other Fourier methods has been well evidenced in the literature across all vanilla and exotic options (see e.g. [29], [30], [32]). We perform these comparisons for various spot price \(S_0\) and maturity \(T\), and display the results in figure 2: it is clear that the series prices perfectly match the PROJ prices across all moneyness. Let us note that we choose to perform the summations up to order \(n_1=n_2=n_3=80\), in the series prices but, actually and as will be demonstrated in the next subsections, a truncation to a far lower order would ensure excellent numerical accuracy in most situations.

7.2 Performance study↩︎

In this section, we want to demonstrate the practical tractability and efficiency of the series formula we obtained, by displaying the computation time needed to attain a predetermined precision, and by comparing it with the PROJ computation time. We perform this study for the most general double-sided TS model and some of its particular cases (BG and one-sided TS), demonstrating a very competitive convergence speed when compared to state-of-the-art Fourier pricing; as we will see, the convergence is particularly accelerated in the BG and one-sided TS cases.

First, we note that a major advantage of our pricing method is the total absence of a parameter to set, while such a parameter (or even several) is in general mandatory to implement Fourier pricing methods and can be delicate to select. This parameter can be a control parameter, conditioned by square integrability in the case of Carr-Madan pricing (it is known as a damping parameter in this case and is often denoted by \(\alpha\)), or be related to truncation intervals (for instance, in the case of the PROJ pricing it is related to the integration support size of the projection component). Even if some heuristic rules can be used to choose this parameter, notably some moment-based rules, these heuristics fail when it comes to high precision pricing, and manual increase of this parameter must therefore be performed, as illustrated in tables 1, 2 and 3. In particular, in the double TS case (table 1), it is clear that neither the “default" value of \(\alpha\) computed with the heuristic, nor several other attempts for this paremeter (\(\alpha\in\{25,30,40\}\)), allow to achieve a \(10^{-8}\) precision, even by increasing the discretization parameter \(N_{\text{PROJ}}\) (up to \(N_{\text{PROJ}}=2^{20}\)) and the number of functions in the Riesz basis - actually a manual increase to \(\alpha=50\) is necessary to attain this precision; in contrast, the series formula ?? allows for a direct computation of prices without having to worry about such parametrization, which of course is a very desirable feature in practice. We note, however, that, once the appropriate truncation parameter is chosen, the PROJ method remains approximately 100 times faster than the Mellin pricing formula ?? (due to the presence of 3 summation indices); the series display nonetheless a very satisfactory performance, with only 0.09 second needed to reach a \(10^{-9}\) precision.

Under bilateral Gamma and one-sided tempered stable processes, the pricing series for the European call are obtained from the digital option prices provided in propositions 12 and 14 respectively, and are performed over a single index only (instead of three summation indices in the general TS case). This allows for considerably faster computation, as shown in tables 2 and 3. For example, in the bilateral Gamma case, if high precision (\(10^{-8}\sim10^{-9}\)) is desired, the Mellin pricing series is about 10 times faster than PROJ (even after choosing a hyperparameter \(\alpha=50\)); a similar overperformance is observed for the one-sided TS model.

This accelerated convergence behavior in the bilateral Gamma and one-sided TS cases is also clearly visible in figure 3, where relative error (in dotted lines) and computation time (in plain lines) are displayed as functions of the maximal summation index \(N\), for both these models as well as for the most general double-sided TS model. In terms of error, we can see that the one-sided pricing formula reaches a \(10^{-10}\) precision after 20 summation terms only, and roughly reaches machine precision after 25 terms, with a total computation time of \(10^{-5}\sim 10^{-4}\) seconds only (in line with the results of table 3). For the bilateral Gamma and double-sided TS cases, the precision increases progressively and, actually, almost linearly as a function of \(N\), which constitutes a surprising and very interesting behavior, as it allows to reliably control the accuracy (at least from an empirical point of view) of the pricing results; we also note that, for these two models, we can reach a \(10^{-10}\) precision for \(N=80\) approximately, but the overall computation time need to reach this precision remains around \(10^{-4}\) seconds for the bilateral Gamma model (which, again, is extremely fast), while it is of order \(10^{-1}\sim1\) seconds for the general doublie-sided case. Of course and as already mentioned, this is because only \(N\) terms are needed to perform the bilateral Gamma summations in proposition propositions 12, while \(N^3\) terms are needed in the double-sided TS series ?? .

7.3 Comparison to other Fourier methods↩︎

Finally, we perform a final check by providing a short numerical comparison to other standard Fourier based pricing techniques, namely the PROJ [29], Lewis [37], Carr-Madan [27], Gil-Pelaez [63] and Hilbert [60] methods. We perform this check for several classical jump models covered by the TS family:

1. TS, KoBoL and CGMY, that pertain to the double-sided TS family, and whose Mellin pricing series is given by ?? ;

2. BG and VG, that are within the bilateral Gamma family (itself a subfamily of the TS class with \(\beta_\pm=0\)) and whose Mellin pricing series is given by proposition 12;

3. TS\(^+\), that is, the spectrally positive TS distributions, that are also a special case of TS distributions (with a Lévy measure supported by the positive real axis, i.e. \(\alpha_-=\lambda_-=\beta_-=0\)) and whose Mellin pricing series is given by proposition 14.

The results are summarized in table 4 and feature very good agreement; without suprise, Lewis, PROJ and Mellin display best precision.

8 Conclusion↩︎

In this paper, we have proved Mellin-Barnes representations for TS density functions and derived their associated series expansions. We have also provided Mellin-Barnes representations for the prices of digital and European options when the log-price process is a TS process, and we were able to express these integrals as series expansions whose terms are straightforward to compute. These formulas become remarkably simple when considering particular cases such as bilateral gamma or one-sided TS processes. Moreover, our technique has the major advantage of being hyperparameter free.

We have also presented a detailed numerical analysis of the pricing series expansions, including performance studies and comparison with other approaches. This analysis demonstrates that our method is competitive with state-of-the-art pricing techniques. All numerical implementations can be found in our public repository TS-Pricing (written in Python).

Future work should consider adapting our methodology to TS processes with \(\beta_\pm \in(1,2)\), these processes having a different characteristic function. Other extensions should include derive similar series expansions for additive processes and their particular cases (e.g. Sato processes).

9 Computational details for the proofs of propositions 8 and 9↩︎

In this appendix, we provide detailed computations of the residues associated to the poles (P1), (P2) and (P4) arising in the proofs of proposition 8 and 9.

9.1 Residues associated to (P1)↩︎

Let \(S_1\) define the series of residues induced by the singularities in the subset (P1). We have: \[\begin{gather} S_1 := e^{(\gamma-r)T} \sum_{\boldsymbol{n}\in\mathbb{N}^4} \frac{(-1)^{n_{1+4}}}{n_{1,4}!} \frac{ \Gamma(1-n_1+\beta_+n_2+\beta_-n_3) \Gamma(n_1-\beta_-n_3) }{ \Gamma(1+\beta_+n_2) \Gamma(-\beta_+n_2) \Gamma(-\beta_-n_3) (-n_1+\beta_+n_2+\beta_-n_3-n_4) } \\ \times (a_+T)^{n_2}(a_-T)^{n_3}\lambda_+^{n_4} \underaccent{\bar}{\lambda}^{n_1} (-k)^{n_1-\beta_+n_2-\beta_-n_3+n_4}. \end{gather}\]

Some algebra yields: \[\begin{align} S_1 &= \sum_{\boldsymbol{n}\in\mathbb{N}^4} \frac{(-1)^{n_{1+4}}}{n_{1,4}!} \frac{(-\beta_-n_3)_{n_1}}{(-n_1+\beta_+n_2 + \beta_-n_3 - n_4)}\\ &\quad\quad\times\frac{\Gamma(1-n_1+\beta_+n_2+\beta_-n_3)}{\Gamma(1+\beta_+n_2)\Gamma(-\beta_+n_2)} (a_+T)^{n_2}(a_-T)^{n_3} \underaccent{\bar}{\lambda}^{n_1} (-k)^{n_1-\beta_+n_2-\beta_-n_3+n_4}\lambda_+^{n_4}\\ &= -\sum\limits_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{1+3}}}{n_{1,3}!} (-\beta_-n_3)_{n_1}\frac{\Gamma(1-n_1+\beta_+n_2+\beta_-n_3)}{\Gamma(1+\beta_+n_2)\Gamma(-\beta_+n_2)} (a_+T)^{n_2}(a_-T)^{n_3} \underaccent{\bar}{\lambda}^{n_1} \\ &\quad\quad\times\lambda_+^{-n_1+\beta_+n_2+\beta_+n_3}\sum_{n_4\in\mathbb{N}}\frac{(-1)^{n_4}}{n_4!(n_1-\beta_+n_2 - \beta_-n_3 + n_4)}(-\lambda_+k)^{n_1-\beta_+n_2-\beta_-n_3+n_4}\\ &= -\sum\limits_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{1+3}}}{n_{1,3}!} (-\beta_-n_3)_{n_1}\frac{\Gamma(1-n_1+\beta_+n_2+\beta_-n_3)}{\Gamma(1+\beta_+n_2)\Gamma(-\beta_+n_2)} (a_+T)^{n_2}(a_-T)^{n_3} \underaccent{\bar}{\lambda}^{n_1} \\ &\quad\quad\times\lambda_+^{-n_1+\beta_+n_2+\beta_-n_3}\gamma_{\text{inc}}(n_1-\beta_+n_2-\beta_-n_3, -\lambda_+k) \end{align}\] where the last line is obtained by recognizing the series expansion for the incomplete gamma function (see [56]). It comes that: \[S_1 = -e^{(\gamma-r)T}\sum_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{1+3}}}{n_{1,3}!} (a_+T)^{n_2}(a_-T)^{n_3} c_{\boldsymbol{n}}^{(1)}(k;\lambda_+).\]

9.2 Residues associated to (P2)↩︎

Let \(S_2\) define the series of residues induced by the singularities in the subset (P2). We have: \[\begin{gather} S_2 := \sum_{\boldsymbol{n}\in\mathbb{N}^4} e^{(\gamma-r)T}\frac{(-1)^{n_{1+4}}}{n_{1,4}!} \frac{ \Gamma(-1-n_1-\beta_+n_2-\beta_-n_3) \Gamma(1+n_1+\beta_+n_2) }{ \Gamma(1+\beta_+n_2) \Gamma(-\beta_+n_2) \Gamma(-\beta_-n_3) (-1-n_1-n_4) } \\ \times (a_+T)^{n_2}(a_-T)^{n_3}\lambda_+^{n_4} \underaccent{\bar}{\lambda}^{1+n_1+\beta_+n_2+\beta_-n_3} (-k)^{1+n_1+n_4}. \end{gather}\]

Some algebra yields: \[\begin{align} S_2 &= \sum_{\boldsymbol{n}\in\mathbb{N}^4}\frac{(-1)^{n_{1+4}}}{n_{1,4}!} \frac{(1+\beta_+n_2)_{n_1}}{-1-n_1-n_4} \frac{\Gamma(-1-n_1-\beta_+n_2-\beta_-n_3)}{\Gamma(-\beta_+ n_2)\Gamma(-\beta_-n_3)} (a_+T)^{n_2}(a_-T)^{n_3}\\ &\quad\quad\times\underaccent{\bar}{\lambda}^{1+n_1+\beta_+n_2+\beta_-n_3} (-k)^{1+n_1+n_4} \lambda_+^{n_4}\\ & = -\sum_{\boldsymbol{n}\in\mathbb{N}^3}\frac{(-1)^{n_{1+3}}}{n_{1,3}!} (1+\beta_+n_2)_{n_1} \frac{\Gamma(-1-n_1-\beta_+n_2-\beta_-n_3)}{\Gamma(-\beta_+ n_2)\Gamma(-\beta_-n_3)} (a_+T)^{n_2}(a_-T)^{n_3}\\ &\quad\quad\times\underaccent{\bar}{\lambda}^{1+n_1+\beta_+n_2+\beta_-n_3} \sum_{n_4\in\mathbb{N}}\frac{(-1)^{n_{4}}}{n_{4}!}(-k)^{1+n_1+n_4}\frac{1}{1+n_1+n_4} \lambda_+^{n_4}\\ &= -\sum_{\boldsymbol{n}\in\mathbb{N}^3}\frac{(-1)^{n_{1+3}}}{n_{1,3}!} (1+\beta_+n_2)_{n_1} \frac{\Gamma(-1-n_1-\beta_+n_2-\beta_-n_3)}{\Gamma(-\beta_+ n_2)\Gamma(-\beta_-n_3)} (a_+T)^{n_2}(a_-T)^{n_3}\\ &\quad\quad\times\underaccent{\bar}{\lambda}^{1+n_1+\beta_+n_2+\beta_-n_3} \lambda_+^{-1-n_1} \sum_{n_4\in \mathbb{N}}\frac{(-1)^{n_{4}}}{n_{4}!}(-\lambda k)^{1+n_1+n_4}\frac{1}{1+n_1+n_4}\\ &= -\sum_{\boldsymbol{n}\in\mathbb{N}^3}\frac{(-1)^{n_{1+3}}}{n_{1,3}!} (1+\beta_+n_2)_{n_1} \frac{\Gamma(-1-n_1-\beta_+n_2-\beta_-n_3)}{\Gamma(-\beta_+ n_2)\Gamma(-\beta_-n_3)} (a_+T)^{n_2}(a_-T)^{n_3}\\ &\quad\quad\times\underaccent{\bar}{\lambda}^{1+n_1+\beta_+n_2+\beta_-n_3} \lambda_+^{-1-n_1}\gamma_{\text{inc}}(1+n_1,-\lambda_+ k) \end{align}\] where the last line is obtained by recognizing the series expansion for the lower incomplete gamma function (see [56]). It comes that: \[S_2 = -e^{(\gamma-r)T}\sum_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{1+3}}}{n_{1,3}!} (a_+T)^{n_2}(a_-T)^{n_3} c_{\boldsymbol{n}}^{(2)}(k;\lambda_+).\]

9.3 Residues associated to (P4)↩︎

Let \(S_4\) define the series of residues induced by the singularities in the subset (P4). We have: \[\begin{gather} \label{S495residue} S_4 := \sum_{\boldsymbol{n}\in\mathbb{N}^3} e^{(\gamma-r)T}\frac{(-1)^{n_{1+3}}}{n_{1,3}!} \frac{ \Gamma(n_1-\beta_-n_3) \Gamma(1+n_1+\beta_+n_2) \Gamma(-n_1-\beta_+n_2) }{ \Gamma(1+\beta_+n_2) \Gamma(-\beta_+n_2) \Gamma(-\beta_-n_3) }\\ \times (a_+T)^{n_2}(a_-T)^{n_3}\lambda_+^{n_1+\beta_+n_2} \underaccent{\bar}{\lambda}^{-n_1+\beta_-n_3}. \end{gather}\tag{36}\] We first notice that this can be rewritten with Pochhamer’s symbols: \[S_4 := e^{(\gamma-r)T}\sum_{\boldsymbol{n}\in\mathbb{N}^3} \frac{(-1)^{n_{2+3}}}{n_{1,3}!} (-\beta_-n_3)_{n_1} (a_+T)^{n_2}(a_-T)^{n_3}\lambda_+^{n_1+\beta_+n_2} \underaccent{\bar}{\lambda}^{-n_1+\beta_-n_3}.\] since \(\Gamma(-n_1-\beta_+n_2)/\Gamma(-\beta_+n_2) = (-1)^{n_1}/(1+\beta_+n_2)_{n_1}\). By first summing over \(n_1\), one recognize the hypergeometric series \(_1F_0(-\beta_+n_3;;\lambda_+/\underline{\lambda})\), whose sum is known in closed form (\(_1F_0a;;z)=(1-z)^{-a}\)) and we therefore have: \[S_4 := e^{(\gamma-r)T}\sum_{(n_2,n_3)\in\mathbb{N}^2} \frac{(-1)^{n_{2+3}}}{n_{2,3}!} (a_+T)^{n_2}(a_-T)^{n_3}\lambda_+^{\beta_+n_2} \underaccent{\bar}{\lambda}^{\beta_-n_3} (1-\lambda_+/\underaccent{\bar}{\lambda})^{\beta_-n_3}\] Noticing that \((1-\lambda_+/\underaccent{\bar}{\lambda})^{\beta_+n_3} = \lambda_-^{\beta_-n_3}\) (definition 8 ) and that the two remaining series are simply exponential expansions, we have: \[S_4 := e^{(\gamma-r)T}e^{-a_-T\lambda_-^{\beta_-}-a_+T\lambda_+^{\beta_+}} =e^{-rT} ,\] where the last equality stems from the definition of \(\gamma\) in 8 .

9.4 Residues associated to (P4) in the AN case↩︎

In the AN case, the residues 36 induced by the (P4) subset become: \[\begin{gather} S_4^{\text{AN}} := \sum_{\boldsymbol{n}\in\mathbb{N}^3} Ke^{k+(\gamma-r)T}\frac{(-1)^{n_{1+3}}}{n_{1,3}!} \frac{ \Gamma(n_1-\beta_-n_3) \Gamma(1+n_1+\beta_+n_2) \Gamma(-n_1-\beta_+n_2) }{ \Gamma(1+\beta_+n_2) \Gamma(-\beta_+n_2) \Gamma(-\beta_-n_3) }\\ \times (a_+T)^{n_2}(a_-T)^{n_3}(\lambda_+-1)^{n_1+\beta_+n_2} \underaccent{\bar}{\lambda}^{-n_1+\beta_-n_3}. \end{gather}\] Performing the same computations as in section 9.3 yields: \[S_4^{\text{AN}} := Ke^{k+(\gamma-r)T}\sum_{(n_2,n_3)\in\mathbb{N}^2} \frac{(-1)^{n_{2+3}}}{n_{2,3}!} (a_+T)^{n_2}(a_-T)^{n_3}(\lambda_+-1)^{\beta_+n_2} (\lambda_-+1)^{\beta_-n_3}.\] Recognizing the exponential series yields: \[S_4^{\text{AN}} = Ke^{k+(\gamma-r)T}e^{-a_+T(\lambda_+-1)^{\beta_+}-a_-T(\lambda_-+1)^{\beta_-}} = Ke^{k+(\zeta-r)T}\] where the last equatlity follows from the definition of \(\gamma\) in 8 and of the martingale adjustment \(\zeta\) in 13 .

10 ATM option price under bilateral gamma process↩︎

In [5], the authors give a closed formula for the ATM digital options in the BG model. Let us prove that the CN component in our ATM pricing formula ?? is the same the one obtained by [5] (same proof can be made for the AN option, and therefore for the European option). To this aim, we consider, as in [5], the risk-free rate and dividend yield to be \(r=q=0\). With these settings, the formula obtained by [5] for the CN call is: \[\frac{\lambda_+^{\alpha_+}\lambda_-^{\alpha_-}\Gamma(\alpha_++\alpha_-) }{\Gamma(\alpha_+)\Gamma(1+\alpha_-)\lambda_+^{\alpha_++\alpha_-}} ~_2F_1(\alpha_++\alpha_-,\alpha_-;1+\alpha_-,-\lambda_-/\lambda_+).\] Introducing \(\Lambda:=\lambda_-/\lambda_+\), this equation can be simplified in: \[\frac{\Lambda^{\alpha_-}\Gamma(\alpha_++\alpha_-) ~_2F_1(\alpha_++\alpha_-,\alpha_-;1+\alpha_-,-\Lambda)}{\Gamma(\alpha_+)\Gamma(1+\alpha_-)}\] and applying [56] we have: \[\begin{gather} \label{eq:annex-interm-2f1} _2F_1(\alpha_++\alpha_-,\alpha_-;1+\alpha_-,-\Lambda) = \frac{\Gamma(\alpha_+)\Gamma(1+\alpha_-)}{\Gamma(\alpha_++\alpha_-)(1+\Lambda)^{\alpha_-}} ~_2F_1(\alpha_-,1-\alpha_+;1-\alpha_+;(1+\Lambda)^{-1})\\ +\frac{\Gamma(-\alpha_+)\Gamma(1+\alpha_-)}{\Gamma(\alpha_-)\Gamma(1-\alpha_+)(1+\Lambda)^{\alpha_-+\alpha_+}} ~_2F_1(\alpha_++\alpha_-,1;1+\alpha_+;(1+\Lambda)^{-1}). \end{gather}\tag{37}\] Substituting the identity: \[_2F_1(\alpha_-,1-\alpha_+;1-\alpha_+;(1+\Lambda)^{-1}) = ~_1F_0(\alpha_-,(1+\Lambda)^{-1}) = \left(\frac{\Lambda}{1+\Lambda}\right)^{-\alpha_-}\] in 37 yields: \[\begin{gather} \frac{\Lambda^{\alpha_-}\Gamma(\alpha_+ + \alpha_-)}{\Gamma(1+\alpha_-)\Gamma(\alpha_+)}~_2F_1(\alpha_++\alpha_-,\alpha_-;1+\alpha_-,-\Lambda) = 1\\ +\frac{\Lambda^{\alpha_-}\Gamma(-\alpha_+)\Gamma(\alpha_++\alpha_-)}{\Gamma(\alpha_-)\Gamma(1-\alpha_+)\Gamma(\alpha_+)(1+\Lambda)^{\alpha_-+\alpha_+}} ~_2F_1(\alpha_++\alpha_-,1;1+\alpha_+;(1+\Lambda)^{-1}). \end{gather}\] Noticing that \(\Gamma(-\alpha_+)/\Gamma(1-\alpha_+) = -\alpha_+\) and \(\alpha_+\Gamma(\alpha_+) = \Gamma(1+\alpha_+)\), we have: \[\begin{gather} \frac{\Lambda^{\alpha_-}\Gamma(\alpha_+ + \alpha_-)}{\Gamma(1+\alpha_-)\Gamma(\alpha_+)}~_2F_1(\alpha_++\alpha_-,\alpha_-;1+\alpha_-,-\Lambda) = 1\\ -\frac{\Lambda^{\alpha_-}\Gamma(\alpha_++\alpha_-)}{\Gamma(\alpha_-)\Gamma(1+\alpha_+)(1+\Lambda)^{\alpha_-+\alpha_+}} ~_2F_1(\alpha_++\alpha_-,1;1+\alpha_+;(1+\Lambda)^{-1}). \end{gather}\] With the following relation: \[\frac{\Lambda^{\alpha_-}}{(1+\Lambda)^{\alpha_++\alpha_-}} = \frac{\lambda_-^{\alpha_-}\lambda_+^{\alpha_+}}{(\lambda_++\lambda_-)^{\alpha_++\alpha_-}},\] the desired result follows.

References↩︎