VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA232 10.15559/23-VMSTA232 Research Article Variance Gamma (nonlocal) equations https://orcid.org/0000-0001-5123-5916 ColantoniFaustofausto.colantoni@uniroma1.it∗ Department of Basic and Applied Sciences for Engineering, Sapienza University of Rome, Italy Corresponding author. 2023 259202310441342420220239620231282023 © 2023 The Author(s). Published by VTeX2023 Open access article under the CC BY license.

Some equations are provided for the Variance Gamma process using the definition other than that based on a time-changed Brownian motion. A new nonlocal equation is obtained involving generalized Weyl derivatives, which is true even in the drifted case. The connection to special functions is in focus, and a space equation for the process is studied. In conclusion, the convergence in distribution of a compound Poisson process to the Variance Gamma process is observed.

 Variance Gamma process fractional calculus Gamma subordinator nonlocal equations 60G51 34A08 Introduction

The Variance Gamma process is a famous Lévy process used in mathematical finance especially for option pricing (see [19, 22]), also known in physics as Laplace motion [15]. Another financial model involving this process has been considered quite recently in [2].

Further generalizations of the Variance Gamma process have been developed in [16, 10], while an application of such a process is given in [20].

The classic Variance Gamma can be obtained by considering a Brownian motion with a random time given by an independent Gamma subordinator [18, Eq. (7)]. The concept of subordination has been introduced by Bochner [5] and, as for the other subordinated processes, we can associate to the Variance Gamma process the Phillips operator [21] as generator for the governing equation. Recently, some new equations for the Variance Gamma process and the Gamma subordinator have been provided in [3], involving time-operators, differently from classic theory on space operators.

At the current stage, the nonlocal equation for the Variance Gamma process, defined as difference of two independent Gamma subordinators [18, Eq. (8)], is not considered and also a deeper analysis concerning this process and the equations for the modified Bessel functions is not taken into account. In order to close such a gap, we focus on these equations and the knowledge of these allows to understand the process from a new point of view, also for simulations. Then, we continue our dissertation by examining the similar version of the nonlocal equation for the Variance Gamma process with a drift. This is possible because the definition by difference of two independent Gamma subordinators holds in presence of drift also. At the end, we consider the compound Poisson process, related to the Gamma subordinator (the construction for any subordinator is developed in [25, Proposition 3.3]), and its convergence (in distribution) to the Variance Gamma process.

 Preliminaries

A subordinator is a nondecreasing Lévy process [4, Chapter III] and it is characterized by its Laplace transform [24, Theorem 5.1]. Let Φ : ( 0 , ∞ ) ↦ ( 0 , ∞ ) be a Bernstein function which is uniquely defined by the so-called Bernstein representation (Lévy–Kintchine representation) Φ ( λ ) = ∫ 0 ∞ ( 1 − e − λ z ) Π ( d z ) , λ > 0 , where Π on ( 0 , ∞ ) with ∫ 0 ∞ ( 1 ∧ z ) Π ( d z ) < ∞ is the associated Lévy measure. In general Bernstein functions can have a drift and a killing rate, but in this paper we are assuming them to be zero.

We also recall that Φ ( λ ) λ = ∫ 0 ∞ e − λ z Π ‾ ( z ) d z , λ > 0 , where Π ‾ ( z ) = Π ( ( z , ∞ ) ) is termed tail of the Lévy measure.

