Computation of Greeks and Multidimensional Density Estimation for Asset Price Models with Time-Changed Brownian Motion

Abstract The main purpose of this article is to propose computational methods for Greeks and the multidimensional density estimation for an asset price dynamics model defined with time-changed Brownian motions. Our approach is based on an application of the Malliavin integration-by-parts formula on the Gaussian space conditioning on the jump component. Some numerical examples are presented to illustrate the effectiveness of our results.


Introduction
The Malliavin integration-by-parts formula was established in the historic paper of Bismut (1981), and served as a key tool to reveal the relation between the Hörmander condition on the hypoelliptic problem and the existence of smooth densities for solutions to stochastic differential equations. Moreover, in his book (Bismut, 1984), he studied the logarithmic gradient of the fundamental solution to the heat equation on a Riemaniann manifold. In mathematical finance, Malliavin calculus on the Wiener space is applied in Fournie et al. (1999) to the sensitivity analysis for the Black-Scholes asset price dynamics model. Moreover, the so-called Malliavin-Thalmaier integration-by-parts formula was introduced in Malliavin and Thalmaier (2006) with a view towards the computation of multidimensional probability density functions, of which numerical issues are studied in Kohatsu-Higa and Yasuda (2009). For stochastic differential equations driven by Brownian motion, the technique based on the integration-by-parts formula from Malliavin calculus has become standard.
It is then a natural question whether a similar approach can be obtained in the case of computation of Greeks and the multidimensional density for jump processes. Concerning Malliavin calculus for jump processes, the existence of weights for the logarithmic gradient dates back to the work of Bichteler et al. (1987). For practical considerations, it is important to derive explicit weights, in order to design an efficient Monte Carlo evaluation. In particular, Davis and Johansson (2006) and Cass and Friz (2007) studied the Malliavin Greeks for jump diffusion processes. Their approach consists of conditioning by the jump component and then performing the Malliavin calculus techniques with respect to the Brownian motion component. In other words, their models are required to be a superposition of independent jump and diffusion components, where the diffusion must be nondegenerate. On the other hand, El- Khatib and Privault (2004) applied the calculus focused on the Poisson arrival times of Carlen and Pardoux (1990), while Bally et al. (2007) took a unified approach considering the derivatives with respect to both the Poisson arrival times and the amplitudes of the jumps. Takeuchi studied in 2009 the same problem of solutions to stochastic differential equations with jumps via the martingale approach, without using the approach introduced in Bichteler et al. (1987). Takeuchi (2008, 2009) derive Greeks formulas for asset price dynamics described by gamma processes, based upon a scaling property of gamma processes with respect to the Esscher transform parameter. It seems that there is still no unified approach to the sensitivity analysis and also the multidimensional density estimation for random processes beyond the pure diffusion cases.
The main purpose of this article is to obtain an expression for the calculation of Greek and the multidimensional density for an asset price dynamics model defined with time-changed Brownian motions. As will be illustrated later, our framework encompasses some Le vy process models and stochastic volatility models. Our approach is based on an application of the Malliavin integration-by-parts formula on the Gaussian space conditioning on the jump component. We consider this approach best suited for actual financial models as we reflect in the examples that we have considered. This approach has been taken in Cass and Friz (2007) and Davis and Johansson (2006). Nevertheless they use the representation of the model as a stochastic differential equation while we use an expression as subordinated Wiener processes. In the Le vy processes framework with an expression as a time-changed Brownian motion, our results improve those of Bally et al. (2007) and El-Khatib and Privault (2004) in the sense that our Le vy processes are of more realistic infinite activity type. When restricted to the variance gamma model, our results can be recovered from those of Takeuchi (2008, 2009).
The rest of this article is organized as follows. We first introduce in Section 2 the timechanged Brownian motion and give three examples of time-changed Brownian motions of practical interest. In Section 3, we apply the Malliavin integration-by-parts formula on the Gaussian space to derive unbiased estimators of various Greeks in our model setting with a view towards efficient Monte Carlo evaluation. We next present in Section 4 that the Malliavin-Thalmaier integration-by-parts formula can also be applied to the multidimensional density estimation in our framework. Finally, Section 5 concludes.

