Bismut–Elworthy–Li formula for subordinated Brownian motion applied to hedging financial derivatives

Abstract The objective of the paper is to extend the results in Fournié, Lasry, Lions, Lebuchoux, and Touzi (1999), Cass and Fritz (2007) for continuous processes to jump processes based on the Bismut–Elworthy–Li (BEL) formula in Elworthy and Li (1994). We construct a jump process using a subordinated Brownian motion where the subordinator is an inverse -stable process with . The results are derived using Malliavin integration by parts formula. We derive representation formulas for computing financial Greeks and show that in the event when , we retrieve the results in Fournié et al. (1999). The purpose is to by-pass the derivative of an (irregular) pay-off function in a jump-type market by introducing a weight term in form of an integral with respect to subordinated Brownian motion. Using MonteCarlo techniques, we estimate financial Greeks for a digital option and show that the BEL formula still performs better for a discontinuous pay-off in a jump asset model setting and that the finite-difference methods are better for continuous pay-offs in a similar setting. In summary, the motivation and contribution of this paper demonstrates that the Malliavin integration by parts representation formula holds for subordinated Brownian motion and, this representation is useful in developing simple and tractable hedging strategies (the Greeks) in jump-type derivatives market as opposed to more complex jump models.


PUBLIC INTEREST STATEMENT
This paper entitled Bismut-Elworthy-Li Formula for Subordinated Brownian Motion Applied to Hedging Financial Derivatives provides pricing and risk management methods usable by quantitative risk analysts. It provides knowledge on how one can design and determine prices of new derivative products in a jump-type financial market. In addition, market risk can be managed using Greek representations to determine for instance, proportions in one's investment portfolio, make informed decision based on the behaviour of volatility and/or potential extreme shocks in the market. It can be argued that simpler methods could be used instead to achieve the above objectives. Agreed, however, care must be taken to ensure that simple model of choice represents the market as closely as possible. This current paper is geared towards achieving exactly that. Stochastic volatility models with similar modelled jumps could be further innovations.

Introduction
The problem of computating the Greeks of derivatives with smooth pay-off functions has been extensively studied. The problem where the pay-off function is not necessarily regular poses a different level of difficult and requires a different approach. See for instance Fournié, Lasry, Lions, Lebuchoux, and Touzi (1999) and Cass and Fritz (2007). Notice that existing techniques avoid differentiating the pay-off function by introducing a weight function. The Bismut-Elworthy-Li (BEL) representation formula (Elworthy & Li, 1994) is one scenario of such innovations. In this paper we show that the known relationship between the Malliavin derivative and the first variation process still holds for an alphastable subordinated Brownian motion and results in an explicit martingale weight factor. This allows for an extension of the BEL formula to subordinated Brownian motion and as a result, Greeks can be easily computed in jump-type and/or emerging markets. The subordinator belongs to the Lévy family of four parameter alpha-stable distributions. Price dynamics of almost all financial instruments in the different financial markets are observed to deviate from the Gaussian distribution. Various models in literature have been developed to closely estimate the dynamics of these markets. The rich and robust family of alpha-stable distributions has proven successful over most traditional models in capturing skewed and heavy tailed distributions. An application to estimate the densities of solutions to stochastic differential equation with subordinated Brownian motion under Malliavin framework is discussed in Kusuoka (2010). A rather different approach is discussed in Wyłomańska (2012) where the authors, in addition to investigating the densities of subordinated Brownian motion, also discuss some properties related to transforms and averaged mean squared displacements of the process. They consider and compare both cases of stable processes and their inverses. In addition they provide some parameter estimation techniques. An intuitive study related to the work in this current paper is Elworthy and Li (1994). However, this is limited to the Delta. The authors derived derivatives of solutions of diffusion equations and demonstrate that these derivatives exhibit and allow for the estimation of a diffusion equations' smoothing properties. In addition, they use their results to study the logarithmic gradient of heat kernels. Their results can be extended to deriving representation formulas for more Greeks. A less related but still interesting work is by Song and Vondraček (2003) where properties of a killed subordinated Brownian motion by an ∕2-stable are compared with those of an -stable subordinated Brownian motion. They show possible comparability in their killing measures and propose bounds on the Green function and a jumping kernel of a subordinated ( ∕2) process. The work by Zhang (2012) on the derivative and gradient estimate for SDEs driven by stable processes is another intuitive literature vital to our study. Other related works that are recent include Takeuchi (2010), Khedher (2012) and Kawai and Kohatsu-Higa (2010). There is more existing literature on stable distributions that we cannot exhaustively discuss. Interested readers can refer to Sun and Xie (2014) for a comprehensive literature on stable distributions.
In summary, the motivation and contribution of this paper demonstrates that Malliavin integration by parts representation formula holds for subordinated Brownian motion and, this representation is useful in developing simple and tractable hedging strategies (the Greeks) in jump-type derivatives market as opposed to more complex jump models.
The rest of the paper is organised as follows: In Section 2 we review some traditional techniques for computing the Greeks. In Section 3 we derive some useful results about stable processes and Malliavin Calculus with regards to subordinated Brownian motion. In Section 4 we analyse the differential calculus of our choice of process and show that the integration by parts formula exists. In Section 5 we present the BEL formula for subordinated Brownian motion and use the result to derive the main results of the paper. In Section 6 we discuss some applications. Section 7 concludes.

