Computation of Greeks using Malliavin’s calculus in jump type market models

We use the Malliavin calculus for Poisson processes in order to compute sensitivities for European and Asian options with underlying following a jump type diﬀusion. The main point is to settle an integration by parts formula (similar to the one in the Malliavin calculus) for a general multidimensional random variable which has an absolutely continuous law with diﬀerentiable density. We give an explicit expression of the diﬀerential operators involved in this formula and this permits to simulate them and consequently to run a Monte Carlo algorithm.


Introduction
In the last years, following the pioneering papers [9], [8] a lot of work concerning the numerical applications of the stochastic variational calculus (Malliavin calculus) has been done. This mainly concerns applications in mathematical finance: computations of conditional expectations (which appear in the american option pricing, for example) and of sensitivities (the so called Greeks). The models at hand are usually log-normal type diffusions and then one may use the standard Malliavin calculus. But nowadays people are more and more interested in jump type diffusions (see [5]) and then one has to use the stochastic variational calculus corresponding to Poisson point processes. Such a calculus has already been developped in [2] concerning the noise coming from the amplitudes of the jumps and in [4] concerning the jump times. More recently [3] gives a unified approach using the language of the Dirichlet forms. Another point of view, based on chaos decomposition may be found in [7], [1], [12], [14] and [11].
The aim of our paper is to give a concrete application of the Malliavin calculus approach for sensitivity computations (Greeks) for jump diffusion models. We give two examples: in the first one, we use the Malliavin calculus with respect to the jump amplitudes. In the second one, we add a Brownian part and differentiate with respect to both the jump amplitudes and to the Wiener increments.
The Malliavin integration by parts formula used in this paper is an elementary one. Notice that any numerical scheme used in a Monte Carlo algorithm appears as a function of a finite number of random variables H 1 , . . . , H n . It turns out that, if the law of these random variables is absolutely continuous with respect to the Lebesgue measure and has a smooth density, then the strategy developed in the Malliavin calculus may be implemented for H 1 , . . . , H n and an integration by parts formula is derived. This is not specific for the Brownian motion or for a Poisson point process but represents an elementary abstract calculus.
The paper is organized as follows. In Section 1 we present the problem: we give the stochastic equation verified by the underling asset which is a one dimensional jump diffusion (one may consider multi-dimensional diffusions as well but we stay in the one dimensional setting just to avoid technical difficulties). Then we precise how an integration by parts formula may be used in the Monte Carlo algorithm to compute sensitivities. In Section 2 we give the integration by parts formula. In Section 3 we use the Malliavin's calculus with respect to the amplitudes of the jumps in order to compute sensitivities. In section 4, we consider a Merton's model in which stochastic integrals driven by a Brownian motion and by a compound Poisson process appear. We have the choice of using the Malliavin's calculus with respect to the Brownian motion, to the amplitudes of the jumps or to both of them. We perform the three algorithms and we compare the empirical variance of each of them. It turns out that the most performant algorithm (with the smallest variance) is that one based on both the Brownian motion and the amplitudes of the jumps. So the numerical experiments (there is no theorical result) indicate that one has to use the integration by parts formula with respect to all the noise which is available in the model. The numerical results are contained in section 5. We compare the Malliavin method to the finite difference method. For European and Asian call options, the results are rather similar but for digital options, the results given by the Malliavin method are much better than these given by the finite difference method. This can be explained by the larger irregularity of the payoff function.

The problem
We compute the sensitivities of an option with payoff φ, where the asset (S t ) t≥0 follows a one-dimensional jump diffusion process driven by a compound Poisson process. Let us precise our problem. Let N be a compound Poisson process. We denote by (T i ) i∈N the jump times of t → N t , and by (J t ) t≥0 the corresponding standard Poisson process with parameter λ > 0, that is J t counts the number of jumps up to t and for all n ≥ 1, T n − T n−1 has exponential distribution with mean λ. For all i ∈ N * , we define ∆ i = N T i − N T i−1 , so ∆ i represents the jump amplitude at time T i . The random variables (∆ i ) i∈N * are independent and identically distributed. We assume that the law of ∆ i is absolutely continuous with respect to the Lebesgue measure and we denote by ρ its density, that is ∆ i ∼ ρ(a) da, for all i ∈ N * .
We deal with two models. In the first one, the asset S t is a pure jump diffusion solution of the SDE In the second one, we add a Brownian part : Our aim is to compute , that is the Delta of an European (respectively Asian) option of payoff φ. In several papers (such as [8], [9], [10]), a Malliavin calculus approach is used. We will follow a similar strategy in our frame, using Malliavin calculus for jump type diffusions. In the case of European option for example, we write Assume now that J T = n = 0. Then, using an integration by parts formula on {J T = n}, we write : where H n is a weight involving "Malliavin derivatives" of S T and of ∂S T ∂β . These differential operators are similar to those in [2], but the frame here is much more simple, since there are no accumulation of small jumps. So we will derive this integration by parts formula using elementary arguments. We obtain In order to employ this formula in a Monte-Carlo algorithm, we proceed as follows: we simulate the jump times (T i n ) n∈N , i = 1, . . . , M and the jump amplitudes (∆ i n ) n∈N , i = 1, . . . , M , and we compute the corresponding In the case of Merton's model we proceed as follows: we first simulate the jump times T i , i ∈ N.
Once they are fixed we construct an Euler scheme in which the times T i , i ∈ N are included in the discretization grid. Then we simulate the amplitudes of the jumps and the Brownian increments and we thus perform the Monte-Carlo algorithm.