Preliminaries
Let us begin with general notations which will be used throughout the text. We let N be the set of positive integers. We denote by R d the d-dimensional Euclidean space with the norm Á k k and the inner product by hÁ; Ái. We fix ð; f; PÞ as our underlying Computation of Greeks for Time-Changed Brownian Motion 303 probability space. We mean by the time-changed Brownian motion the following stochastic process where 2 R, 2 R, > 0 and fW t : t ! 0g is a standard Brownian motion in R and fY t : t ! 0g is a non-decreasing stochastic process in R, with Y 0 ¼ 0, a.s., independent of fW t : t ! 0g. For convenience, we define the characteristic exponent e ' Y t ðyÞ :¼ E e iyY t Â Ã : Intuitively, one may regard the original clock ft : t ! 0g as calendar time, while the new random clock fY t : t ! 0g as business time. A more active business day implies a faster business clock and vice versa.

Examples of Time-Changed Brownian Motion
If fY t : t ! 0g is a subordinator, that is, an one-sided Le vy process, the above timechanged Brownian motion is a Le vy process, which is often called Le vy process of type G (see, e.g., Rosin ski, 1991). Moreover, unless it is a deterministic subordinator (that is, Y t ¼ at, a > 0), the time-changed Brownian motion becomes a pure-jump Le vy process. In Examples 2.1 and 2.2, we give two such examples.
Example 2.1 (Variance gamma Le vy processes) The stochastic process (2.1) is called a variance gamma (Le vy) process of Madan et al. (1998), if we set fY t : t ! 0g to be a gamma process, whose characteristic exponent is given by ' Y t ðyÞ ¼ Àat ln 1 À iy=b ð Þ , while its marginal density is given in closed form by The model has attracted the attention of finance practitioners, and thus appears often in the computational finance literature (see, e.g., Carr and Madan (1999) and Fu (2007)). The characteristic function of Moreover, it is also shown in Madan et al. (1998) that the variance gamma process can be expressed as the difference of two independent gamma processes, where the gamma processes fY t;p : t ! 0g and fY t;n : t ! 0g can be characterized by Equation (2) with ða; bÞ ¼ ð À1 ; ð p Þ À1 Þ and ða; bÞ ¼ ð À1 ; ð n Þ À1 Þ for some > 0, where the parameters are given by p ¼ 1 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 þ 2 2 = p þ =2 and n ¼ 1 respectively. & Example 2.2 (Normal inverse Gaussian Le vy processes) The normal inverse Gaussian process can be expressed as a time-changed Brownian motion by setting fY t : t ! 0g to be a Le vy process induced by the inverse Gaussian distribution, whose marginal at time t is identical in law to the first time that a Brownian motion with drift v reaches the positive level t, that is, Then, the normal inverse Gaussian process See Barndorff-Nielsen (1998), Carr et al. (2003) and Schoutens (2003) for more details. & It has been a general consensus that return volatilities vary stochastically over time. The concept of 'time change' can also be applied to capture the evidence of stochastic volatility. That is, by the Gaussian scaling property, random changes in volatility can be related to random time changes. The following is an illustrative example of the business clock, which is successfully applied in Carr et al. (2003) to control the clock of pure-jump Le vy processes.
Example 2.3 (Brownian motion, time changed via integrated Cox-Ingersoll-Ross (CIR) square-root process) Let fz t : t ! 0g be a mean-reverting positive process defined via the stochastic differential equation where z 0 , , h, are positive constants with 2h > 2 . The parameter h serves as the long run rate of time change, is the mean reversion rate, while l is the volatility of the time change. This is often called the CIR square root process (Cox et al., 1985). Taking the process fz t : t ! 0g as the instantaneous rate of time change, we define an increasing process Y t :¼ R t 0 z s ds. It is well known that the characteristic function of Y t is given by t ðyÞ :¼ E e iyY t Â Ã ¼ e Aðt;yÞz 0 Bðt; yÞ; where :¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 À 2l 2 iy q , Aðt; yÞ :¼ 2iy=ð þ cothðtÞ=2Þ and Bðt; yÞ Using the above fact, the characteristic function of where the computation is performed given Y t . For more details, we refer, for example, to Carr et al. (2003) and Schoutens (2003). &