Sensitivity analysis
The theory of risk-neutral valuation asserts that given a complete filtered probability space (Ω,  , ( t ) t≥0 , ) with a subjective probability measure and a filtration ( t ) t≥0 satisfying the usual conditions, we can construct a pay-off function Φ of an option under an equivalent martingale measure whose price is given by where Φ = e −rT Φ is a discounted pay-off function by a risk free interest rate, r and X T is the value of the underlying at maturity time T. A Greek is a derivative of (1) with respect to a certain model parameter p which could be the initial value X 0 of the underlying, the volatility parameter , time to maturity = T − t, the option strike, E or the interest rate, r i.e.
Clearly, (2) poses a problem if Φ (X T ) is not differentiable. Different techniques in literature have been explored to resolve this challenge. For instance the Likelihood method (Fournié et al., 1999) suitable for known distribution of the underlying. It takes the form where denotes the density of the underlying. Another widely known approach is the method of Malliavin calculus (see Bavouzet & Messaoud, 2006;Bayazit & Nolder, 2013;Bermin, 2000;Di Nunno, Øksendal, & Proske, 2008;Nualart, 1995). It eliminates differentiation of the pay-off function by introducing a weight term composed of Malliavin derivative and the Ornstein-Uhlenbeck operator terms, i.e.
where the weight factor  consists of Malliavin derivatives of random variables X and G belonging to some space with nice properties e.g. L 2 . This approach is more flexible than the previous one in the sense that the distribution of the underlying is irrelevant in computing the Greeks. However, it is computationally expensive. Another approach is the Bismut-Elworthy-Li representation formula Elworthy and Li (1994), Cass and Fritz (2007) and Baños, Duedahl, Meyer-Brandis, and Proseke (2015), where p = x, the underlying initial spot price and a t is some bounded function satisfying The Bismut-Elworthy-Li Formula (5) applies to continuous diffusion processes but can be adapted to finite (see Cass & Fritz, 2007) and infinite (see Chen, Song, & Zhang, 2015) jump processes. The usefulness of this formula is its allowance for an explicit representation of the Delta of a financial derivative. (1) Greek = V t p . (3) T ∫ 0 a s ds = 1 e.g. a = 1 T .
For instance by employing the theory of Malliavin calculus, the weight  can be obtained explicitly for different Brownian motion functionals. Malliavin calculus for both continuous and jump diffusion processes has been extensively discussed in literature and there exist enormous applications on the subject (see Baños et al., 2015;Bavouzet & Messaoud, 2006, Bayazit & Nolder, 2013Cass & Fritz, 2007;Di Nunno et al., 2008;Elworthy & Li, 1994;Kusuoka, 2010). The focus of this paper is to compute the Greeks for a wide range of pay-offs irrespective of their structure by employing the BEL formula in the framework of subordinated Brownian motion by ∕2-stable distributions, with ∈ (0, 2).