Malliavin calculus for simple functionals
On a probability space (Ω, F, P ) we consider a sequence of independent random variables (V n , n ∈ N * ). We suppose that for all n ≥ 1, V n has moments of any order.
Hypothesis 3.1 For every n ∈ N * , the law of V n is absolutely continuous with respect to the Lebesgue measure and has the density p n which is continuously differentiable on R and such that for all k ∈ N, lim y→±∞ |y| k p n (y) = 0. We also assume that ∂ y ln(p n (y)) = ∂ y p n (y) p n (y) has at most polynomial growth.
We introduce some notations. For k ≥ 1, we denote by C k ↑ (R n ) the space of the functions f : R n → R which are k times differentiable and such that f and its derivatives up to order k have at most polynomial growth. For a multi-index α = (α 1 , . . . , α k ) we denote ∂ k α f = In the case k = 0, C 0 ↑ (R n ) denotes the set of the functions f : R n → R which are continuous with respect to each argument and have at most polynomial growth.
We define now the simple functionals and the simple processes. A random variable F is called a simple functional if there exists some n ∈ N * and some measurable function f : R n → R such that We denote by S (n,k) the space of the simple functionals such that f ∈ C k ↑ (R n ). A simple process of length n is a sequence of random variables We denote by P (n,k) the space of the simple processes of length n such that u i ∈ C k ↑ (R n ), i = 1, . . . , n. Note that if U ∈ P (n,k) then U i ∈ S (n,k) , i ∈ N * .
We define now the differential operators which appear in the Malliavin's calculus.
The Malliavin derivative D : The Skorohod integral δ : P (n,1) → S (n,0) : Remark 3.1 For every p ∈ N and for every F ∈ S (n,1) and U ∈ P (n,1) , one has This comes from the following hypothesis: E |V i | p < ∞, for all i and f , u i , and ∂ ln(p i ), i = 1, . . . , n, have at most polynomial growth.
The operator δ is the adjoint of D (this is the reason of being of θ i ): Proposition 3.1 Let F ∈ S (n,1) and U ∈ P (n,1) then where < ., . > is the scalar product in R n .
Proof. One writes the last equality being obtained by integration by parts, where we have used that for all k ∈ N, Coming back to expectations we obtain (3).
We introduce now the Ornstein Uhlenbeck operator L = δ(D) : S (n,2) → S (n,0) : Suppose that F, G ∈ S (n,2) then, as an immediate consequence of the duality relation (3) one obtains Moreover, the standard differential calculus rules give the following chain rule. 1) and (ii) Let F, G ∈ S (n,2) then Finally, given F = (F 1 , . . . , F d ), F i ∈ S (n,1) , we introduce the Malliavin covariance matrix We are now able to state the integration by parts formula. 1) . We assume that the matrix σ F is invertible and denote γ F := σ −1 F . We also suppose that Then for every smooth function φ : with H i (F, G) ∈ L 1 and Proof. Using the chain rule we obtain for each j = 1, . . . , d Then the duality relation gives By hypothesis (8), the above expectations are finite, so the proof is complete.

Integration by parts with respect to the jump amplitudes
In this section we use the integration by parts formula (9) with respect to ∆ i ∼ ρ(a) da, i ∈ N * , which are the amplitudes of the jumps of (N t ) t≥0 . Let (S t ) t≥0 be solution of the SDE We work under the following hypothesis.