We focus only on the Laplace symbol Φ ( λ ) = a ln 1 + λ b = a ∫ 0 ∞ 1 − e − λ y e − b y y d y , λ > 0 , a > 0 , b > 0 . Thus, in this case, the Lévy measure is Π a , b ( d y ) = a e − b y y d y and the associated Gamma subordinator H = { H t , t ≥ 0 }, starting from zero, is such that E 0 [ e − λ H t ] = e − t Φ ( λ ) , λ > 0 , where P x denotes the probability measure for a process started from x at time t = 0 and E x the mean value with respect to P x . The interested reader can consult [4, Chapter III] for more details on subordinators. Since Π a , b ( 0 , ∞ ) = ∞, from [13, Theorem 21.3] we have that H has infinite activity, i.e. strictly increasing sample paths with jumps; indeed the symbol Φ does not admit any drift. We use the notation P x ( H t ∈ d x ) = h ( t , x ) d x . In addition, it is well known that, ∀ t > 0, h ( t , x ) = b a t Γ ( a t ) x a t − 1 e − b x , x > 0 , 0 , x ≤ 0 , trivially verifies ∫ 0 ∞ e − λ x h ( t , x ) d x = b a t ( λ + b ) a t = e − t a ln 1 + λ b , λ > 0 , t > 0 , which coincides with formula (3). We observe that the continuity of the function h ( t , · ), when x → 0, depends on the time variable t. Indeed, the Gamma subordinator has the time dependent property (see [13, Chapter 23]). From this, the Variance Gamma process inherits also continuity problems for its probability density function.

Let B : = { B t , t ≥ 0 } be the one-dimensional Brownian motion starting from zero, independent of H, and g ( t , x ) = e − x 2 / 4 t / 4 π t be its probability density function. The Variance Gamma process X : = { X t , t ≥ 0 } can be defined as X = B H : = B ∘ H, i.e. it is a Brownian motion time-changed with a random clock given by an independent Gamma subordinator, and we know from [13, Theorem 30.1] that X is still a Lévy process. Its probability density function is p ( t , x ) = ∫ 0 ∞ g ( s , x ) h ( t , s ) d s , and the Lévy symbol is Φ ( ξ 2 ), as we see from E 0 [ e i ξ X t ] = p ˆ ( t , x ) = ∫ − ∞ ∞ e i ξ x p ( t , x ) d x = ∫ − ∞ ∞ e i ξ x ∫ 0 ∞ g ( s , x ) h ( t , s ) d s d x = ∫ 0 ∞ e − s ξ 2 h ( t , s ) d s = e − t Φ ( ξ 2 ) = e − a t ln ( 1 + ξ 2 b ) . For the Variance Gamma process we know an explicit representation for p ( t , x ). From [11, formula 3.478] we have ∫ 0 ∞ x ν − 1 exp { − β x q − α x − q } d x = 2 q α β ν 2 q K ν q 2 α β , q , α , β , ν > 0 , where K ν is the modified Bessel function, so we get that ∫ 0 ∞ g ( s , x ) h ( t , s ) d s = ∫ 0 ∞ e − x 2 4 s 4 π s b a t Γ ( a t ) s a t − 1 e − b s d s = = b a t 4 π 1 Γ ( a t ) ∫ 0 ∞ s a t − 1 s exp − x 2 4 s − 1 − b s d s . Then we use (8) by choosing ν = a t − 1 2 , q = 1, α = z 2 4 , β = b and we obtain p ( t , x ) = ∫ 0 ∞ g ( s , x ) h ( t , s ) d s = b a t 4 π 1 Γ ( a t ) 2 x 2 4 b 1 2 ( a t − 1 2 ) K a t − 1 2 2 x 2 4 b = b a t π 1 Γ ( a t ) 1 2 b ( a t − 1 2 ) | x | a t − 1 2 K a t − 1 2 ( | x | b ) .

Since we are dealing with a subordinate stochastic process (see [6, 14]), we can use the Phillips representation (or the Bochner subordination) to generate a new operator through subordination. It is well known that the Phillips operator − Φ ( − Δ ) [21] is the generator of X, when it is restricted to D o m ( Δ ) ([13, Theorem 32.1] or [12, Theorem 4.3.5]), and it can be written as − Φ ( − Δ ) u ( x ) = a ∫ 0 ∞ ( T y u ( x ) − u ( x ) ) e − b y y d y , where Φ is defined in (2) and T is the semigroup associated to the Brownian motion on R, such that the characteristic symbol is T ˆ t = e − t ξ 2 . This leads us to ∫ − ∞ ∞ e i ξ x − Φ ( − Δ ) u ( x ) d x = a ∫ 0 ∞ ( e − y ξ 2 − 1 ) e − b y y d y u ˆ ( ξ ) = − a ln 1 + ξ 2 b u ˆ ( ξ ) , where we indicate by u ˆ the Fourier transform of u, for u in the Schwartz space S ( R ). If we combine (7) and (12), we obtain the well known result ∂ ∂ t p ( t , x ) = − Φ ( − Δ ) p ( t , x ) , t > 0 , x ∈ R , with p ( t , 0 ) = δ ( x ) in the sense that the Fourier transform of the fundamental solution of (13) is satisfied by p.