Stable distributions
In this section we shall adapt some of the definitions and notations from Kateregga, Mataramvura, and Taylor (2017).
Definition 3.1 A stable distribution is a four-parameter family of distributions denoted by S( , , , ) where ∈ (0, 2] is the characteristic exponent responsible for the shape of the distribution.
> 0 is the scale parameter (it narrows/extends the distribution around ).
∈ ℝ is the location parameter (it shifts the distribution to the left or the right).
In this paper we are only interested in positive, non-decreasing, stable cádlág processes denoted by (S t ) t≥0 whose density functions are defined on ℝ + . For simplicity, we shall denote by  the set of all such processes. More so, we shall denote the inverse of S t by L t , t ≥ 0. The purpose is to introduce jumps in asset returns through subordinated Brownian motion and in turn, provide pricing and hedging (the Greeks) representation formulas.
Essentially, densities of stable processes are estimated using characteristic functions through Fourier transforms. 1 The characteristic function of S t is obtained using its domain of attraction and the Lévy-Khinchine representation formula (see Applebaum, 2004) and it's given by Consequently, the density of S t is given by Meerschaert & Straka, 2013): For l ∈ [0, T], it is easy to see that the following equivalence relation holds: The process L t is also interpreted as the first passage time of S. Moreover Let h L (l;t) denote the density function of L t . Using relation (10) we deduce We can therefore approximate h L (l;t) by estimating the integral in (12) using the density h S (t;l) in (8) and its characteristic function (7).

Subordination
Let B t t≥0 , denote standard Brownian motion and (L t , ∈ (0, 1]) be an inverse stable process independent of B. Then the process B L t is referred to as stable subordinated Brownian motion.
starting at zero. We define a joint product probability space where Ω = ([0, L T ]) endowed with the natural filtrations: We further introduce a separable Hilbert space defined by to obtain a complete abstract probability space (, Ω, , ( t ) t≥0 ).
The first and second moments of subordinated Brownian motion follow by applying semigroup properties. Literature on semigroups and infinitesimal generators is detailed in Song and Vondraček (2003) and Applebaum (2004).
The transition density p(x, y; t) of B is given by and its semigroup (P t ) t≥0 given by where f is a nonnegative Borel function on ℝ.
The operator (Q t ) t≥0 has a transition density q(x, y; t) defined by Lemma 3.6 The mean and variance of B L t are computed by Proof Suppose f in Definitions 3.4 and 3.5 is such that f (z) = z. Using (17) and (18) we immediately observe that Equation (21) follows immediately by definition of variance of B L l . That is ✷ Note that for ∈ (0, 1], [L t ] exists and can be computed. If L t ≡ t in (20) and (21) we recover standard Brownian motion with mean 0 and variance t. In general, the k-th moment of L t is given by Lemma 3.7 The covariance of B L t is given by Since for all t ∈ ℝ + , B L t has independent increments with zero mean, we have Similarly for L t ≤ L s we have the covariance as [L t ]. Then, we write .

Malliavin derivative in the direction of jump processes
In this section we explore the differential calculus of B L t using directional derivatives. Little has been done on this, a few references include Kusuoka (2010) and Zhang (2012).
Lemma 4.1 Following Definition 3.3, we let f be an ( t )-adapted right-continuous process with left limits satisfying Then its ( t )-martingale stochastic integral exists, it's well defined and can be expressed as an where S t is the inverse stable process of L t .
Proof This follows from the standard change of time. ✷ Following Definition 3.3, we denote by D the Malliavin derivative operator defined on  such that ḣ represents differentiation of h with respect to L t . In addition, D h will denote Malliavin differentiation in direction h.
be subordinated Brownian motion associated with (, Ω, ( t ) t≥0 , ). Then For the second part of the lemma we notice that since L t is of bounded variation, the contribution of its small jumps is almost negligible and the number of jumps is finite. We partition [0, T] The result follows immediately, i.e. ✷

Lemma 4.3 Suppose f is a right continuous function with left limits, then
where S t is the inverse stable process of L t . Proof This follows a similar argument as in the previous Lemma. ✷