The deterministic equation
We introduce the following deterministic equation. We fix some 0 < u 1 < . . . < u n < T , and we denote u = (u 1 , . . . , u n ). We also fix a = (a 1 , . . . , a n ) ∈ R n . To u and a we associate the equation We use this deterministic equation in order to express S t as a simple functional: So in order to compute the derivatives of S t , we have to compute those of s t . Under hypothesis 4.2, we may perform these computations by differentiating with respect to a j , j = 1, . . . , n in (12). We obtain and for i > j, we obtain ∂ 2 a j ,a i s t by symmetry.
(ii) Using (13), ∂ an s T = ∂ a c(t n , a n , s t − n ) + T tn ∂ x b(r, s r ) ∂ an s r dr, so

Integration by parts formula
We come back to the problem developed in section 1 and we deal with equation (9).

Merton Process and Euler scheme
In this section we deal with the Merton model: where ρ, b and c satisfy the previous hypothesis 4.1, 4.2 and 4.3. Moreover, we assume The function x → σ(u, x) is twice continuously differentiable and there exists two constants C > 0 and > 0 such that: We will present two alternative calculus for this model. The first one is based on the Brownian motion only and the second one is based on both the Brownian motion and the jump amplitudes. Suppose that the jump times T 1 < . . . < T n are given (this means that we have already simulated T 1 , . . . , T n in the Monte-Carlo algorithm). We include them in the discretization grid: so we consider a time grid 0 = t 0 < t 1 < . . . < t m < . . . < t M = T and we assume that T i = t m i , i = 1, . . . , n for some m 1 < . . . < m n . For t > 0, we denote m(t) = m if t m ≤ t < t m+1 . Then the corresponding Euler scheme is given bŷ This corresponds to the deterministic equation : where we have denoted by Moreover we have : The first derivatives ofŝ t satisfy the following equations : For the derivatives of higher order one may derive similar equations. Now we have the choice of using the integration by parts formula from the section 2, using ∆ i W or both ∆ i and ∆ i W . And in each case we have different forms for the differential operators. For example on the set The other differential operators will change in a similar way. Note that ∆ i W ∼ N (0, t i+1 − t i ) so that the corresponding L W i is given by Then the Ornstein-Uhlenbeck operator will be This allows us to use the integration by parts formula corresponding to the Brownian motion only, or to both Brownian motion and jump amplitudes (the first case leads to the same calculus as in [6] and [13]). It is more delicate to prove the non degeneracy condition (8) if we only use the jumps.

Numerical results
We compute here the Delta of two European and Asian options (call option with payoff φ(x) = (x − K) + and digital option with payoff φ(x) = 1 x≥K ) with underlying (S t ) t≥0 following a jump type diffusion model. Let N be a compound Poisson process with jump intensity λ. For all i ∈ N * , we denote ∆ i the jump amplitude of N at the jump time T i . We suppose that ∆ i ∼ N (0, 1), i ≥ 1.
We deal with two different jump diffusion models. The first one is motivated by Vasicek's model for interest rates (but we consider a jump process instead of a Brownian motion): The second one is of Black-Scholes type: We compare the different Monte Carlo estimators of the Delta, obtained by using

Localization method
For European and Asian call options, we use the same variance reduction method as the one introduced in [9]. For European options, sensitivity analysis using Malliavin calculus leads to terms such as φ(S T ), H n (S T , ∂S T ∂β ) (take I T for S T in the case of Asian options), which may have a large variance. It is possible to avoid this problem by using a localizing function which vanishes out of an interval [K − δ , K + δ], for some δ > 0. Let us introduce some notations.
Note that B δ is the derivative of G δ . We define Since F δ vanishes out of [K − δ, K + δ], the value of the second expectation does not blow up as H n increases.

Finite Difference method
Arbitrage theory gives an expression for the price u( , ) of an option, with underlying S and payoff φ, as the following expected value : To compute the Delta, the finite difference method makes a differentiation using the Taylor expansion of the price with respect to S 0 . Indeed, we shift S 0 with and compute the new price u(0, S 0 + ), then the first term of the Taylor expansion of the price around S 0 is given by: We choose the symmetric estimator and we use the same simulated paths in the two "shifted expectation" in order to reduce the variance.

Malliavin estimator
We recall that we have computed the Delta estimator using Malliavin integration by parts formula (see section 2): We compute here the Malliavin weights H i J T (S i T , ∂ β S i T ).
• We first study the diffusion process defined by (25). We have an explicit expression of S T on {J T = n}: In order to compute the Malliavin derivatives of S T involved in the weight H J T , we differentiate with respect to the jump amplitudes in (27). Thus, for all 1 ≤ i ≤ n, one has and the covariance matrix is defined by : .
Since ∂ ln ρ(∆) ∂∆ = −∆, one has Finally, using equation (10), • We study the jump diffusion process defined by (26). On {J T = n}, we have so, for all 1 ≤ i ≤ n, Notice that if (1 + σ ∆ i ) = 0, then S T = 0. So we employ the convention 0 0 = 0. This is just a matter of notation. Let us define then we get, for all 1 ≤ i ≤ n