If we change perspective and consider the operator in time, an interesting equation, presented in [3, Remark 3.3], for the Variance Gamma process is ∂ 2 ∂ x 2 p ( t , x ) = b p ( t , x ) − p t − 1 a , x , x ∈ R , with p ( t , 0 ) = δ ( x ) and a t > 1.

The condition a t > 1 is not surprising, since h ( t , · ) is continuous when a t > 1.

 Main results Nonlocal equations

In the preceding section we have seen how in (13), from the definition of time-changed process, we can associate to the Variance Gamma process a nonlocal operator. However, it is not the only possible definition. From (7), we see that e − a t ln 1 + ξ 2 b = 1 + ξ 2 b − a t = 1 − i ξ b − a t 1 + i ξ b − a t . Then, in the Wiener–Hopf factorization [17, Remark 3.3], the Bernstein function associated to the Gamma subordinator (2) is used. Therefore, we have another definition for our process X = G − L where G and L are two independent Gamma subordinators both with parameters a and b . From a financial point of view this representation has the meaning of the difference between independent ‘gains’ and ‘losses’ (see [18]). Formula (14) means that X has trajectories of locally bounded variation, as it is the difference of two increasing functions, and this leads us to investigate the connection between the Variance Gamma process and the nonlocal operators of the Gamma subordinator.

Let us introduce the following generalized Weyl derivatives for x ∈ R: D a , b + u ( x ) : = ∂ ∂ x ∫ − ∞ x u ( s ) Π ‾ a , b ( x − s ) d s , D a , b − u ( x ) : = − ∂ ∂ x ∫ x ∞ u ( s ) Π ‾ a , b ( s − x ) d s , respectively defined for function u such that u ( s ) Π ‾ a , b ( x − s ) ∈ L 1 ( − ∞ , x ) and u ( s ) Π ‾ a , b ( s − x ) ∈ L 1 ( x , ∞ ) , ∀ x ∈ R .

We could also refer to (15) and (16) as generalized Riemann–Liouville derivatives on the whole real axis (see [23, Formula (5.6)]).

From [7, Eq. (37)] we know that Π ‾ a , b ( x ) = a E 1 ( b x ), where E 1 is the exponential integral E 1 ( x ) : = ∫ x ∞ e − z z d z , x > 0 . This allows us to write (15) and (16) in a compact form: D a , b + u ( x ) : = a ∂ ∂ x ∫ − ∞ x u ( s ) E 1 ( b ( x − s ) ) d s , D a , b − u ( x ) : = − a ∂ ∂ x ∫ x ∞ u ( s ) E 1 ( b ( s − x ) ) d s .

From a simple change of variables, by assuming that u is at least absolutely continuous, our definitions (15) and (16) concur with [25, Definition 2.8] and we achieve D a , b + u ( x ) : = ∫ 0 ∞ ∂ ∂ x u ( x − s ) Π ‾ a , b ( s ) d s , D a , b − u ( x ) : = − ∫ 0 ∞ ∂ ∂ x u ( x + s ) Π ‾ a , b ( s ) d s . The operators D a , b + and D a , b − are important due to the fact that we can easily compute the characteristic function essential to the study of Lévy processes. Indeed, by recalling (1), for u ∈ S ( R ) we have ∫ − ∞ ∞ e i x ξ D a , b + u ( x ) d x = ( − i ξ ) ∫ − ∞ ∞ e i x ξ ∫ 0 ∞ u ( x − s ) Π ‾ a , b ( s ) d s d x = ( − i ξ ) u ˆ ( ξ ) ∫ 0 ∞ e i s ξ Π ‾ a , b ( s ) d s = a ln 1 − i ξ b u ˆ ( ξ ) . Similarly we obtain ∫ − ∞ ∞ e i x ξ D a , b − u ( x ) d x = a ln 1 + i ξ b u ˆ ( ξ ) . These results coincide with [25, Lemma 2.9]. We now analyze the connection between these operators and the Variance Gamma process by exploiting (14). In particular, we will show that the Fourier transform of its probability density function solves the following equation.

