Computation of option greeks under hybrid stochastic volatility models via Malliavin calculus

This study introduces computation of option sensitivities (Greeks) using the Malliavin calculus under the assumption that the underlying asset and interest rate both evolve from a stochastic volatility model and a stochastic interest rate model, respectively. Therefore, it integrates the recent developments in the Malliavin calculus for the computation of Greeks: Delta, Vega, and Rho and it extends the method slightly. The main results show that Malliavin calculus allows a running Monte Carlo (MC) algorithm to present numerical implementations and to illustrate its effectiveness. The main advantage of this method is that once the algorithms are constructed, they can be used for numerous types of option, even if their payoff functions are not differentiable.


Introduction
In finance, pricing an option is the main issue of managing a trade. However, once an option is settled, its price does not remain constant. Instead, it follows a dynamic path during its survival time. Therefore, the market participants should protect themselves against the unexpected price changes by managing the variations in the option price.
The price risk and its management is always inseparably associated with the Greeks, which are derivatives of an option price with respect to its certain underlying parameters. The information that the Greeks contain is used to measure the unexpected option price changes based on specific risk factors. From this point of view, the Greeks are commonly used to construct a replicating portfolio to protect the main portfolio against some of the possible changes related to the certain risk factors. Thus, the computation of Greeks is more important than obtaining the price of an option [13], and it becomes a fundamental research area in mathematical finance. Here, it is worth to emphasize that the change in the underlying asset price is a very significant risk factor that affects the option price virtually. Hence, among other Greeks, Delta, which measures the change in the option price for a unit change in the price of underlying asset price, is an essential indicator to determine the balance of underlying asset and options which is a hedging ratio.
Over the past two decades, there has been an increasing attention focused on the accurate pricing of hybrid products that are based on a combination of underlying assets from different classes [14]. In light of the recent developments in the hybrid products, it has become difficult to ignore the computation of Greeks from which the underlying asset evolves. Thus, the study considers an extension of stochastic volatility models, namely Hybrid Stochastic Volatility (HSV) models. Within these models, the interest rate evolves from a stochastic process rather than from a constant or a deterministic function to have a more flexible multi-factor stochastic volatility (SV) model. However, the computation of Greeks under these assumptions becomes complicated. It is because these suggested models do not have an explicit distribution or because in most cases the option payoffs are not differentiable. Therefore, the likelihood method and the pathwise method are not appropriate ones in computation of Greeks for this kind of underlying assets. On the other hand, one may use the finite difference method since it relies on the approximation of a derivative as the change in a dependent variable over a small interval of the independent variable, and it can be written using a small set of difference operators. However, it is computationally expensive and problematic for discontinuous option payoff functions. In this respect, the computation of Greeks in hybrid products is not straightforward as in the Black-Scholes model (BSM).
After the pioneering studies [11,12], the integration by parts formula in the context of the Malliavin calculus has come to be considered as one of the main tools in the computation of Greeks. Since then, a considerable amount of literature has been published in this research area including some numerical applications. In most of these studies, the market dynamics are assumed to follow the BSM assumptions. However, recently, there has been an increasing interest in the computation of Greeks under the assumptions of SV models [3,9,15] and jump-diffusion type models [2,5,6].
The theory of the Malliavin calculus is attractive among both theoreticians and practitioners for three main reasons. First, it allows researchers to derive explicit weights to design an efficient MC algorithm. Second, it requires neither any differentiability condition on payoff functions (unlike the pathwise method) nor probability density functions of underlying assets (unlike the likelihood method). Third, it is more flexible than other methods in the sense that different perturbations of the Brownian motion yield different weights, which leads to an a priori interest in picking weights with small variances. Hence, it may be used for all types of option, whether they have a continuous or discontinuous payoff function.
The driving goal of this study is to introduce an expression for the computation of Greeks via Malliavin Calculus under the HSV model dynamics, and then, to use an MC algorithm for the expected values and their differentials that correspond to the Greeks. The study shows their accuracy by comparing the findings with the findings of finite difference method in all variations. The novel part of this study is to provide applicable formulas using the technical and theoretical results in the context of the Malliavin calculus. To derive general formulas for the Greeks, first, an HSV model is defined in such a way that it allows researchers to obtain generalized Greeks for all types of hybrid models using fundamental tools in the Malliavin calculus, namely integration by parts and the Bismut-Elworthy-Li formulas on the Gaussian space. The main results reveal that the Malliavin calculus gives the Greeks as a product of an option payoff function and an independent weight function called the Malliavin weight. Moreover, the numerical illustrations reveal that the computation cost for the Malliavin calculus method is less than the computation cost for the finite difference method in all variations.
This article consists of 4 sections. In Section 2, it introduces a general formula for the HSV models and the perturbed price processes. Then, it presents the main results which extend the existing theoretical Greeks formulas. In Section 3, the formulas are derived for the Heston stochastic volatility model with a stochastic interest rate, namely, the Vasicek model. Also in this section, numerical illustrations are presented for Delta, Rho, and Vega of a European call option to verify the efficiency and to compare findings with the outcomes of the finite difference method in all variations. Finally, Section 4 is the conclusion.

