Residue Sum Formula for Pricing Options under the Variance Gamma Model

We present and prove a triple sum series formula for the European call option price in a market model where the underlying asset price is driven by a Variance Gamma process. In order to obtain this formula, we present some concepts and properties of multidimensional complex analysis, with particular emphasis on the multidimensional Jordan Lemma and the application of residue calculus to a Mellin–Barnes integral representation in C3, for the call option price. Moreover, we derive triple sum series formulas for some of the Greeks associated to the call option and we discuss the numerical accuracy and convergence of the main pricing formula.


Introduction
The pricing of financial derivatives, such as options, is one of the pivotal tasks of mathematical finance, yet it can be an arduous task to develop a model that is consistent with the empirical evidence, soluble, and where its numerical estimation is neither erroneous nor time consuming. One of the first attempts to solve this quandary was the Gaussian model that was first introduced by Fischer Black and Myron Scholes in [1] and later expanded by Robert Merton in [2], aptly named the Black-Scholes model, where the stochastic process driving the underlying asset price is modeled by a geometric Brownian motion. Its simplicity and the admission of a close formula for the option price are the main reasons why, to this day, it remains the most frequently used model by market practitioners. Still, the model fails to account for sudden price drops or increases that can be expressed as discontinuous price jumps; moreover, it assumes the volatility to remain constant for changes in relation to the strike price and time to maturity, contrary to the evidence derived from empirical data and, furthermore, the distributions of asset returns have been shown to be negatively skewed and exhibit fat-tails that are not captured by the symmetric Gaussian model.
The Black-Scholes model has been generalized in many ways. Let us mention two of the most important classes of models used in these generalizations: stochastic volatility models and jump models (see [3]). The model that we consider in this paper is a particular jump model, which assumes that the underlying asset price dynamics are described by a Lévy Process, namely the Variance Gamma process, which was first proposed by Dilip Madan and Eugene Seneta in [4]. The descriptive power of models that are based on Lévy processes for accurately portraying financial markets (not displaying the aforementioned problems that are present in the Black-Scholes model) has been known since the works of Benoît Mandelbrot [5] and Eugene Fama [3], and they have been gaining traction in recent decades with the advent of technology and computer development. Yet, the Black-Scholes model remains mostly ubiquitous. The main reason for this state of affairs is that pricing data presented in [23], we test the accuracy of the Variance Gamma formula and its Greeks. The last section is dedicated to some concluding remarks.

Multidimensional Residue Calculus
We start by enumerating, without proving, some of the definitions and results of multidimensional complex analysis. A more in-depth theoretical introduction to multidimensional complex analysis can be found in the textbooks [24,25].

Definition 1 (Grothendieck Residue).
Let h and f i , for any index i ∈ {1, . . ., n}, be functions in C n , where h is holomorphic. Consider the meromorphic differential n-form which has the singularities D j = {z ∈ C : f j (z) = 0}, such that the intersection n j=1 D j is discrete. The Grothendieck residue on a singularity a ∈ n j=1 D j is defined as where C a = {z ∈ U a : f j (z) = , j = 1, . . ., n} is a cycle in a small neighborhood U a of the singularity a with the orientation d(arg f 1 ) ∧ . . . ∧ d(arg f n ) ≥ 0.
Before proceeding we will formalize the concept of a multidimensional polyhedron in C n . The two following definitions will underpin most of the theorems from Sections 2 and 3.
Definition 2 (Polyhedron). Consider a proper (the inverse images of a compact set are compact) holomorphic map g : C n → G where G = G 1 × . . . × G n is a domain (connected open subset of a finite-dimensional vector space) where, for each j = 1, . . ., n, G j ⊂ C is a domain with piecewise smooth boundary. We define a polyhedron Π as the inverse image Π := g −1 (G), (3) and for a multi-index K = {k 1 , . . ., k p } ⊂ {1, . . ., n} we define the polyhedron's faces as σ K := {z : g k (z) ∈ ∂G k for k ∈ K, g j (z) ∈ G j for j ∈ K C }. (4) Definition 3 (Compatible divisors). Consider the polyhedron Π and the family of divisors {D i } i∈{1,...,n} , they are said to be compatible if for any i ∈ {1, . . ., n} we get Analogously to the one dimensional integral on the real axis, we may want to compute an integral σ ω where ω is the meromorphic form (1) and σ is the boundary of a polyhedron Π. For an unbounded polyhedron we need the integrand to vanish as it goes to infinity. To achieve this goal let us introduce the auxiliary functions ρ j = f j 2 f 2 , for any j ∈ {1, . . ., n}, where f 2 = | f 1 | 2 + . . . + | f n | 2 . Using the functions (6) we define the differential (n, p − 1)forms as where J = {j 1 , . . ., j s } ⊂ {1, . . ., n}, for 1 < s ≤ n, is a multi-index, (j, J) is the position of j in set J and ∂ρ J [j] = ∂ρ j 1 ∧ . . . [j]. . . ∧ ∂ j s . We now define the multidimensional condition analogous to lim R→∞ S R f (z)dz = 0.
Definition 4 (Jordan condition). Consider the sphere S R = {z ∈ σ : z = R}, where σ = σ 12...n is the boundary of the polyhedron Π. A differential form ξ J satisfies the Jordan condition on face σ J o , with J o = {1, . . ., n} \ J, if exists a sequence of positive real numbers R k that goes to infinity, such that Note that, for n = 1, there exists only one form ξ = ω and thus the condition (8) corresponds to the unidimensional condition lim z→∞ S R f (z)dz = 0. For the multidimensional case, consider the set N = {z ∈ C n : f (z) = 0} = n i=1 D i . Thus, we have the following theorem, which allows us to compute the integral of the meromorphic form (1) on the boundary of the polyhedron as a sum of residues, and it is the main theoretical tool used in the proof of Theorem 2.