Let p ( t , x ) be the probability density function of the Variance Gamma process X. Then we have ∂ ∂ t p ( t , x ) = − D a , b + p ( t , x ) + D a , b − p ( t , x ) , t > 0 , x ∈ R , with p ( 0 , x ) = δ ( x ).

The initial value can be easily checked, since the Fourier transform of the δ distribution is 1. The characteristic function of the left-hand side, with respect to x, is given by ∫ − ∞ ∞ e i x ξ ∂ ∂ t p ( t , x ) d x = ∂ ∂ t e − a t ln 1 + ξ 2 b = − a ln 1 + ξ 2 b e − a t ln 1 + ξ 2 b . For the right-hand side, by using (18) and (19), we have − ∫ − ∞ ∞ e i x ξ D a , b + p ( t , x ) + D a , b − p ( t , x ) d x = − a ln 1 − i ξ b e − a t ln 1 + ξ 2 b + a ln 1 + i ξ b e − a t ln 1 + ξ 2 b = − a ln 1 + ξ 2 b e − a t ln 1 + ξ 2 b , so that our claim holds true.  □

The fact that, in the last theorem, there is a sum of nonlocal derivatives is not new in the theory of probability and nonlocal operators. For example, for the Brownian motion time-changed with an independent stable subordinator, we know that the one-dimensional Riesz derivative can be written as sum of Marchaud derivatives (see [9, page 12]).

From the financial definition of X, the effect of the difference between independent gains and losses should be studied through the two nonlocal operators (15) and (16). Indeed, the last result transforms this definition in a nonlocal equation for the process.

In Theorem 1 we took advantage of the possible definition of the process X as difference of two independent Gamma subordinators. We observe that this information on the process brings us to the next result for the operators.

Let u ∈ S ( R ). Then the following equivalence is satisfied − Φ ( − Δ ) u ( x ) = − D a , b + u ( x ) + D a , b − u ( x ) , x ∈ R .

We easily show this result by using the characteristic functions. On the left-hand side we have, from (12), ∫ − ∞ ∞ e i ξ x − Φ ( − Δ ) u ( x ) d x = − a ln 1 + ξ 2 b u ˆ ( ξ ) . On the right-hand side, from (18) and (19), we have − ∫ − ∞ ∞ e i x ξ D a , b + u ( x ) + D a , b − u ( x ) d x = = − a ln 1 − i ξ b u ˆ ( ξ ) + a ln 1 + i ξ b u ˆ ( ξ ) = − a ln 1 + ξ 2 b u ˆ ( ξ ) , which concludes the proof.  □

Adding the drift. We now focus on the Variance Gamma process with drift and we prove that the nonlocal equation is still true. Let B θ : = { B t θ , t ≥ 0 } be the drifted Brownian motion on R starting from zero, independent of H, and g θ ( t , x ) = e − ( x − θ t ) 2 / 4 t / 4 π t , for the drift θ ∈ R, be its probability density function. The drifted Variance Gamma process X θ : = { X t θ , t ≥ 0 } is X θ = B H θ : = B θ ∘ H. Its characteristic function turns out to be E 0 [ e i ξ X t θ ] = 1 − i ξ θ b + ξ 2 b − a t . We see that, again, the process can be written as difference of two independent Gamma subordinators as suggested by [19, Eq. (8)]. Indeed, the following holds 1 − i ξ θ b + ξ 2 b − a t = 1 − i ξ θ 2 4 + b − θ 2 − a t 1 + i ξ θ 2 4 + b + θ 2 − a t . Thus we have X θ = G θ − L θ , where G θ and L θ are two independent Gamma subordinators, the first one with parameters a and θ 2 4 + b − θ 2 and the second one with parameters a and θ 2 4 + b + θ 2 . As well as for the Variance Gamma process, X θ admits paths of locally bounded variation and we show the next result.