Computation of the generalized greeks
Consider a fixed filtered probability space (Ω, F , F t∈[0,T ] , Q) which is rich enough to accommodate a Brownian motion of dimension three for the computation convenience. Note that in this representation F t is a filtration generated by three independent standard Brownian motions W i t for i = 1, 2, 3 and Q is the risk-neutral probability measure.
Throughout the study, it is assumed that the underlying asset price evolves from the SDE system where (Z t ) i t∈[0,T ] 's are correlated Brownian motions with correlation coefficients ρ ij ∈ (−1, 1) for i, j = 1, 2, 3. These correlated Brownian motions may also be represented by a combination of three independent Brownian motions (W t ) i t∈[0,T ] such as Throughout the study, it is assumed that the correlation coefficients ρ ij are chosen in such a way that µ 3 is a real number. The solutions S t , V t , and r t represent an underlying asset, volatility and interest rate processes with initial values S 0 , V 0 , and r 0 , respectively. Here, it is assumed that σ, u, v, f , and g are continuously differentiable functions with bounded derivatives of order at least two. Moreover, σ, v, and g are assumed to be adapted and not equal to zero everywhere in the domain. Now, suppose the SDEs given in equations (1), (2), and (3) can be merged into a three-dimensional SDE system and it is represented as where, Notice that X t is a Markov process and (W t ) t∈[0,T ] is a three-dimensional standard Brownian motion. In this representation, it is convenient to assume that β and a both are at least twice continuously differentiable functions with bounded derivatives and adapted for the sake of computations. Moreover, to ensure the existence of a unique strong solution to equation (4), it is assumed that both β and a satisfy the Lipschitz and polynomial growth conditions. Under these assumptions, the study also guarantees that X t is a Markov process, the trajectories of this solution are almost surely continuously differentiable for all t up to explosion time [17].
Furthermore, in order to have an appropriate solution, it should also be guaranteed that the diffusion matrix a satisfies the uniform ellipticity condition [15], Here, it is worth to emphasize that if the conditions on µ i,j and the functions σ, v, g are satisfied, the diffusion matrix a satisfies the uniform ellipticity condition. This implies that a is a positive definite matrix, and its eigenvalues are greater than a small positive integer ǫ ∈ R. Hence, a is an invertible matrix, and its inverse is bounded. This condition also implies that for any bounded function γ : [0, T ]×R 3 → R 3 , a −1 γ is bounded, and further, (a −1 γ)(X t ) lies in the Hilbert space L 2 ([0, T ] × Ω) [11].
The first variation process of (X t ) t∈[0,T ] , which is the derivative of (X t ) t∈[0,T ] with respect to its initial value x, plays an important role in the computation of Greeks. Hence, the study reminds the definition of the first variation process.
Definition 2.1 (First variation process). Let X t be a process given by equation (4). Then, the first variation of this process is defined by where β ′ and a ′ i are the Jacobian of β and ith column vector of matrix a with respect to x, respectively. Here, 1 3×3 is the identity matrix of R 3 , and Y t = D x X t . Note that β and a are assumed to be at least twice continuously differentiable functions with bounded derivatives. Moreover, the diffusion matrix a satisfies the uniform ellipticity condition (5), and X t has continuous trajectories. Now, it is the time to introduce the following lemma for the integrity of the study.
, then X t is Malliavin differentiable and the Malliavin derivative of X t can be written as follows One can calculate the components of the first variation process as in the following proposition. Proposition 2.3. Let X t and the first variation process Y t for t ∈ [0, T ] be defined by equations (4) and (6), respectively. Then, Y ij t = 0 and Y 23 Furthermore, Y 12 t and Y 13 t satisfy Proof. In order to find the entries of the system of SDEs which is satisfied by Y , one has first to calculate the Jacobian of β and the ith column vector, a i , of the diffusion matrix a for i = 1, 2, 3. These are By inserting these equations into equation (6), one obtains the first variation process as with an initial value of the identity matrix in R 3 . Then, one can easily deduce that Y ij t = 0 and Y 23 t = 0 a.s. if i > j for i, j = 1, 2, 3 by applying the Itô calculus. Moreover, for t ∈ [0, T ] with the initial values Y 11 0 = Y 22 0 = Y 33 0 = 1. Remark 2.4. As an immediate conclusion of equation (8), Y 11 t may be written in terms of the asset price, i.e.
For the computation purposes, it is convenient to consider a continuous time trading economy with a finite time horizon t ∈ [0, T ]. Suppose that the uncertainty in this economy is idealized by a fixed filtered probability space (Ω, F , F t∈[0,T ] , Q) and the information evolves according to the filtration Furthermore, consider an option in this market, which has a square integrable payoff function denoted by Φ = Φ(X t1 , . . . , X tn ), where Φ is a payoff function that is continuously differentiable with bounded derivatives. In other words, the option is Now, following studies [5,9], the option price, p, at time t = 0 with a maturity T < ∞ is traditionally calculated as the expected value of the discounted payoff at maturity conditionally to the present information, which is described by the σ-algebra F 0 . Now, it is possible to define the price process in the following definition.
Definition 2.5. Price, p(x), of an option with an underlying X and a payoff function Φ, is defined by where, 0 = t 1 , . . . , t n = T is a partition of the finite time horizon [0, T ], and p(x) denotes today's fair price of the option. Namely, in this study E[.] is the expected value under the risk-neutral probability.
Note that the objective of the computation of Greeks is to differentiate the price, p, of an option with respect to model parameters.
Remark 2.6. In the computation of Greeks, it is fundamentally assumed that the option payoff function, Φ, is continuously differentiable with bounded derivatives. However, they are generally not smooth in genuine markets. In this case, if the law of the underlying asset is absolutely continuous, and the option payoff function, Φ, is Lipschitz, it is possible to derive explicit Malliavin weights for the Greeks by Proposition A.4.
In the computation of Greeks via the Malliavin calculus, a weight function, which is independent of the option payoff function, is obtained. To obtain a valid computation result, one has to guarantee that the Malliavin weights do not degenerate with probability one. Hence, [11] demonstrates the set of square integrable functions Proposition 2.7. Suppose that the functions β and a in equation (4) are continuously differentiable with bounded derivatives, and the diffusion matrix a satisfies the uniform ellipticity condition (5). Moreover, the option payoff function, Φ, is square integrable, and continuously differentiable with bounded derivatives. Now, consider the price process given in Definition 2.5 with a maturity T < ∞, then where ∇, and π denote the gradient and Malliavin weight, respectively. Here, where α ∈ Γ n , and Y t is the first variation process.