Malliavin Calculus: Some Notations and Properties
Let us next present some basic facts on the Malliavin calculus. For details and notation, we refer to the book of Nualart, (2006). Fix T > 0, d 1 2 N, d 2 2 N, and let fW t : t ! 0g be a standard Brownian motion in R d 1 , independent of fY t : t 2 ½0; Tg, Þ 0 . Throughout this study, the time-changing process fY t : t 2 ½0; Tg is one-dimensional. Let ðf t Þ t2½0;T and ðg t Þ t!0 be the natural filtrations generated, respectively, by fY t : t 2 ½0; Tg and fW Y t : t ! 0g. Let F :¼ ðF 1 ; . . . ; F d 2 Þ 0 be a random vector in R d 2 , measurable with respect to g T :¼ ðW Y t ; t 2 ½0; TÞ. We denote by D ðkÞ Á the Malliavin derivative on the Gaussian space, conditionally on the time changing process, with respect to the k-th component of the underlying Brownian motion. Define for each k ¼ 1; . . . ; d 2 , Moreover, we denote by k the Skorohod integral over (0,Y T ] given f T with respect to the k-th component of the underlying Brownian motion, that is, for a suitable smooth stochastic process fG t : t ! 0g in R d 2 , where the multivariate Skorohod integral is taken componentwise. Both the operators D and are linked by the equality; for each k ¼ 1; . . . ; d 2 , Moreover, taking the conditional expectation, we get the duality; for each Here, we have assumed that F 2 D 1;2 and G 2 DomðÞ. We next give an integration-by-parts formula in one dimension (the multidimensional case is introduced in Section 4), which is our key tool to integrate the derivative in the following expectation for F 2 C 1 b ðR; RÞ, where X, Z and H are suitable random variables. We follow the argument of Montero and Kohatsu-Higa (2003). By the chain rule of the Malliavin derivative, we have that given f T and for t 2 ½0; Y T , where fh t : t 2 ½0; Y T g is a suitable 'smooth' stochastic process such that the terms appearing below are well defined. Integration over ½0; Y T leads to Using the duality Equation (4), we get We will often refer this relation to the (conditional) integration-by-parts formula in what follows.

Computation of Greeks
Fix T > 0 and assume that fY t : t 2 ½0; Tg is such that We first consider an one-dimensional asset dynamics driven by a time-changed Brownian motion, where the expectation is henceforth taken with respect to a risk neutral probability. With fY t : t 2 ½0; Tg in Examples 2.1, 2.2 and 2.3, the discounted stock price model (7) is clearly a martingale with respect to the risk neutral measure under the integrability condition (6) due to the Markovian property of fY t : t 2 ½0; Tg.
If the business clock fY t : t 2 ½0; Tg is deterministic, then our framework simply reduces to the most basic Black-Scholes model. To avoid such triviality, we rule out the deterministic clock a priori. Until Section 3.5, we work on the above onedimensional model and thus suppress the scripts in the Malliavin derivative D (k) and in the Skorohod integral k .

Delta
We first derive an unbiased estimator of the so-called delta, which is the sensitivity index with respect to the initial value.
where the interchange of the derivative and the expectation is justified by the dominated convergence theorem. In fact, as " # 0, and by the Taylor theorem, which is clearly uniformly bounded in ".
Next, given f T , we can apply the conditional integration-by-parts formula (5) with h t ;1 as where the last equality holds since given f T , Taking the expectation on the both sides, we get the desired result for F 2 C 2 b ðR; RÞ. It remains to remove the regularity assumption on F. To this end, let us come back to F such that E½FðS T Þ 2 is locally uniformly bounded in S 0 . We can always find a sequence fF n g n2N of continuously differentiable functions, from R to R, with compact support such that lim n"þ1 E½jF n ðS T Þ À FðS T Þj 2 ¼ 0. Hence, by the Cauchy-Schwartz inequality with E½jW Y T =Y T j 2 ¼ E½Y À1 T < þ 1, we have that for each S 0 308 R. Kawai and A. Kohatsu-Higa as n " þ1. Then, we obtain that Then, by taking limits with respect to n, we obtain that E½Fðð1 þ "=S 0 ÞS T Þ is continuous in S 0 . In a similar fashion, we can prove that E½FðS T ÞW Y T =Y T is continuous in S 0 . Finally, by taking limits with respect to n in Equation (10) and dividing by ", we obtain that E½FðS T Þ is differentiable with respect to S 0 and that the desired formula holds. The proof is complete. & Remark 3.2 The sensitivity with respect to the risk-free rate r and the drift parameter can also be derived via the identity (8), respectively, as where the interchanges of the derivative and the expectation can be justified as before. Moreover, note that the proof of Proposition 3.1 does not require the Markovian property of the time-changing process fY t : t 2 ½0; Tg. &