Hence, for European options
For Asian options, on {J T = n} we have with the convention T 0 = 0 and T n+1 = T . For all t ∈ [T j , T j+1 [, r S u du, so S t = S T j e r (t−T j ) and we obtain Then So for all 1 ≤ i ≤ n, we have Then we get Hence, for Asian options where K = n j=1 K 2 j,T and K j,T = After simulating trajectories of the specified diffusion, we can compute the Malliavin weights by (28), (33) and (35). Note that we only need to know the amplitudes and the times of the jumps to state a Monte-Carlo estimator of the Delta. Remark 6.2 For this model, ∂ a c(t, a, x) = σ x so hypothesis 4.3 is not satisfied. But, using the localization method (that is S T ∈ [K − δ, K + δ] and so I T ∈ [K − δ, K + δ]), the non degeneracy condition (8) still holds true: (32) and (34) give

Merton Process
Here we focus on the Merton model (2): where W is a Brownian motion independent on the compound Poisson process N . We suppose that the jump amplitudes ∆ i are independent, identically and Gaussian distributed. We have the explicit solution: so we do not need to use Euler's scheme.
The covariance matrix is

Straightforward computations give
where A µ is given by (29). We put these terms together to get the Malliavin weight: where B µ and C µ are defined by (30) and (31) respectively.
Recently in [6] and [13], the Delta of an European option is computed by using the Malliavin calculus with respect to the Brownian motion only. Notice that if we use our integration by parts formula, but we just take into account the derivatives with respect to the Brownian motion, we find H = W T S 0 σ T , which is exactly the weight obtained in [13] (as well as in Black-Scholes). So the difference between our algorithm and the one of [13] comes from the supplementary part (corresponding to the derivatives with respect to the jump amplitudes) which appears in our Malliavin weight H. In figure 2, we compare the two algorithms. Moreover, in table [1], we give the quotient between the empirical variances of the two algorithms. It turns out that the variance of the Brownian-jump algorithm (presented here) is smaller than the variance of the pure Brownian algorithm (presented in [13]). Moreover, the difference becomes more important as the number of jumps up to T increases: this happens when the maturity T is larger or when the intensity λ of the Poisson measure is larger. We conclude that more noise one uses in the integration by parts formula, better the algorithm works (there is no theoretical result in this sense, but only numerical evidence).       These figures confirm that we can numerically compute the Greeks for European options with a pure jump underlying process. We obtain numerical results similar to those in the Wiener case ( [9] and [8]). On the one hand, for European and Asian Call options, the Malliavin estimator has more variance than the finite difference one (see figures 4, 6 and 7): the finite difference method approximates the first derivative of the payoff, whereas the Malliavin estimator contains a weight independent on the payoff, which may increase the variance. Using a localization technique permits to reduce the variance of the Malliavin estimator and to make it closer to the finite difference one. On the other hand, the Malliavin estimator of a Digital option has less variance than the finite difference one (see figures 3 and 5) and so does not need to be localized: in this case, the first derivative of the payoff is a Dirac and, contrary to the finite difference method, the Malliavin calculus permits to avoid this strong discontinuity. Finally, notice that for both call and digital options, the finite difference method requires to simulate twice more samples of the asset than the Malliavin method does : the finite difference method uses the samples starting from S 0 and those starting from S 0 + . So the Malliavin method is less time consuming.

Conclusion
In order to develop a Malliavin calculus for jump processes, we worked here on a space of simple functionals of a finite set of random variables representing the source of randomness. First, we applied the Malliavin calculus to pure jump processes, by differentiating with respect to the jump amplitudes. Second, in the case of Merton process, we used the total randomness coming from the jump amplitudes and the Brownian increments. Then, for the sensitivity analysis, we have set up (under some non degeneracy conditions) an integration by parts formula, which permits to replace the derivative of the payoff by a Malliavin weight.
The Numerical experiments show that using Malliavin approach becomes extremely efficient for a discontinuous payoff. Some specific techniques can be used to reduce the variance of the Malliavin estimator. For the Merton process, we compare our Delta estimator with the one where one omits the presence of the jump process (that is we do not differentiate with respect to the jump amplitudes). The results show that it is better to use the two sources of randomness especially when we have more jumps.