Discrete multiple stochastic integral
In this section we derive the integration by parts formula associated with our choice of random process B L t . We denote by k the time between the k-th and (k + 1)-th jumps of B L t .
We start with the triplet (, Ω, ), a joint probability space introduced in Definition 3.3 with a real seperable Hilbert space , with a scalar product ⟨⋅, ⋅⟩  and the norm for an element g ∈ , denoted by ‖g‖  , Ω is the completion of  and is the extension to the Borel −algebra of Ω of a cylindrical measure. We define independent stable random variables k : = B L k − B L k−1 which are canonical projections from Ω to ℝ. We assume L k is a form of subordinator introduced in Janicki and Weron (1993, p. 33) which is an ∕2-stable totally skewed Lévy motion with increasing sample paths ( ∈ (0, 2), = 1). This is a symmetric alpha-stable process (S S) with positive Poisson jumps. Therefore, B L k belongs to a class of S S Lévy motion processes with its jumps only at the jump times of L k . As a consequence, we use Charlier polynomials to define the multiple stochastic integrals with respect to our process.
Definition 4.5 The Charlier polynomials are defined as They can also be expressed explicitly as The Charlier polynomials form an orthogonal basis of  2 (Ω,  , ) with respect to the Poisson measure (dx) = x x! e dx. Moreover, we have: where nm = 0 when n ≠ m and nm = 1 for n = m. Therefore any function F ∈  2 (Ω,  , ) can be uniquely represented as with the corresponding norm given by ‖f ‖ 2 = ∑ n≥0 �f n � 2 n n!.
Next, we construct the discrete multiple stochastic integral using the Wick product following similar arguments in Privault (1990).
Suppose  * is a set of all functionals of the form Q( 0 , … , n−1 ) where Q is a real polynomial and n ∈ ℕ, we regard  * as an algebra generated by C ( ) n ( k ):k, n ∈ ℕ and define the Wick product in the following manner.
Definition 4.6 The Wick product of two elements F, G ∈  * denoted by F ⋄ G is defined (relaxing for simplicity) as: where for a ∈ ℕ d , a! = a 1 ! … a d ! and n = (n 1 , … , n d ), m = (m 1 , … , m d ) and k 1 ≠ ⋯ ≠ k d .
Let  = l 2 (ℕ) be a space of square-summable sequences. There exists a discrete chaotic decomposition of  2 (Ω, ) whose elements F can each be represented as a sum of multiple stochastic integrals of kernels of  •n = l 2 (ℕ) •n . That is, .
where f n ∈  •n , n ∈ ℕ and I n (f n ) is the discrete multiple stochastic integral of symmetric functions of discrete variable.
The stochastic integral of f ∈ l 2 (ℕ) is an isometry from  = l 2 (ℕ) to  2 (Ω,  , ) as follows: Proposition 4.7 Let f ∈ l 2 (ℕ) and define I 1 (f ): = Proof Consider a partition 1 < ⋯ < m−1 where i , i = 1, … , m − 1 are the jump times of B L t and let We notice that for j 1 < j 2 on each [ i , i+1 ), we have which is arrived at by conditioning with respect to  t i, j 2 and applying the tower property. Meanwhile for j 1 = j 2 = j we have The last equation follows from the law of total expectation. The result follows immediately by combining both cases. ✷ The discrete multiple stochastic integral I n (f n ), f n symmetric in l 2 (ℕ n ) with finite support can be defined directly using the Wick product.
Definition 4.8 The symmetric tensor product f 1 • ⋯ •f n is defined as where f 1 , … , f n ∈  and Σ n is the set of all permutations of {1, … , n}. Moreover, suppose g 1 , … , g n ∈ l 2 (ℕ) with finite supports, we have .
The definitions above suggest the results (Privault, 1990, Prop. 2 and 3) for the Poisson process, also hold for our choice of process B L t with similar proofs.
⟨I n (f n ), The annihilation operator defined in (26) has an equivalent in the discrete chaotic decomposition of  2 (Ω,  , ) given by Moreover, the following lemma holds.
Lemma 4.11 Suppose  denotes a dense set of elements of u ∈  2 (Ω) ⊗ l 2 (ℕ) such that u k = k h k , k ∈ ℕ where h:ℕ → Q has finite support in ℕ and define an operator : →  2 (Ω) by Then for any F ∈ Dom(D) and u ∈ Dom( ), we have where Dom represents domain and the operators D and are also closable and adjoint to each other.
Proof See Privault (1990). ✷ The above results can be extended  2 (Ω) ⊗ l 2 (ℕ) to  2 (Ω) ⊗  2 (ℝ + ) yielding: This in turn leads to the following duality relation: The duality formula is an important tool for finding alternative representations of derivatives of expectations of irregular functions. We discuss this later.
Note that for the rest of the article we shall use abbreviations SDE and SSDE to represent stochastic differential equation and subordinated stochastic differential equation respectively.