Vega
In our framework, we may naturally regard the parameter as the volatility. The next formula is the so-called vega, which is the sensitivity of the premium with respect to the volatility parameter.
Proof. As before, it suffices to show the result with F 2 C 2 b ðR; RÞ. We have where the interchange of the derivative and the expectation can be justified by the dominated convergence theorem in a similar manner to the proof of Proposition 3.1.

Computation of Greeks for Time-Changed Brownian Motion 309
Next, in view of the delta formula in Proposition 3.1, it remains to work on E F 0 ðS T ÞS T W Y T ½ . Given f T , we apply the conditional integration-by-parts formula (5) with h t ;1 as where the last equality holds by Equation (9). Finally, for the arguments to remove the smoothness on F, we need which concludes. &

Gamma
Next, we consider the second-order sensitivity with respect to the initial value, that is, the delta of the delta. In the derivation, we literally take 'delta' of the delta formula of Proposition 3.1.
Let F : R ! R be a measurable function such that E½FðS T Þ 2 < þ 1 is locally uniformly bounded in S 0 > 0. Assume further that E½Y À4 T < þ 1. Then, we have Proof. As before, it suffices to show the result with F 2 C 2 b ðR; RÞ. In view of Proposition 3.1, we have Then, by a similar argument to the previous proofs, we have where the interchange of the derivative and the expectation can be justified by the dominated convergence theorem. Finally, for the arguments to remove the smoothness on F, we need 310 R. Kawai and A. Kohatsu-Higa which completes the proof. &

Delta of Asian
Notice that by setting Y t ¼ t, our formulas recover the counterparts in the Black-Scholes framework. We may say, and indeed as our intuition, that in the derivation of the logarithmic derivatives, it does not matter how the clock is governed, but matters only where the Brownian motion is at the end. A delta formula of the Asian options is derived below. Interestingly, even such a path-dependent payoff yields a similar form to its counterpart in the Black-Scholes, which can be found in Montero and Kohatsu-Higa (2003).
Proposition 3.5 Fix T > 0 and define Let F : R ! R be a measurable function such that E½FðF 1 =TÞ 2 is locally uniformly bounded in S 0 > 0. Assume further that fY t : t 2 ½0; Tg is a non-decreasing Le vy process with Le vy measure defined on R þ and that there exist p . 1 and s . 1 with qðsÞ :¼ ð1 À 1=sÞ À1 such that Z z>1 e 4z ðdzÞ< þ 1; where :¼ ð4p 2 s 2 2 þ psð0 _ ÞÞ _ ð2p 2 qðsÞ 2 2 þ pqðsÞð0 _ ÞÞ and E Y À4pqðsÞ " e 8p 2 qðsÞ 2 2 Y " h i < þ1; for some " 2 ð0; TÞ: Then, we have Proof. As so far, we proceed with the proof in the case F 2 C 2 b ðR; RÞ. First, observe that where the interchange of the derivative and the expectation is justified by the dominated convergence theorem as before with E½F 2 1 < þ1. To prove this, observe that Computation of Greeks for Time-Changed Brownian Motion 311 where we have used the Cauchy-Schwarz inequality, the Fubini theorem and the infinite divisibility of Le vy processes ' Y t ¼ t' Y 1 . To continue our proof, given f T , note first that where the interchanges of the Malliavin derivative and the Lebesgue integral hold true. (See, e.g., Property P2 of Fournie et al., 1999.) By using the identity The desired formula now follows from due to the conditional integration-by-parts formula (5) with h t ;1.