Computation of Delta
Delta of an option is a measure of variations in its price with respect to initial underlying asset price, and it determines a hedging ratio. It can be computed by using Proposition 2.8.
Proposition 2.8. Suppose the functions β and a are both as in equation (4). Moreover, they are continuously differentiable functions with bounded derivatives and the diffusion matrix a satisfies the uniform ellipticity condition (5). Now, consider an option with payoff Φ, which is a continuously differentiable function with bounded derivatives. Then, Delta of the option with the price function (11) is where ∆ MW is the Maliavin weight of Delta and it is Proof. Using the inverse of the diffusion matrix a An immediate conclusion of Proposition 2.7 with the special choice of α(t) = 1 T is given by Here one may choose α(t) = 1 T since European options are priced at maturity and therefore t i = T . Then, with Remark 2.4, the final result may be obtained easily.
As it is seen in the above formula, the gradient of the option price is denoted by ∇p(x) = ( ∂p ∂S0 , ∂p ∂V0 , ∂p ∂r0 ) ⊤ , where ⊤ denotes the transpose. The first row of the solution corresponds to the option's Delta. The remaining two rows correspond to the changes in the price with respect to the initial volatility and initial interest rate, respectively. Hence, as a consequence of this result, one can present the following two remarks. Remark 2.9 (Vega Vt ). The sensitivity of an option to its initial volatility is Remark 2.10 (Rho rt ). The sensitivity of an option to the initial interest rate is