Malliavin derivative of solutions to subordinated SDEs
In the following, K is defined as a Hilbert space, D 1, 2 (K) as a Sobolev space of K-valued functions associated with the H-derivative,  2 (H; K) as a total set of a K-valued linear operator of Hilbert-Schmidt class on H,  1, 2 (dB L t ;K) as the total set of ( t )-predictable (ℝ × K)-valued functions such that for (t, X) ∈ D 1, 2 (ℝ × K), t ∈ [0, T] with where h is given in (15). We denote (dt;K) as the total set of ( t )-predictable (ℝ × K)-valued func- Let (Ω,  , ) represent a joint probability space introduced in the previous sections and consider the following stochastic differential equation: A scenario of (14) with Lipschitz coefficients and standard Brownian motion has been discussed in Fournié et al. (1999) and Bayazit and Nolder (2013) to compute financial Greeks. The case of non-Lipschitz coefficients with subordinated Brownian motion is handled in Bavouzet and Messaoud (2006), Sun and Xie (2014) and much less related in Cass and Fritz (2007) and Di Nunno et al. (2008). For non-Lipschitz with standard Brownian motion, (see for instance Baños et al., 2015). In the current paper we are interested in a model with Lipschitz coefficients and subordinated Brownian motion.
Proposition 4.12 Suppose in the stochastic differential equation (42) has Lipschitz coefficients and assume (t, X) ∈ D 1, 2 (ℝ × K) and b(t, X) ∈ D 1, 2 (ℝ × K). Then the solution X t to (42)  Then apply the product rule and (26) to the second term on the RHS of (46) to obtain (43).
If h is an orthogonal basis, we can express as where h = ‖h‖ĥ and ĥ is a unit vector of h. For h orthonormal gives (44). ✷ If ≡ 1 the last term of (45) is zero and by the Grönwall's inequality, we have: Also note from (46) that the first variation process can be deduced as: Combining (48) and (49) results into the following useful relation Alternatively,

BEL formula for subordinated stochastic differential equations
Bismut-Elworthy-Li formula for general Lévy processes is studied in Cass and Fritz (2007). We derive representations for subordinated Brownian motion based on (Elworthy & Li, 1994).
Proposition 5.1 Let X t be the solution to (42) on the horizon [0, T] where L t ≡ t (see for instance Baños et al., 2015;Cass & Fritz, 2007) and let Φ:ℝ → ℝ denote some bounded function. Suppose we can define a functional V t (X T ) of X T by Then the derivative of V with respect to x is given by where a t is some bounded function satisfying Proof Apply the classical chain rule on [ x Φ(X T )] and use the relation (51) followed by the chain rule in the Malliavin sense. Finally apply the duality relation (39), in that order. A similar proof is provided in Sturm (2004) using the identity D t X T = J T J −1 (53) is Bismut-Elworthy-Li formula for Geometric Brownian motion. (48) T ∫ 0 a s ds = 1.
Proposition 5.2 Suppose (46) has Lipschitz coefficients. Let R(t, X t ) denote the right inverse of (t, X t ) where (t, X t ) is elliptic. For any function Φ ∈ C 1 b (ℝ) and h ∈ ℝ, we have (the x argument is relaxed for simplicity) where Proof Working backwards and using the results obtained above we have Next we apply the duality relation (39), Equation (50), Grönwall's inequality and (45) for some arbitrary h not necessarily an orthonormal basis and ≥ 1 in that order, where we have chosen We provide a detailed analysis on the above result in the following section.