Greeks for Baskets
In this subsection, we consider a basket of correlated asset prices, and its sensitivity with respect to the correlation parameter r, and its so-called cross-gamma, that is, 0 ÞÞ. We here restrict ourselves to a basket of only two assets, while our results can easily be extended to a higher dimension.
Fix r 2 ðÀ1; 1Þ and set hðrÞ : where S ð1Þ 0 and S ð2Þ 0 are positive constants. Set F :¼ a 1 F 1 þ a 2 F 2 , for some real a 1 and a 2 . Then, we have the following.
Proposition 3.6 Let F : R7 !R be a measurable function such that E½FðF Þ 2 is locally uniformly bounded in S ð1Þ 0 > 0 and S ð2Þ 0 > 0. Then, it holds that If E½Y À2 T < þ1, then it holds that Proof. In a similar manner to the previous proofs, it suffices to assume F 2 C 2 b ðR; RÞ. For the sensitivity with respect to the correlation parameter, we proceed where the interchange of the derivative and the expectation at the first equality can be justified by standard uniform integrability arguments, where the second equality can be justified by D s FðF Þ ¼ F 0 ðF Þa 2 F 2 2 hðrÞ and where the third equality holds by the conditional integration-by-parts formula (5) with h t ;1.
For the cross-gamma, we start with deriving an expression for ð@=@S ð2Þ 0 ÞE½FðF Þ, which can be done by where the last equality holds by D s FðF Þ ¼ F 0 ðF Þa 2 F 2 2 hðrÞ. With the identity D s FðF Þ ¼ F 0 ðF Þða 1 1 F 1 þ a 2 2 rF 2 Þ, we have s FðF Þds 0 @ 1 A : Using this identity, we proceed which yields the desired equality by the conditional integration-by-parts formula (5) with h t ;1. The proof is complete. (In our setting, regardless of the order of the partial derivatives and the choice of the Skorohod integrals, we arrive at the same formula.) & Remark 3.7 For all the formulas of this section, we have imposed the assumptions on the L 2 ðÞ-integrability of the payoff and on the finite negative moment E½Y Àp T < þ1, for some p > 0. As a matter of course, with a higher order integrability of the payoff and with the use of the Hö lder inequality (instead of the Cauchy-Schwarz inequality) just as in the proof of Proposition 3.5, we may decrease the order p of the negative moment to a certain extent. In particular, in the case of uniformly bounded payoffs, such as digital options (to be tested in our numerical experiments below), the order p can be taken arbitrarily small. Moreover, the order p has been increased in order to remove the smoothness of the payoff function F. Hence, we can again decrease the order p when dealing with smooth Fs. For the reader's convenience, let us here investigate the negative moments for the three candidates given in Examples 2.1, 2.2 and 2.3. For the gamma process, by using the density function (2), we can easily derive that E½Y Àp T < þ1 if and only if p 2 ð0; aTÞ. Meanwhile, it is well known that the inverse Gaussian distribution has the density Computation of Greeks for Time-Changed Brownian Motion 315 where parameters are the same as given in Example 2.2. Hence, E½Y Àp T < þ1 for every p 2 ð0; þ1Þ. Finally, for the integrated CIR square-root process, we use the identity where the functions are as given in Example 2.3. We can show that AðT; isÞ, À 2 ffiffiffiffi ffi 2s p =l; BðT; isÞ,e À ffiffiffi 2s p hT=l and thus E½e ÀsY T ,e À ffiffiffi 2s p ð2z 0 þhTÞ=l as s " þ1. This implies that the negative moment is well defined for every order. Interestingly, the above observations indicate that the inverse Gaussian process and the integrated CIR square root process cannot stay very low for long. &

Numerical Illustration
We test our formulas on a digital payoff FðxÞ ¼ 1ðx > KÞ with the normal inverse Gaussian model described in Example 2.2. In the numerical results given in Figure 1, 'FD' indicates the Monte Carlo convergence of the finite difference counterparts. For the delta, the vega, the gamma and the Asian delta, we set model parameters S 0 ¼ 100, For the baskets, we use the same inverse Gaussian process and set S 0 . We fix T ¼ 1 and r ¼ 0:05 throughout. The figures and the variance ratios evidently indicate a faster convergence of our Greeks formulae, compared to the finite difference approximations. Note that the finite difference approximations are in general biased.