Let p θ ( t , x ) be the probability density function of the drifted Variance Gamma process X θ . Then we have ∂ ∂ t p θ ( t , x ) = − D a , θ 2 4 + b − θ 2 + p θ ( t , x ) + D a , θ 2 4 + b + θ 2 − p θ ( t , x ) , t > 0 , x ∈ R , with p θ ( 0 , x ) = δ ( x ).

We observe that p θ ( t , x ) can be written as (6), where g θ replaces g. By using that 1 − i ξ θ b + ξ 2 b = 1 − i ξ θ 2 4 + b − θ 2 1 + i ξ θ 2 4 + b + θ 2 , the proof is analogous to that of Theorem 1 and the Fourier transform of p θ solves the equation.  □

Trivially, if θ = 0, the last differential equation coincides with that of Theorem 1 as expected.

 Variance Gamma and special functions

In this section we examine how special functions bring us to a new equation for the Variance Gamma process. From (10), we have seen that the modified Bessel function K ν ( x ) appears in the density and this special function solves (see [1, formula 9.6.1]) x 2 ∂ 2 ∂ x 2 u ( x ) + x ∂ ∂ x u ( x ) − ( x 2 + ν 2 ) u ( x ) = 0 , where, in our case, ν is a time function. Another interesting fact on the function K ν ( x ) is that it can be written in terms of the Kummer (confluent hypergeometric) function usually denoted by U, so that it is connected to the Kummer equation. This moves our thinking about differential space equations for the Variance Gamma process.

The following differential equation is satisfied by the density of the Variance Gamma process X, for t > 0 fixed: x ∂ 2 ∂ x 2 p ( t , x ) − ( 2 a t − 2 ) ∂ ∂ x p ( t , x ) − b x p ( t , x ) = 0 , x ∈ R , with p ( 0 , x ) = δ ( x ).

The initial value can be easily checked, since the Fourier transform of the δ is 1. The characteristic function, with respect to x, of (20) is ∫ − ∞ ∞ e i x ξ x ∂ 2 ∂ x 2 p ( t , x ) − ( 2 a t − 2 ) ∂ ∂ x p ( t , x ) − b x p ( t , x ) d x = = − i ∂ ∂ ξ ∫ − ∞ ∞ e i x ξ ∂ 2 ∂ x 2 p ( t , x ) d x − ( 2 a t − 2 ) ( − i ξ ) ∫ − ∞ ∞ e i x ξ p ( t , x ) d x − b ( − i ) ∂ ∂ ξ ∫ − ∞ ∞ e i x ξ p ( t , x ) d x = − i ∂ ∂ ξ ( − ξ 2 ) e − a t ln 1 + ξ 2 b + ( 2 a t − 2 ) ( i ξ ) e − a t ln 1 + ξ 2 b − i b 2 a t ξ b 1 + ξ 2 b − a t − 1 = 2 i ξ 1 + ξ 2 b − a t − 2 i ξ 2 a t ξ b 1 + ξ 2 b − a t − 1 + ( 2 a t − 2 ) ( i ξ ) 1 + ξ 2 b − a t − 2 i a t ξ 1 + ξ 2 b − a t − 1 = − 2 i ξ 2 a t ξ b 1 + ξ 2 b − a t − 1 + 2 i a t ξ 1 + ξ 2 b − a t − 2 i a t ξ 1 + ξ 2 b − a t − 1 = − 2 i ξ 2 a t ξ b 1 + ξ 2 b − a t − 1 + 2 i a t ξ 1 + ξ 2 b 1 + ξ 2 b − a t − 1 − 2 i a t ξ 1 + ξ 2 b − a t − 1 = − 2 i ξ 2 a t ξ b 1 + ξ 2 b − a t − 1 + 2 i ξ 2 a t ξ b 1 + ξ 2 b − a t − 1 + 2 i a t ξ 1 + ξ 2 b − a t − 1 − 2 i a t ξ 1 + ξ 2 b − a t − 1 = 0 , so that the Fourier transform of p solves (20), as required.  □

 Compound Poisson process convergence