Computation of Rho
The computation of Rho is not as straightforward as is the computation of Delta since the interest rate is neither constant nor deterministic. Hence, instead of directly differentiating the option price with respect to the interest rate, one may consider adding a perturbation term ǫ to the drift term and then try to observe the effect of the perturbation on the option. Here, it is necessary to clarify what exactly is meant by a perturbed process X ǫ t to observe the change in the price with respect to change in the drift term.
As in [11,5], the study introduces the perturbed process (X ǫ t ) t∈[0,T ] as follows: where ǫ is a small scalar and γ : [0, T ] × R 3 → R 3 is a bounded function. Furthermore, β and a satisfy the regularity conditions that are discussed above.
To interpret the impact of a structural change in the drift and the price, one should perturb the price process as in the following definition.
Definition 2.11. Suppose X ǫ t is the solution of the SDE system given in (15) for t ∈ [0, T ] and Φ is a continuously differentiable function at least order two with bounded derivatives. Then, the perturbed price process p ǫ (x) is given by Now it is convenient to present the following proposition to show the sensitivity of the option to the parameter ǫ in the point ǫ = 0. Proposition 2.12. Suppose that β, a are continuously differentiable functions with bounded derivatives, and a satisfies the uniform ellipticity condition (5). Then, for any square integrable and continuously differentiable function with bounded derivatives Φ, ǫ −→ p ǫ (x) is differentiable at any x ∈ R 3 and Proof. See the proof in [11].
Proposition 2.13. Suppose β and a are continuously differentiable functions with bounded derivatives. Moreover, a satisfies the uniform ellipticity condition (5), and the square integrable option payoff function Φ is a continuously differentiable function with bounded derivatives. Then, Rho of the option is where Proof. Rho measures the effect of a change in the interest rate on the option price. In Proposition 2.12, there are mainly three sources of perturbation: drift terms of risky asset, volatility, and interest rate processes. The function γ(X t ) can be chosen as any combination of these three sources. Since the study is investigating the effect of the interest rate on the option price, it should perturb the original drift with γ(X t ) = (S t , 0, 0) ⊤ . In this case, γ is a bounded function since t ∈ [0, T ] is a continuous time trading economy with a finite time horizon. Note the fact that the dynamics are given under a risk-neutral probability measure and perturbation. Hence, the discount process becomes e − T 0 (rt+ǫ)dt . First one should find Then, by inserting the equation above into the expectation term, one obtains It is also possible to compute the effect of other parameters on the option by special choices of γ. For instance, as an immediate result of Proposition 2.13, one can present the following two remarks.
Remark 2.14. Suppose that the assumptions in Proposition 2.13 hold and the stock price evolves from the Heston model with a stochastic interest rate. Then, if γ(X t ) = (0, κ, 0) ⊤ , one obtains the sensitivity of the option price with respect to κ. Remark 2.15. Suppose that the assumptions in Proposition 2.13 hold and the interest rate is assumed to follow the Vasicek interest rate model. Then, if γ(X t ) = (0, 0, a) ⊤ , one obtains the effect of "speed of reversion" parameter on the option price.