Theorem 1 (The Jordan Lemma).
Let ω be a meromorphic form with the polar divisors {D i } i∈{1,...,n} compatible with polyhedron Π. If, for every multi-index J, the differential form ξ J satisfies the Jordan condition on σ J o , then we get The Jordan's Lemma proof can be found in [17].

One-Dimensional Mellin-Barnes Integral
We will start by presenting, without proof, some basic concepts and properties for the one-dimensional Mellin-Barnes integral. For a more in depth look at Fourier, Laplace, and Mellin Transforms, as well as their corresponding properties, we recommend the book [26].
Definition 5 (Mellin-Barnes Integral). The Mellin-Barnes Integral is given by a ratio of products of Gamma functions of linear arguments where its characteristic quantity, ∆, is defined by Before proceeding, let us state the Stirling's approximation of the gamma function (see [27]) Equipped with this expression, one can easily see how ∆ governs the behavior of the ratio of Gammas as |z| → ∞, and, therefore, which residues, of the singularities to the left or right of the (z) = 0 strip, one will sum to compute the integral (10): For example, one can express an exponential term e x as the Mellin-Barnes integral where γ > 0.

Three-Dimensional Mellin-Barnes Integral
In this subsection, similarly to what we did previously for the unidimensional case, we will present the triple Mellin-Barnes integral and deduce its formula as a sum of residues, which will be crucial when proving the main result of this paper. In order to achieve this, let us first consider its integral form Henceforth, for brevity, we will denote the three-form integrand of (15) by ω. Its zeroes will be the complex planes L ν j = {(z 1 , z 2 , z 3 ) ∈ C 3 : a j1 z 1 + a j2 z 2 + a j3 z 3 + b j = −ν}, for any ν ∈ N and j ∈ {1, . . ., m}, which represent each singularity given by the gamma functions present in the numerator of the form ω. We will also denote the vectors a j :=   a j1 a j2 a j3   , c k :=   a k1 a k2 a k3   and, most importantly, define the characteristic vector as Suppose that ∆ is a non-zero vector. In this case, we can define the plane P ∆ where its real part intersects the point γ and has ∆ as its normal vector, i.e., P ∆ := {z ∈ C 3 : ( ∆, z ) = ( ∆, γ )}, and thereupon we can define the admissible-polyhedra, Π ∆ , as the real volume "below" P ∆ , i.e., Π ∆ := {z ∈ C 3 : ( ∆, z )) ≤ ( ∆, γ )}.
From (17), we can ascertain that its only vertex is σ 1,2,3 = γ and that n 1 , n 2 , and n 3 are the normal vectors of the faces σ 1 , σ 2 , and σ 3 of the polyhedron: If the polyhedron was providently constructed, then we can balkanize the singularities, L ν j , into three distinct sets, such that they are compatible with the polyhedron, i.e.: Theorem 2 (Residue formula for the Triple Mellin-Barnes integral). Let ω be the three-form integrand of (15) with the characteristic vector ∆ = 0 and divisors D 1 , D 2 , and D 3 , as defined in (21), being compatible with the admissible polyhedron Π ⊂ Π ∆ . Subsequently, the sum formula holds: where the series on the right-hand side converges absolutely for any t ∈ U, for U, defined as, U = {t ∈ (C \ 0) 3 : |arg t 1 | < π, |arg t 2 | < π, |arg t 3 | < π, arg t < (π/2)α} (23) and α is a constant that is given by, where S 1 = {y ∈ R 3 : |y| = 1} is the unit sphere in R 3 .
For a proof of Theorem 2, please see Appendix A.