The convergence of compound Poisson processes to subordinators has been extensively studied. In [8, Theorem 1] M. D’Ovidio shows the connection to a difference of α-stable subordinators and in [25, Proposition 3.3] we have the construction for any subordinator. In this section, we exploit these mentioned results to obtain the convergence of a compound Poisson process to the Variance Gamma process.

We consider the i.i.d. random variables Y j , Y j ∼ Y (distributed as Y), with probability density function ν Y ( y ) = e − b y y 1 E 1 ( b γ ) 1 y ≥ γ , γ > 0 , where E 1 is defined in (17). Let ϵ j ∼ ϵ be i.i.d. (centered) Rademacher random variables, with the distribution P ( ϵ = + 1 ) = 1 2 , P ( ϵ = − 1 ) = 1 2 . Now, we define Y ∗ = ϵ Y with probability density function ν Y ∗ = 1 2 ν Y ( − y ) + 1 2 ν Y ( y ) . We are ready to provide the following convergence result.

Let N ( t ), t ≥ 0, be a homogeneous Poisson process with parameter 1, independent of the i.i.d random variables Y j ∗ ∼ Y ∗ , with law (21). We have that ∑ j = 0 N ( t a E 1 ( b γ ) ) Y j ∗ → l a w X t 2 as γ → 0 .

Since we are dealing with a compound Poisson process, we know that E 0 exp i ξ ∑ j = 0 N ( t a E 1 ( b γ ) ) Y j ∗ = exp t a E 1 ( b γ ) ( E 0 exp ( i ξ Y ∗ ) − 1 ) = exp t a E 1 ( b γ ) ( E 0 exp ( i ξ ϵ Y ) − 1 ) = exp t a E 1 ( b γ ) 1 2 E 0 exp ( i ξ Y ) + 1 2 E 0 exp ( − i ξ Y ) − 1 2 + 1 2 = exp t a E 1 ( b γ ) 1 2 E 0 [ exp ( i ξ Y ) − 1 ] + 1 2 E 0 [ exp ( − i ξ Y ) − 1 ] = exp t a 2 ∫ γ ∞ ( e i ξ y − 1 ) e − b y y d y + a 2 ∫ γ ∞ ( e − i ξ y − 1 ) e − b y y d y . If γ → 0, we get that exp t a 2 ∫ γ ∞ ( e i ξ y − 1 ) e − b y y d y + a 2 ∫ γ ∞ ( e − i ξ y − 1 ) e − b y y d y → exp − t 2 a ln 1 − i ξ b + a ln 1 + i ξ b , hence the claim holds by using (14) and the classic Lévy continuity theorem.  □

This result may be important for the simulations: numerical analysis for compound Poisson process is simpler than that for Bessel functions of (10).

 Acknowledgement