Computation of Vega
For the computation of Vega, one need to define a new perturbed process as in the computation of Rho, but in this case, the perturbation will occur in the diffusion term. However, in the end, it is necessary to calculate the Skorohod integral. Hence, the result of Proposition A.6 will be used.
The perturbation approach in this section is based on the approach used in [5,11]. First, consider the perturbed asset price process where ǫ is a small scalar, γ is a 3 × 3 matrix valued continuously differentiable function with bounded derivatives. Furthermore, β and (a+ǫγ) satisfy the aforementioned regularity conditions. Here it is necessary to introduce a variation process with respect to ǫ, which is the derivative of X ǫ t with respect to the parameter ǫ, To avoid degeneracy, the setΓ n of square integrable functions in R, (t)dt = 1, ∀i = 1, . . . , n , is defined in [11]. The following Proposition tells how sensitive the price of an option on the perturbed process is to ǫ in the point ǫ = 0.
Proposition 2.16. Suppose that a satisfies the uniform ellipticity condition (5) and . . , n, there exists a −1 (X)Y B ∈ Dom(δ). Then, for any square integrable option payoff function, Φ, with continuously differentiable and bounded derivatives, holds. Here,B for t 0 = 0 andα ∈Γ n . Moreover, if B is Malliavin differentiable, the Skorohod integral is calculated according to Remark A.10 and it is Proof. The proof can be found in [5].
Proposition 2.17. Consider the three-dimensional SDE (4) and its perturbed process (18). Assume that β and a are continuously differentiable functions with bounded derivatives and, moreover, a satisfies the uniform ellipticity condition (5). Then, Vega of an option with a square integrable payoff function Φ, which is continuously differentiable with bounded derivatives, is where Vega MW is the Maliavin weight of Vega, and Proof. First obtain a perturbed process by perturbing the original diffusion matrix with γ, where it is chosen as Then, using the fact that V t and r t do not depend on S t , one can deduce that the variation process Z ǫ t has vanishing components; so in Z ǫ=0 t the components Z 2 t and Z 3 t are almost surely zero. Then, as in Proposition 2.16, define the vector Here, Then, substituting these two equation into B ti , one obtains On the other hand, from equation (19) it is known that Z ti satisfies the following dynamics for t ∈ [0, T ]. From this setting, one can write With the Itô formula, one can easily find the solution as Thus, using equation (21) and Remark 2.4, one has for t i ∈ [0, T ]. According to Proposition 2.16, the Skorohod integral δ(a −1 (X)YB) remains to be calculated. Here, B · is Malliavin differentiable and its Malliavin derivative is Then, one obtains the trace .

As a result, by choosingα
dt .
Remark 2.18. If the option payoff depends on only the maturity T , then dt .

Numerical illustration
This section is devoted to numerical illustrations of the Greeks of a European call option with a strike price K and an option payoff function Φ = max{S T − K, 0}, where S T is the price of the underlying asset at maturity T < ∞. It is also worth to emphasize that a European call option has a Lipschitz payoff function and it belongs to the space of locally integrable functions denoted by L 2 . Further, the space C ∞ c of infinitely differentiable functions having a bounded compact support, where c is an arbitrary compact subset of R 3 , is dense in L 2 . Hence, there exists a sequence of functions Φ n ⊂ C ∞ c that converges to the main payoff function Φ ∈ L 2 , see A.4. Therefore, one can apply the formulas introduced in the previous section to European options.
The formulas introduced in Section 2 are for a general case since the functions in the SDEs (1)-(3) are given with closed forms. Therefore, it is possible to find the Greeks for all stochastic volatility models through the special choice of functions that are introduced in SDEs (1)- (3). In this study, these functions are chosen according to the well-known Heston stochastic volatility model with a stochastic interest rate, namely the Vasicek model for the simulation purposes. Under these special choices, the SDEs (1)-(3) become where the initial values are S 0 , V 0 and r 0 , respectively. Here, it is assumed that the coefficients κ, θ, σ, a, b, and k are all positive numbers. Under this setting, the functions β and a in (4) are continuously differentiable, and satisfy the Lipschitz condition. Moreover, a satisfies the uniform ellipticity condition (5). These assumptions are enough to compute the Greeks in the BSM framework. However, one needs more assumptions in the Heston stochastic volatility model framework because the square root function is not differentiable and not globally Lipschitz. The Novikov condition for the Heston stochastic volatility model, κθ ≥ σ 2 , guarantees that the volatility process is always positive. Hence, it is assumed that the Novikov condition is satisfied, and the initial volatility V 0 is positive. Moreover, in the studies of [1,10], it is proved that the Heston stochastic volatility model is Malliavin differentiable under the Novikov condition. In the rest of the paper, it is assumed that the Heston stochastic volatility model satisfies the Novikov condition and it is Malliavin differentiable.
Using the above given values, the study presents simulations for the Greeks Delta, Rho and Vega of a European call option after conducting successfully the computations in the previous section. Figures 1, 2 and 3 allow comparing the finite difference method in all variations and the Malliavin calculus on sample sizes. The computed Delta, Rho and Vega values of both methods are very stable and quite good, even for a low number of MC simulations. Furthermore, if the number of simulation increases, the Greeks values become more stable for each method. Therefore, one should increase the number of simulations for both methods to have a more accurate value for the Greeks.