Option Pricing Driven by a Variance Gamma Process
Equipped with the results of the previous section, we will now center our attention on deducing the main result of this paper, the formula for the price of a European call option under the Variance Gamma model. Similar to the derivation presented in [20], we will arrive at this result in two steps. First, we will derive the Mellin-Barnes representation for the aforementioned call option price and, secondly, we will use residue calculus to derive the triple sum series formula.

Mellin-Barnes Representation for a Call Option
Before introducing the call option price, let us recall that we can define the Variance Gamma process as the difference between two independent gamma processes (see [6]), where τ is the time variable, G 1 τ ∼ Gamma(Cτ, 1/M) and G 2 τ ∼ Gamma(Cτ, 1/G) with parameters C > 0, M > 0, and G > 0. The probability density function of a gamma process G ∼ Gamma(α, β) is given by Under these definitions, one can, by direct computation, easily arrive at the characteristic function for the variance gamma process X VG (τ; C, G, M), where u ∈ R is the main variable of the characteristic function and τ is the time variable.
Let us now briefly discuss the standard problem of option pricing. For more details about this problem, we refer to [28]. Consider that we have an underlying risky asset with price process denoted by the stochastic process {S t } t∈[0,T] , where T denotes the maturity of the European option that we want to evaluate. In our market model, we assume that the continuously compounded risk-free interest rate is denoted by the parameter r and q represents the dividend yield. In the standard Black-Scholes model, it is well known that the market is arbitrage-free and complete. Moreover, we know from the fundamental theorems of asset pricing that an equivalent martingale measure or risk-neutral measure Q exists if and only if the market is arbitrage-free, and the measure Q is unique if and only if the market is complete. Therefore, under the standard Black-Scholes model, the equivalent martingale measure Q exists and it is unique. The measure Q is equivalent to the real (physical) probability measure P and the process S t e −(r−q)t is a martingale under Q. The underlying risky asset price at time T, under the measure Q, is given by where τ = T − t and B Q τ is a standard Brownian motion under measure Q. When considering Lévy Processes (like the Variance Gamma process), the market is, in general, incomplete and many equivalent martingale measures exist. Using the mean correcting martingale measure is one standard way to choose just one equivalent martingale measure. Under this particular measure, in the market model that is driven by the Variance Gamma process, we have that the risky asset price at time T is given by where µ := log(φ VG (−i, 1; C, G, M)) is the mean correcting parameter. For more details on the mean correcting martingale measure and the mean correcting parameter, see Chapter 6 of [23]. Let us recall that a European call option with maturity T and strike price K > 0 is a contract that gives the buyer of the contract the right, but not the obligation, to buy the underlying risky asset at time T by the price K. The payoff of the call option is given by (S T − K) + . Using the well-known risk-neutral valuation formula for option pricing, we can write the price at time t of a European call option with maturity T and strike price K > 0 as the discounted expected value where the variables are the underlying risky asset price at time t (denoted by S = S t ), the strike price of the option K, the risk-free interest rate r, the dividend yield q, and the time until maturity τ = T − t.
Using the notation introduced in the previous paragraphs, we can now state and derive the Mellin-Barnes integral representation in C 3 for the European call option price under the Variance Gamma model.
). Consider the polyhedra P 1 , P 2 ⊂ C 3 , which is defined by Subsequently, the price of a European call option that is driven by the Variance Gamma process is given by the formula where, for any c 1 ∈ P 1 and c 2 , ∈ P 2 , we define I 1 VG and I 2 VG by Proof. For the mean correcting martingale measure Q, the price of a European Call Option under the Variance Gamma process (25) is given by (see (30)) By definition, G 1 τ ∼ Gamma(Cτ, 1/M) and G 2 τ ∼ Gamma(Cτ, 1/G); hence, from (26), their probability density functions are given, respectively, by In order for the term (e [log]+x−y − 1) + to be different from zero, we must have y ≤ x + [log]. Additionally, notice that, by definition, both x and y are non-negative; hence, the values for which the integrand in (37) is not zero are given by the set {(x, y) ∈ R 2 : −[log] < x, y < x + [log]}, with the removal of the values in the set {(x, y) ∈ R 2 : −[log] < x < 0, y < x + [log]} for the cases where [log] > 0. Thus, C VG can be expressed as the sum where I 1 VG and I 2 VG are given by By the Mellin-Barnes representation of the exponential that is given by Equation (14), where c 11 , c 12 < 0, we can represent both terms e −Mx and e −Gy as the integrals c 11 +i∞ , respectively. Implementing these new representations on the integral of Equation (39) results in Applying integration by parts over the y variable to Equation (41) produces The Mellin-Barnes exponential representation for e [log]+x−y is given by the integral This substitution of terms in (42) will yield where c 1 was extended to the third dimension, i.e.,