Multidimensional Density Estimation
This section treats the problem of obtaining a Monte Carlo estimation of the value of the density for a multidimensional time-changed Brownian motion model. Let us begin with some additional notations on the Malliavin calculus. Fix T > 0, d 1 2 N and d 2 2 N. Let fW t : t ! 0g be a standard Brownian motion in R d 1 and define by F : ¼ ðF 1 ; . . . ; F d 2 Þ 0 a random vector in R d 2 , measurable with respect to g T . We denote by F the Malliavin matrix of F given f T , whose ðk 1 ; k 2 Þ-entry is The following is a key tool due to Malliavin and Thalmaier (our formulation follows from Bally and Caramellino, 2006). We illustrate, in Example 4.2, how the algebra actually proceeds with numerical Results.
Then, the density p F ðÁÞ of F is given by The idea of the above result is to use the fact that in the sense of generalized functions, one has that (see, for example, Evans, 1998) ÑQ d 2 ðy À xÞ ¼ x ðyÞ: Therefore, the above result can be understood as an application of the integration-byparts formula.
Example 4.2 Consider the two-dimensional model investigated in Section 3.5. Set Computation of Greeks for Time-Changed Brownian Motion 317 Hence, the Malliavin matrix F of the random vector F given f T is while its inverse almost surely exists and is given by Given f T , using the conditional integration-by-parts formula (5) with h t ;1, we get Y T À 2 hðrÞY T ; ifðl; mÞ ¼ ð2; 2Þ: Hence, we get Y T À 2 Y T hðrÞ F 2 2 Y T hðrÞ :

R. Kawai and A. Kohatsu-Higa
Hence, we get for x :¼ ðx 1 ; x 2 Þ 0 2 R 2 , Y T À 2 Y T hðrÞ F 2 2 Y T hðrÞ 2 4 3 5 : We can show that the above formula is well defined when there exists p > 1 and q > 0 such that Note that as in the case of the Greeks formulas of European type, the Markov feature of the time control process is not required in the derivation of the above estimator and we need only the sample of the marginal Y T , but not the entire sample paths fY t : t 2 ½0; Tg. We present in Figures 2 and 3 simulation results, respectively, of the variance gamma model and the normal inverse Gaussian model. & For two reasons, the above density estimation formulae outperforms the standard kernel density estimation, for example, where K is a suitable kernel and where h is the bandwidth. First, the kernel density estimation is only asymptotically unbiased as h # 0, while in reality h can never be taken zero. On the other hand, our density estimators are intrinsically unbiased. Second, in the kernel density estimation, Monte Carlo summands provide almost no contribution to the convergence of Monte Carlo simulation unless realizations are 0 ¼ 100; r ¼ 0.05, T ¼ 1, 1 ¼ 2 ¼ 0.12, 1 ¼ 2 ¼ 0.1 and ¼ 1/6. very close to x. This problem turns out to be very serious when we set the bandwidth h to be extremely small. In our formula, meanwhile, all the Monte Carlo realizations make an equal contribution.
A drawback of our formula is its infinite variance due to the singularity of the denominator F À x k k 2 . To get around this issue, we may employ an approximation technique studied in Kohatsu-Higa and Yasuda (2009) also in our framework, in exchange for the loss of the unbiasedness. In such a case, the method is similar to the kernel density estimation but the rates of bias and variance are smaller than in the classical method due to the use of the conditional integration-by-parts formula.

Concluding Remarks
In this article, we have computed formulas for the Greeks and the multidimensional density estimation for an asset price dynamics model defined with time-changed Brownian motions. Our approach is based on an application of Malliavin duality formulas on the Gaussian space conditioning on the jump component. We have noticed that the Markov property of the time-changing process is not required and that our formulas may recover the counterparts with the degenerate time clock.
On the sensitivity analysis for asset price models related to jump processes, a variety of different approaches have been proposed in the literature. In some cases, we may derive different formulas for a single object. For example, the variance gamma model can also be in the framework of Takeuchi (2008, 2009), where different formulas are obtained for delta, gamma and vega. It would be interesting to somehow investigate how those distinct formulas should be compared, for example, in terms of the Monte Carlo variance.