The author wishes to thank Mirko D’Ovidio and the anonymous reviewers for their helpful advice. The author also thanks Sapienza University of Rome and the group INdAM-GNAMPA for their support.

 References Abramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, vol. 55. US Government Printing Office (1964). MR0167642 Ascione, G., Mehrdoust, F., Orlando, G., Samimi, O.: Foreign exchange options on Heston-CIR model under Lévy process framework. Appl. Math. Comput. 446 (2023). MR4550246 Beghin, L.: Geometric stable processes and related fractional differential equations. Electron. Commun. Probab. 19, 1–14 (2014). MR3174831. https://doi.org/10.1214/ECP.v19-2771 Bertoin, J.: Lévy Processes, vol. 121. Cambridge University Press, Cambridge (1996). MR1406564 Bochner, S.: Diffusion equation and stochastic processes. Proc. Natl. Acad. Sci. 35(7), 368–370 (1949). MR30151. https://doi.org/10.1073/pnas.35.7.368 Bogdan, K., Byczkowski, T., Kulczycki, T., Ryznar, M., Song, R., Vondracek, Z., Song, R., Vondraček, Z.: Potential theory of subordinate brownian motion. In: Potential Analysis of Stable Processes and its Extensions, 87–176 (2009) MR2569321. https://doi.org/10.1007/978-3-642-02141-1 Colantoni, F., D’Ovidio, M.: On the inverse gamma subordinator. Stoch. Anal. Appl. 41(5), 999–1024 (2023). MR4627143. https://doi.org/10.1080/07362994.2022.2108450 D’Ovidio, M.: Continuous random walks and fractional powers of operators. J. Math. Anal. Appl. 411(1), 362–371 (2014). MR3118491. https://doi.org/10.1016/j.jmaa.2013.09.048. Ferrari, F.: Weyl and Marchaud derivatives: A forgotten history. Mathematics 6(1), 6 (2018). https://doi.org/10.3390/math6010006 Gajda, J., Wyłomańska, A., Kumar, A.: Generalized fractional Laplace motion. Stat. Probab. Lett. 124, 101–109 (2017). MR3608219 Gradshteyn, I.S., Ryzhik, I.M.: Table of Integrals, Series, and Products, 7th edn., p. 1171. Elsevier/Academic Press, Amsterdam (2007). MR2360010 Jacob, N.: Pseudo Differential Operators and Markov Processes. Vol. I, p. 493. Imperial College Press, London (2001). Fourier analysis and semigroups. MR1873235 Ken-Iti, S.: Lévy Processes and Infinitely Divisible Distributions. Cambridge University Press (1999). MR3185174 Kim, P., Song, R., Vondraček, Z.: Potential theory of subordinate Brownian motions revisited. In: Stochastic Analysis and Applications to Finance. Interdiscip. Math. Sci., vol. 13, pp. 243–290. World Sci. Publ., Hackensack, NJ (2012). MR2986850 Kotz, S., Kozubowski, T.J., Podgórski, K.: The Laplace Distribution and Generalizations. Birkhäuser Boston, Inc., Boston, MA (2001). MR1935481. https://doi.org/10.1007/978-1-4612-0173-1 Kozubowski, T.J., Meerschaert, M.M., Podgórski, K.: Fractional Laplace motion. Adv. Appl. Probab. 38(2), 451–464 (2006). MR2264952 Kwaśnicki, M.: Spectral analysis of subordinate Brownian motions on the half-line. Stud. Math. 206(3), 211–271 (2011). MR2860308 Madan, D.B., Seneta, E.: The variance gamma (vg) model for share market returns. J. Bus. 63(4), 511–524 (1990). https://doi.org/10.1086/296519 Madan, D.B., Carr, P.P., Chang, E.C.: The variance gamma process and option pricing. Rev. Finance 2(1), 79–105 (1998). https://doi.org/10.1023/A:1009703431535 Meerschaert, M.M., Kozubowski, T.J., Molz, F.J., Lu, S.: Fractional Laplace model for hydraulic conductivity. Geophysical Research Letters 31(8) (2004) Phillips, R.S.: On the generation of semigroups of linear operators. Pac. J. Math. 2(3), 343–369 (1952). MR50797 Reed, W.J.: Brownian-Laplace motion and its use in financial modelling. Commun. Stat., Theory Methods 36(1-4), 473–484 (2007). MR2405311 Samko, S.G., Kilbas, A.A., Marichev, O.I.: Fractional Integrals and Derivatives. Gordon and Breach Science Publishers, Yverdon (1993). Theory and applications. MR1347689 Schilling, R.L., Song, R., Vondraček, Z.: Bernstein Functions, 2nd edn. De Gruyter Studies in Mathematics, vol. 37. Walter de Gruyter & Co., Berlin (2012). MR2978140 Toaldo, B.: Convolution-type derivatives, hitting-times of subordinators and time-changed C 0 -semigroups. Potential Anal. 42(1), 115–140 (2015). MR3297989