This change of variables will result in the expression
Replacing (44), in the original integral (43) we obtain Similarly, for the integral By replacing (46) on the integral (45) and subsequntly inserting the resulting term on the original expression (39), we finally obtain Formula (34). This integral converges if all of the arguments of the Gamma functions in the numerator are positive. This happens Conversely, we can charter the same steps for the integral (40) of I 2 VG that we did for (39) of I 1 VG ; apply the Mellin-Barnes representation of the exponential term to both e −Mx and e −Gy . Subsequently, use integration by parts over the variable y, and again apply the Mellin-Barnes representation of e [log]+x−y , and finally apply the change of variables x := ([log] + x)s, in order to obtain Notice that the integral (47) Replacing (48) on the integral (47) and subsequently implanting the resulting expression, we obtain Formula (35). This integral formula converges if all of the arguments of the Gamma functions in the numerator are positive. This happens when (t),

Residue Summation Formula for a Call Option
Now, we will state and prove the main result of this paper: the triple representation formula for the European call option under the Variance Gamma model.

Theorem 3 (European Call Option
Price under the Variance Gamma Process). The price for the European Call-Option under the Variance Gamma process X VG (τ; C, G, M) is given by the formula where C 1 VG , C 2 VG and C 3 VG are defined by Proof. The Formula (34) for I 1 VG can be written as where c 1 is a three-dimensional point   c 11 c 12 c 13   ∈ P 1 and ω 1 VG is a complex differential threeform that is defined by By (16) . Thus, the admissible half-space Π ∆ 1 is the space that is located under this plane, Given the half-space Π ∆ 1 that is defined by (55), we will now construct an admissible polyhedron Π 1 = σ 1 := g −1 1 (G), as in the case (17), where G is the first octant, which will be uniquely defined by the linear function Under the linear function (56), the polyhedron Π 1 will be admissible, which is, and its faces σ 1 1 , σ 1 2 , and σ 1 3 and vertex σ 1 {1,2,3} will be Finally, we can group (where we used the notation t = (t, t x , t y )) into three families that are compatible with the polyhedron Π 1 , i.e., σ 1 thus verifying all of the conditions needed to apply Theorem 2. Note that this theorem is a version of the Jordan Lemma (Theorem 1) that allows us to represent the integral (53), which is an integral that appears in the formula for the call option price (33), as a sum of residues of the complex differential three-form ω 1 VG in (54). Before applying the residue summation, notice that the form ω 1 VG can be considered as having two sets of discontinuity points under the polyhedron Π 1 ; the first set is defined as S 1 := {t ∈ C 3 : t = −k, t x = −n, t y = −m, (k, n, m) ∈ N 3 }, which represents the singularity points that are given by the functions Γ(t), Γ(t x ), and Γ(t y ), and the second set is defined as S 2 := {t ∈ C 3 : t = −k, t y = −m, −1 − 2Cτ + t + t x + t y = −n, (k, n, m) ∈ N 3 }, which represents the discontinuity points that are given by the functions Γ(t), Γ(t x ) and Γ(−1 − 2Cτt + t x + t y ).
Given this delineated partition, we can now express equation (53) as where we define the terms C 1 VG and C 2 VG , respectively, by The computation of the residues for the first set S 1 present in the first series C 1 VG of (66) is straightforward Embedding the result (67) on the first equation of (66) will produce the sum Formula (50). On the other hand, for the computation of the residues of the set S 2 presented in the second series of (66), let us consider the variable change u := t, u y := t y , and u x := −1 − 2Cτ + t + t x + t y . If we apply these variables changes to the expression (54), the form ω 1 VG will be written as Subsequently, the residues of the second series Analogously to what we did for I 1 VG , the expression (35) of I 2 VG can be written as where c 2 is a three-dimensional point   c 21 c 22 c 23   ∈ P 2 and ω 2 VG is a complex differential threeform defined by Just like in the previous case, we use (16) . Thus, we conclude that the admissible half-space Π ∆ 2 is located under this plane, Similarly to what we did for Π ∆ 1 , given the expression (72) for Π ∆ 2 , we will now construct an admissible polyhedron Π 2 = σ 2 := g −1 2 (G), as in the case (17), where G is the first octant and g 2 is the linear function where I is the identity matrix. Under the linear function prescribed in (73), the polyhedron Π 2 will be admissible, and its faces σ 2 1 , σ 2 2 , and σ 2 3 and vertex σ 2 {1,2,3} will be Finally, we will group the divisors of ω 2 VG into three sets: that are compatible with polyhedron Π 2 ⊂ Π ∆ 2 , i.e., σ 2 1 ∩ D 1 = ∅, σ 2 2 ∩ D 2 = ∅, and σ 2 3 ∩ D 3 = ∅, thus satisfying all of the conditions needed to apply Theorem 2.
Unlike the ω 1 VG case, the form ω 2 VG under the polyhedron Π 2 has only one set of residues S 3 = {t ∈ C 3 : t = −k, t x = −n, t y = −m, (k, n, m) ∈ N 3 } that result from the functions Γ(t), Γ(t x ), and Γ(t y ). Therefore, Equation (35) can be expressed as The computation of the residues of set S 3 in the series of (82) is straightforward: Swapping the term in (83) with Equation (82) will yield the sum Formula (52). Because the derived expressions (50), (51), and (52) ascertain the equalities I 1 VG , this, in turn, proves expression (49).
The price formula for the European call option that is given by the expression (49) entails an easily derivable price for the European put option through the use of Put-Call parity:

The Greeks
Given the simple formula for the European call option deduced previously, one may inquire about the availability of an equally simple measure for risk exposure. The Greeks quantify the sensibility of the option price to changes in the model parameters. In this section, we will show the existence of a series of formulas for the ∆, Γ, ρ, and Θ functions, which will be obtained by a differentiation of (49) on the appropriate parameter.
Theorem 4 (The Greeks). The delta, gamma, rho, and theta functions for a European option under the Variance Gamma process X VG (τ; C, G, M) are given by:

Numerical Results
The formulas of the previous section can be heuristically observed to be sound. In order to do this, we start by comparing the results that were obtained by Formula (49) with both a Monte Carlo simulation for the European call option (under the Variance Gamma process) and the actual values that were observed on the market. We will also verify that it is well behaved, which is, its price for any initial stock value and its implied volatility smile are similar to the expected behavior observed in any stock. We also take these results to observe the speed of convergence of the new method. To conclude, we will study the behavior of the Greek functions (85), (89), (93), and (94), which were derived in the previous section, and then compare them to the ones in the Black-Scholes model.
We compare the option prices that were obtained from Formula (49) with the prices obtained from the Black-Scholes model and the Monte Carlo simulation of the Variance Gamma model. Our goal with this comparison is to test the viability and accuracy of Formula (49) against these well-known benchmark models and standard approaches to option pricing. We could also compare the option prices that were obtained from our formula with the ones obtained under the Variance Gamma model with other numerical techniques, such as the closed form formula in [6], the Fast Fourier Transform method of [7], the Fourier-cosine series proposed in [8], the finite differences method in [9], or the multinomial method of [10]. However, such a comparison of the numerical results from our formula with the results that were obtained from all of these numerical methods is out of the scope of this paper and would require much more numerical analysis and computational work.

Variance Gamma Formula Values and Behavior
For the aforementioned comparison, we will use the values of European call options with the S&P500 as the underlying asset, bought at the close of the market on 18 April 2002. According to [23], at the close of the market on 18 April 2002, we had a risk-free rate of return r = 1.9%, a dividend rate of q = 1.2%, and the stock price closed at S 0 = 1124.47, with a volatility parameter, for the Black-Scholes model, of σ = 0.1812 and risk neutral parameters C = 1.3574, G = 5.8704, and M = 14.2699 for the Variance Gamma model.
We will take advantage of 75 actual recorded values that are presented in Table 1, and make an error estimation for each model, by calculating their respective root mean square error, which is given by the formula RMSE = ∑ n i=1 (market price i − model price i ) 2 /n. Under this metric, Table 2 presents the deviations from the observed results.
Therefore, not only is the Formula (49) more expedient due to much lower computational time, but it also outperforms the Monte-Carlo method (and, consequently, the Black-Scholes by a wide margin). From Figure 1, one can also observe that the formulas for the European options under the Variance Gamma (49) and (84) present typical behavior. For example, a change in the initial asset price influences the price formula as one would expect, yet with some substantial differences from the Black-Scholes model.  Let us note that, in [23], the author estimates the model parameters and calculates the option prices for the same option data and the same Variance Gamma model while using the Fast Fourier transform method of [7], and it obtains a RMSE value of 3.56, which is a little smaller than the one we obtained. We have also calculated the RMSE by the same numerical method and we obtained a RMSE of 3.74. These small discrepancies are perhaps a result of small differences in the precise computation of the maturities of some options.
More importantly the Formula (49) also displays the volatility smile typically present in most assets, with values relatively close to the empirical observed on present asset, as can be seen from Figure 2.