Conclusion
In this study, a very general model skeleton is considered in equations (1), (2), and (3) for the computation purposes, and the formulas with detailed proofs for the well-  known Greeks under the assumptions of stochastic hybrid volatility models are derived. As in many other studies, the Greeks are obtained as an expectation of product of two terms: a payoff function and a weight called the Malliavin weight which is independent of the payoff function. This result indicates that the efficiency of the Malliavin calculus in the computation of the Greeks does not depend on the type of the payoff function. In the computations, it is assumed that the payoff function is continuously differentiable. However, in the case of mathematical finance, payoff functions are not globally differentiable, so it is better to mention that explicit ex-pressions given in this paper are easily extended to payoff functions which are not continuously differentiable. Hence, once the formulas are obtained for one option, one can use the same formula for all kinds of options by changing only the payoff function. Moreover, by substituting the necessary functions into the Greeks formulas, one can obtain the Greeks for all kinds of models, and these formulas can be easily adapted to the special needs of financial engineers working in practice on the computation of Greeks. As an application of the results, the paper also examines a particular case of the Heston stochastic volatility model by assuming the interest rate evolving from the Vasicek model. In order to compare the results, these Greeks are computed with the finite difference method in all variations. It is observed that the formulas which are obtained by using the Malliavin calculus yield results that require a fewer number of simulations than in the finite difference method for Vega and Rho. Moreover, despite the easy implementation of the finite difference method, the duration of computation is higher than in the Malliavin calculus. Since traders need the Greeks for hedging purposes, they have to be computed as fast as possible. It is, therefore, using the Mallaivin calculus is superior to the finite difference method in the computation of Greeks because once the formulas are obtained, they can be used for all types of option and the duration of computation is shorter than in the finite difference method.
A A brief review on Malliavin calculus where ∂f ∂xi is the partial derivative of f with respect to its ith variable. Proposition A.2 (Chain Rule). (Proposition 1.2.3 in [16]) Suppose that F = (F 1 , . . . , F n ) is a random vector whose components belong to the closure of S, D 1,2 , and the function ϕ : R n → R is a continuously differentiable function with bounded partial derivatives. Then ϕ(F ) ∈ D 1,2 and almost surely for t ∈ [0, T ].
Proof. If the function ϕ is smooth the proof can be obtained by the chain rule in classical analysis. Otherwise, the function has to be mollified. In order to mollify ϕ, one can use ρ ǫ (x) = ǫ n ρ(ǫx), where ρ(x) = ce 1 x 2 −1 and c is a chosen coefficient that makes the integral R n ρ(x)dx = 1, to obtain a smooth approximation ϕ * ρ ǫ . Considering the smooth approximations F n of F one obtains ϕ * ρ ǫ (F n ) → ϕ(F ) for min ǫ, n → ∞ in the space L 2 . Then by closedness of the derivative operator D Lemma A.3. Suppose that the sequence F n ∈ D 1,2 converges to F in the space L 2 (Ω, F , P) satisfying sup n E[ DF 2 H ] < ∞. Then, F ∈ D 1,2 and DF n weakly converges to DF in L 2 (Ω × [0, T ]).
A generalization of Proposition A.2 to the functions even not necessarily differentiable is given below.