Main results
This section presents the main results of the paper. The idea is to extend the results by Fournié et al. (1999) to a subordinated stochastic differential equation model by deriving the first-and secondorder derivative representation for the expectation of a function that is not necessarily regular. Specifically, the idea is to by-pass the derivative of the expected (irregular) function by introducing a weight term in form of an integral with respect to subordinated Brownian motion. The results in this section are employed in the following section to estimate the Greeks using MonteCarlo simulations.
In this section, operators and D will be used interchangeably to represent weak derivatives, and J t shall denote the first variation process given by  shall denote a space of bounded integrable functions, (Q t :t ≥ 0) shall denote the semigroup of the solution X t to (46). Lastly, S t t≥0 shall denote a non-decreasing cádlág -stable process and L t t≥0 its inverse with ∈ (0, 1]. We require X t to be complete to enforce some integrability conditions on DX t . Corollary 5.3 Let U:ℝ → (ℝ) with bounded first derivative such that Q t (U):ℝ → (ℝ) and define U X : = (df ) X = Df (X t ). Then the weak derivative of Q t with respect to the initial state x is given by provided the last term exists. Moreover Since X t is non-explosive, the weak derivative Q t is well defined.
Proof This follows directly from the application of a weak chain rule. ✷ Corollary 5.4 Suppose X t ∈ ℝ is non-degenerate and elliptic, there exist an inverse R(t, X t ) of (t, Proof Recall that the result holds for the case of continuous processes (Cass and Fritz, 2007;Elworthy & Li, 1994). We can therefore employ similar arguments of partitioning as in the second part of the Proof of Lemma 4.2 and apply similar steps as in the continuous case but piece-wise, to arrive at the required result. ✷ Theorem 5.5 Let Φ:ℝ → ℝ with its first derivative bounded and continuous: Moreover, for x ∈ ℝ, T > 0, the derivative with respect to x is given by Proof From Corollary 5.4, let T > 0, applying Itô's formula to yields As t → T, and applying the knowledge from Lemmas 4.3 and 4.4 yields Multiplying (66) by a martingale . The required result follows. ✷ Corollary 5.6 Suppose X t ∈ ℝ 2 , t ≥ 0 and indexes 0 ≤ j, k ≤ m, then and and where | ⋅ | denotes the Euclidean norm.
Proof The proof follows directly by applying (5) where h ≡ J 0 . ✷ Theorem 5.7 From Corollaries 5.4 and 5.6, let Φ:ℝ 2 → ℝ be such that its first and second derivatives are bounded and continuous and, Proof From Equation (72) we deduce Suppose L T = L t ∕2 and 0 ≤ ≤ t∕2 then Applying Itô formula to Q t− Φ(X ):0 ≤ L < L t at L = L t ∕2, yields R(S , x)(J j 0 , J k 0 )dB . Next, taking expectations and applying the identity (see Elworthy & Li, 1994) yields the required result. ✷

Applications
This section is dedicated to estimating the Delta, Gamma and Vega from two stochastic models of the asset price namely; the subordinated stochastic differential equation (SSDE) and the Geometric Brownian motion (GBm).

Hedging discontinuous-pay-off type options in Black-Scholes framework
We focus only on the digital option but the analysis could be extended to other discontinuous-payoff type or irregular pay-off options (see Fournié et al., 1999 for instance). The Delta and Gamma follow directly from Theorems 5.5 and 5.7 respectively.
Let X t satisfy the stochastic differential equation with the solution given by Figure 1 shows the dynamics of solutions to SSDE and GBm.
The Greeks are computed by conditioning on L T = , as follows :

Delta
According to Theorem 5.5 Assume a discounted pay-off of a digital option i.e. Φ(X T ) = e −rT

E<X T
where E is the strike price, then we can express the Delta by Figure 2 shows the digital option delta from both the SSDE and GBm models. Observe that the delta from SSDE is slightly higher due to existence of jumps and its convergence is slightly slower. ] and simplying. That is Figure 3 shows digital option gamma from the SSDE and GBm models. Again we observe that the gamma from SSDE is slightly higher than that from GBm.

Vega
The Vega can be deduced similarly using integration by parts. That is Figure 4 shows the Vega from the SSDE and GBm models. Note that the Vega from SSDE is slightly higher. Observe that despite fact the SSDE model has jumps, the convergence rate for the estimation of the Greeks from the model is as good as in the GBm model.  As a matter of interest, we apply the finite-difference method on the subordinated Brownian motion model to estimate the Greeks for a call option from the SSDE model. Recall from Fournié et al. (1999) that the finite-difference method is recommended for computing the Greeks from European options compared to the BEL formula, it performs better in this case. Figure 5 shows the estimation of the Greeks of a European call option using the SSDE model.