Convergence of the Variance Gamma Formula
In order to study the numerical convergence and precision of the new formula, we must first realize that observing the value of each isolated term in the triple sum (49) is fallacious, since an unit increase of, for instance, parameter n will lead to the sum of an extra m × k terms, which may lead to an error of substantially higher magnitude than each individual term. Let us denote each term of the sum (49) by C VG (n, m, k). Given this notation, we can write C VG = ∑ ∞ n=0 ∑ ∞ m=0 ∑ ∞ k=0 C VG (n, m, k). To determine the values n max , m max , and k max , for which the sum C VG converges, for all values K and τ of Table 1, we will apply the Euclidean norm to the 189 resulting values from the three possible double sum series: the series computed by summing the terms of C VG for a fixed n and 0 ≤ m, k ≤ n, i.e., C n const VG (n) = ∑ n m=0 ∑ n k=0 C VG (n, m, k), the series computed by summing the terms of C VG for a fixed m and 0 ≤ n, k ≤ m, i.e., C m const VG (m) = ∑ m n=0 ∑ m k=0 C VG (n, m, k), and the series computed by summing the terms of C VG for a fixed k and 0 ≤ n, m ≤ k, i.e., C k const VG (k) = ∑ k n=0 ∑ k m=0 C VG (n, m, k). Table 3. Convergence of the three double series. From the values shown in Table 3, for any K and τ of Table 1, we can assure convergence with a two-decimal precision, when the sum of the three double sum series has a result lower than 0.005, for instance, n = 22, m = 27, and k = 7.

Double Series Sum Double Series Sum
Let us note that, in [20,22], the authors studied the convergence of their double series for the option price with a 10 −3 precision, while considering the price of a single European call option (with fixed maturity and strike) under the fractional diffusion process and the finite moment Lévy stable process. They have concluded that is enough to sum the terms in the series up to n = 6 and m = 6. The results are quite different, because, in their models the call option price is given by a double series (not a triple series), the convergence was estimated for a single option (we calculate the convergence for options of all maturities and strikes in our data set), and their underlying model is not the Variance Gamma model.

The Greek Formulas Behavior
We terminate this section by visualizing and contrasting the behavior of Greek functions under the two models. The Greek functions for the Variance Gamma model seem similar enough to the ones of the Black-Scholes model, yet they exhibit enough discrepancies to be worthy of note, primarily in the behavior of the Gamma and Theta functions, as can be seen in Figure 3. These discrepancies can have important consequences in establishing hedging strategies and the associated expected profit and loss. We refer to [22] for a detailed discussion of this subject in the case of the finite moment log-stable model. Furthermore, while some greek functions properties, such as ∆ C , Γ C , ρ C , −Θ C > 0, can be difficult to prove analytically from the Formulas (85), (89), (93), and (94), one can easily observe their validity, by generating the graphs of Figure 3 for a variety of different values of S, K, r, q and τ.

Conclusions
We have derived a triple Mellin-Barnes integral representation for the price of a European call option in a market model where a Variance Gamma process drives the price of the underlying asset. Subsequently, we applied the multidimensional residue calculus to the aforementioned integral and computed a triple sum series for the European call option price (49). Triple sum series for the delta, gamma, rho, and theta Greeks were also found by direct differentiation. When numerically tested, Formula (49) exhibited the behavior that was typically observed in the market for European options, and it outperformed the Monte-Carlo simulation method in both accuracy and computational time. The Greeks also displayed their usual behavior.
For practical applications, the simplicity of the aforementioned formulas (such as the lack of necessity of simulations for pricing European options, or of complex numerical schemes to compute the Greeks), coupled with their higher rates of precision and much lower computation time, makes them ideal for financial practitioners. For example, Formulas (85) and (89) can be used directly to generate a portfolio with optimal delta and gamma hedge strategies.
In terms of future research, the most obvious course of action would be to compare the performance of Formula (49), in terms of accuracy and computational time, to other semiclosed formulas such as the Bessel functions representation formula or the Fast Fourier Transform method for the price of a European call option under the Variance Gamma model (as discussed in [6,7]).
In terms of generalization of the theoretical results, the more pressing question will be the ability to use a similar reasoning to that presented in Section 3 to arrive at a sum series for the more general CGMY process and Generalized Tempered Stable process, even if this necessitates a higher dimensional Mellin-Barnes integral. One also might inquire further into the pricing of more complex financial instruments, like American or Barrier options, and especially instruments like Asian options, where the integral that is involved in their definition seems to make them a prime candidate for residue calculus. Proof. In the first part of the proof, we will show that the Mellin-Barnes integral in (15) converges in the set U that is defined by (23). Let us start by noting that, for any z ∈ C, we have the inequality for some constant c 1 , as well as the inequality The last inequality (A2) is verified if and only if the inequality − arg(z) (z) − (z) ≤ − π 2 | (z)| + | (z)| is valid. Inasmuch as we can write s j = re iθ , where r ∈ R + 0 and θ ∈ [−π, π], the inequality can be written as rθ sin θ − r cos θ ≤ − π 2 |r sin θ| + |r cos θ|, which holds for any r ∈ R + and θ ∈] − π, π[. If z does not intersect the set Z − 0 + i{0}, then, as |z| → ∞, we can apply the Stirling Formula (12) and, from the expressions (A1) and (A2), we know there exist constants c 1 and c 2 , such that On the other hand, if we constrict the real value, x, to a compact set K ⊂ R \ Z − 0 , x will be bounded and the gamma function will be continuous in the domain K + iR. Thus, we can denote its supremum as M = sup x∈K |x| < ∞ and infimum as m = inf x∈K |x|. Under this notation, we have e −M < e −|x| and e |x| < e M , additionally as |y| → ∞ we have |x + iy| = x 2 + y 2 ∼ (1 + |y|). Applying these properties to (A3) results in the inequalities for some constants k 1 and k 2 . Taking advantage of the inequalities (A4), we can bound the integrand of expression (15) by ∏ m on the corresponding half-spaces where, ξ 1 {1,2,3} , ξ 2 {1,2,3} , and ξ 3 {1,2,3} are defined by and l 1 , l 2 , l 3 , P {2,3} , P {1,3} , and P {1,2} are subsequently defined by l 1 = {x ∈ R 3 : Further, the computation of the integral ξ 1 over σ {2,3} will be analogous to ξ 2 over σ {1,3} and ξ 3 over σ {1,2} , as will the computation of ξ {2,3} over σ 1 be analogous to ξ {1,3} over σ 2 and ξ {1,2} over σ 3 . Thus, it will suffice to examine just three cases.
Let us consider the sequence of sets Define the surface S k = ∂U k and let ξ be one of the seven integrands of (A11)-(A13) defined on its corresponding σ of (A14) or (A15). It is a well-known property of Lebesgue integrals that there exists a constant c ∈ R, such that: where Z is a set with zero Lebesgue measure. We will prove that, as k → ∞, ξ → 0 at an exponential rate, which results in the integral in (A19), and, consequently, all (A41), (A45), and (A46) being zero. In order to achieve this, we will divide σ ∩ S k into two sets Depending on the σ that we are working with, this separation will yield different results. For instance, for σ 1 , we have As a final tool for our proof, consider the set U δ = {z ∈ C 3 : s j (z) + ν ≥ δ > 0, j = 1, . . ., m; ν ∈ N}, which removes a neighborhood in U around the singularities that are present in the numerator of the ratio of products of gamma functions in ω. Hence, we can write the left most term of (A19) as Our proof will consist of three steps: firstly, to estimate the value of ω | O k ∩U δ ; secondly, to estimate the value of ω | B k ∩U δ ; and thirdly, to extend the two previous results to the value of the integral to the set of neighborhoods U c δ (in the proof of Lemmas 1, 2, and 3).
Step 1: Recall that, for any z ∈ C 3 , we have s j (z) = a j , x + b j + i a j , y , also for any point in the three-dimensional space, x ∈ R 3 , there exist θ ∈ [0, 2π[ and φ ∈ [0, π[, such that x = x x = x (sin φ cos θ, sin φ sin θ, cos φ), hence the real part of s j (z) for z ∈ O k is given by (s j (z)) = R k a j ,x + b j . Having fixed θ and φ, and given that we are working in the space O k ∩ U δ , we will choose the radii R k , such that (s j (z)) will not intersect the singularities of the form ω, (the latter case will be dealt in the lemmas). Thus, the numerator of the form ω can be segregated into two terms, the ones where a j ,x > 0, which we will order as µ + 1 ≤ j ≤ m, and the ones where a j ,x < 0, which we will order as 1 ≤ j ≤ µ. Analogously, we can sort the denominator of ω, where c j ,x > 0, for χ + 1 ≤ j ≤ p, and c j ,x < 0, for 1 ≤ j ≤ χ. In this case, as respectively. Therefore, if we use the relation Γ(s j )Γ(1 − s j ) = π/ sin(πs j ) and then apply the Stirling Formula (12), we get where c 0 is an undetermined constant independent of either k or y values. Because all of the linear equations s j and q j are on the set O k , we have x = R k and y ≤ R k . Given this parametrization, the modulus of the functions s j , 1 − s j , q j , and 1 − q j are bounded below by A 1 R k and above by A 2 R k for some constants A 1 , A 2 > 0. Now, if we take note that R k a i ,x < 0 and −R k c j ≤ R k c j ,x < 0 for 1 ≤ i ≤ µ and χ ≤ j ≤ p, respectively, and R k a i ≥ R k a i ,x > 0 and R k c j ≥ R k c j ,x > 0 for 1 ≤ i ≤ µ and χ ≤ j ≤ p, respectively, we conclude that where c 1 , c 2 , c 3 , and c 4 are constants that we can define without any recourse to the angles θ and ψ, making the upper bound (A23) hold for every z ∈ O k ∩ U δ . Observe that, since we are working in O k ∩ U δ , there does not exist a sequence z i ∈ O k ∩ U δ , such that sin(π( a j , z + b j )) goes to zero; otherwise, inasmuch as sin and a j , z + b j are continuous, we would have z i ∈ U c δ , which is a contradiction, hence 0 < inf z∈O k ∩U δ sin(π( a j , z + b j )) and sup z∈O k ∩U δ 1/ ∏ µ j=1 sin(π( a j , z + b j )) < ∞. On the other hand, for the sin's in the denominator of ω, there exists a constant c, such that χ ∏ j=1 sin(πq j ) = χ ∏ j=1 e iπ (q j ) e −π (q j ) − e −iπ (q j ) e π (q j ) 2i ≤ χ ∏ j=1 e π| (q j )| ≤ e cR k .
A + B = ∆ x < 0, we have −B/A < 1. Consider the function f R (x) = x Ax R Bx , f R (x) has a maximum e − B A R −B/A . Therefore, even if we choose x = x M , the term x A x R B k will increase at a rate lower than e F(y) decreases. We conclude that, as k → ∞, the second integral of (A21) vanishes.
Step 3: We will now undertake the third integral of (A21) where the singularities of the form ω are present. Let us define P j,ν = {x + yi ∈ C 3 : a j , x = −b j − ν, a j , y = 0} and V j,ν ⊃ P j,ν as the open sets that contain the discontinuity that is given by the values in P j,ν . If we fix the angles θ x , φ x , θ y , φ y , for the aforementioned discontinuity, we can restate the previous equations as x a j ,x = −b j − ν and y a j ,ŷ = 0.
For the case where P j,ν ∩ O k = ∅, we have x = R k , and the real part of the singularities will be given by the one dimensional segments T k,j,ν = {x ∈ R 3 : x = R k , a j , x = −b j − ν}. Hence, we can chose the radii R k 's, such that, for each point x = R k 1x = R k (sin φ cos θ, sin φ sin θ, cos φ) belonging the two dimensional segment of T k 1 ,j,ν , when we increase the radius from R k 1 to R k 2 , the points with the same angles θ and φ, but now with the higher length R k 2 , i.e., R k 2x , will not intersect T k 2 ,j,µ , for any µ ∈ N and k 2 > k 1 , and will intersect the remaining T k 1 ,i,µ at most one time, for any µ ∈ N, j = i and k 2 > k 1 . We can construct the radii R k 's in order for this event to occur, because the set containing all of the two dimensional segments T k,j,ν is countable. In fact, we can define its (surjective) enumerating function as: where k ∈ N represents the radii R k , j ∈ {1, . . ., m} represents the Gamma function, and ν ∈ N the zero in the Gamma function. This entails that, as k → ∞, the upper bound (A23) deduced in step 1 will still hold for any angles θ x and φ x . The case where P j,ν ∩ B k = ∅ is slightly more complicated. Let us consider S j,ν = {i ∈ {1, 2, . . ., m} : ∀ z∈V j,ν a i , z / ∈ Z − 0 }, the multi-index set of the gammas that do not have a singularity under s i (z), and Z j,ν = {i ∈ {1, 2, . . ., m} : ∀ z∈V j,ν a i , z ∈ Z − 0 }, the multi-index set of the gammas that have a singularity under s i (z).
Because we are dealing with a case where x ∈ K \ Z − 0 , where K is compact, by the inequalities (A4), the estimate upper bound for the modulus of ωt −z is given by Consider α, as defined in (24). For each t in the domain U that is defined in (23), we have arg t < (π/2)α. If we apply the Cauchy-Schwartz inequality, we obtain: | y, arg t | ≤ π 2 α y for some α < α, Consider the parametrization that is given by x = p 0 1 , y = v 0 λ + v 1 θ + p 0 2 of the plane P j,ν ∩ B k , where (p 0 1 and p 0 2 are points and v 1 and v 2 are vectors). The neighborhood V jν around P jν ∩ B k can be defined as the intersection of the parallel planes {x + iy ∈ C 3 : x = p 1 , y = v 0 λ + v 1 θ + p 2 } with B k , where p 1 and p 2 are in the neighborhood of p 0 The first factor has a zero of order 2s, the second factor has a pole of order s, and the third factor decreases exponentially as R k increases. Therefore, the product of the first and second factors has a zero of order s. This, in conjunction with the previous deductions, completes the proof of lemma 1.

Lemma A2.
There exists a sequence of radii R k , such that R k → ∞, and for which