Nonparametric Malliavin–Monte Carlo Computation of Hedging Greeks

: We propose a way to compute the hedging Delta using the Malliavin weight method. Our approach, which we name the λ -method, generally outperforms the standard Monte Carlo ﬁnite difference method, especially for discontinuous payoffs. Furthermore, our approach is nonparametric, as we only assume a general local volatility model and we substitute the volatility and the other processes involved in the Greek formula with quantities that can be nonparametrically estimated from a given time series of observed prices.


Introduction
In finance, the numerical computation of the option price sensitivities, named Greeks, has attracted much attention in the last few years, as the knowledge of the Greeks and their implementation as risk management tools have paramount importance when trading financial derivatives. However, their numerical computation using the Monte Carlo finite difference method is particularly inefficient in the case when the option payoff is represented by a discontinuous function of the underlying asset, as it is often the case.
In their seminal paper, Fournié et al. (1999) show that the stochastic calculus of variations, also referred as the Malliavin calculus (see Malliavin 1997), is the correct tool for computing hedging Greeks as it allows us to consider the effects of perturbations of the model (e.g., perturbation on the initial data, volatility process, ...) on functionals of random variables expressing the payoff. Fournié et al. (1999Fournié et al. ( ), (2001 propose to exploit the integration by parts formula based on the Malliavin calculus to avoid computing the derivatives of possibly non differentiable payoff functionals as a strategy to get efficient Monte Carlo methods for computing the Greeks. The integration by parts formula can be performed in different ways, producing alternative weights, which are numerically more or less efficient. Even if the different expressions are equivalent from a mathematical point of view, when performing numerical computations they slightly differ, even if the same random series is used. The problem of the optimality of the weight-i.e., providing the minimal variance-is addressed in Kohatsu-Higa and Monteiro (2001) and Benhamou (2003).
Furthermore, the practical implementation of the Malliavin-Monte Carlo method requires the knowledge of the precise form of the diffusion function of the underlying process-i.e., the latent volatility process. This is a crucial point both for pricing and hedging, as demonstrated by the huge amount of literature on determining the correct specification of the underlying model-see Aït-Sahalia (1996). Moreover, the specification of the functional form of the diffusion function is even harder to determine. In this respect, the use of a nonparametric estimator of the volatility can be substantial in terms of avoiding misspecification of the hedging Greeks. Hedging errors due to model misspecification have been investigated from different points of view-see Bouleau (2003); Hayashi and Mykland (2005) and the references therein.
In this paper, we address both issues. We propose a way to compute the weight in the Malliavin weight method that yields a faster convergence rate. Furthermore, our approach is nonparametric, as we only assume a general local volatility model and we substitute the volatility and the other processes involved in the Greek formula with quantities that can be nonparametrically estimated given the time series of observed prices. This result can be achieved by expressing the weight using the rescaled variation which has been introduced in Barucci et al. (2003); Malliavin and Mancino (2002b).
An important advantage of expressing the Malliavin weight in terms of the rescaled variation relies in the fact that this function satisfies an ordinary differential equation instead of a stochastic differential equation. This result has an impact on the numerical efficiency of the Malliavin-Monte Carlo method. In fact, as long as the method is applied only to the case in which the underlying process is a Black-Scholes model, as in Elie et al. (2007);Fournié et al. (1999), the numerical computation of the weight involves only one Itô integral. Nevertheless, this is a very special case: as the underlying SDE is linear, the linearized SDE for the variation process and for the log-normal process are the same. The situation is drastically different if the underlying model is non-linear, as for the local volatility models as well as for stochastic volatility models. Under these models, the expression of the weight involves more than one stochastic integral. The simplest numerical approximation to these stochastic integrals is given by the Euler-Maruyama method. It is well known that this method has strong order of convergence 0.5 and weak order of convergence 1 for Itô stochastic integrals, while it is of order γ = 1 for Riemann integrals. Therefore, when we consider an underlying model which is not log-normal, Fournié et al. (1999) account for numerically computing a double stochastic integral, while our approach requires us only to compute one stochastic integral.
In this paper, we assume an underlying local volatility model,; namely, we model volatility as a level-dependent quantity (i.e., a deterministic function of the asset price). There are different motivations for this kind of model. For instance, in a single asset, market volatility depends on the asset price and, therefore, the market is complete. Moreover, this way to model asset price volatility is well suited to capture the relationship between volatility and asset price-returns (leverage effect). The simplest way to model the negative relation between asset-price returns and volatility is to assume a constant elasticity of variance (CEV) model-see Cox and Ross (1976). Furthermore, level dependent volatility has also been employed to reproduce the implied volatility smile, see Derman and Kani (1994); Dupire (1994); Hobson and Rogers (1998).
The paper is organized as follows. In Section 2, we describe the Malliavin weight approach to Greek computation and its nonparametric application; moreover, the Black-Scholes model and the CEV model are briefly analyzed. In Section 3, we show the computational advantage of computing the weight function by a single stochastic integral instead of a double stochastic integral; in addition, we show the computational efficiency of the method when compared to the standard finite difference Monte Carlo method. Finally, in Section 4 we illustrate an empirical nonparametric implementation of our approach for the EUR/USD exchange rate options on the period from 29 February 2016 to 3 June 2016. Section 5 concludes. The Appendix A resumes the mathematical computation of the rescaled variation and its model-free expression.

Greeks Computation via Malliavin Calculus and the Rescaled Variation
In this section we recall the main ingredients for the computation of the Greeks using the Malliavin calculus, referring to Fournié et al. (1999) and Malliavin and Thalmaier (2006) for details. The underlying asset price is described by a diffusion process (S W (t)) 0≤t≤T whose dynamic is described under the risk neutral probability by the following SDE 1 where a(·) is a deterministic function of the asset price, r the constant risk-free interest rate and W a Brownian motion on a filtered probability space. Consider a contingent claim with a payoff function of the form f (S W (T)) for some function f . The no-arbitrage price of the contingent claim is given by where E is the expectation under the risk neutral probability and the notation S W (T, s 0 ) denotes the solution of (1) at time t = T which has initial condition s 0 . From a mathematical point of view, the Greeks formulae are given by the directional variations of u(s 0 ) in response to a perturbation of the initial datum (in the case of Delta), of the volatility coefficients (in the case of Vega) or of different parameters. In the present paper we focus on the computation of Delta, the first derivative of u(s 0 ) with respect to the initial datum. Consider the variation of the initial condition s 0 → s 0 + ε, for which the Delta is given by If the expression (3) does not have an explicit-closed formula (because in most cases the distribution of S W (t) is not known-see Remark 2), we may think of using Monte Carlo simulations in order to compute the expectation. Moreover, as the function f is usually not smooth, we can resort to the computation of difference quotients (i.e., finite differences) to perform the derivative in (3). Fournié et al. (1999) showed that it is possible to accelerate the rate of convergence of the Monte Carlo method and to overcome the lack of classical differentiability of the payoff function f , if we are able to express the differentiation involved in the Greeks formulae-e.g., (3), as the expectation of the payoff function multiplied by a weight π, which does not depend on the payoff function itself. The authors study numerical applications of the method in the case where the volatility is constant, showing its efficiency. We sketch the idea behind this approach.
In the context of stochastic calculus of variation (see Malliavin (1997)), a perturbation of the initial condition propagates along the trajectory W as where h(t) := t 0 z(s)ds and the process z(s) is the rescaled variation defined as the ratio between the variation process ζ(t) associated to (1) and the diffusion coefficient of (1), namely 1 We emphasize the dependence on the Brownian path, using the subindex W.
Remember that the variation process follows an SDE, which is the linearized equation of the SDE (1) driving the price 2 Then, we have where D h is the Malliavin derivative along the direction h and the last equality is obtained by means of the integration by parts formula (see Malliavin (1997)). The function w(t) can be arbitrarily chosen but it must satisfy T 0 w(t)dt = 1. For simplicity we will consider w(t) = 1/T. Therefore, the weight π is equal to (1/T) T 0 z(t) dW(t) and it does not depend on the payoff function. Furthermore, formula (7) does not require us to differentiate the payoff function. 3 The computation of the weight in (7) is far from being trivial. In fact, it is worth noting that apriori the computation of the weight π involves a double stochastic integral, as the rescaled variation z(t) is the solution to an SDE. However, in Malliavin and Mancino (2002b); Malliavin and Thalmaier (2006) it is proven that the rescaled variation actually evolves in time according to a linear ordinary differential equation. More precisely, the rescaled variation is differentiable with respect to t and it holds that where The general case of this result in stated in Appendix A as Theorem A1. In practice, the Greeks can be calculated in a parametric framework as follows: first determine a parametric model, as a general diffusion process on the underlying process, that can produce the best fit, and estimate the parameters in the model; then, calculate the Greeks using the selected model. We take a nonparametric approach instead. The function λ contains terms, such as the volatility function and its derivatives, which are latent and should be estimated from the data. Following Barucci et al. (2003); Malliavin and Mancino (2002b), it can be expressed in terms of three factors which can be estimated nonparametrically given high frequency asset return data using the Fourier methodology introduced by Malliavin andMancino (2002a, 2009). Therefore, the rescaled variation is expressed by means of quantities which can be empirically estimated using the price process observations.More precisely, the function λ(t) in (9) can be expressed by means of terms which are iterated cross volatilities. The result is the following: Denoting by , the quadratic (co-)variation operation, define the following cross-volatilities: The prime stands here for the first derivative with respect to the level S W (t).

3
It should be noticed that the last integral in (7) is a Skohorod integral, for the duality formula with the Malliavin derivative operator D h . However, if the integrand is adapted to the filtration, then the Skohorod integral and the Itô integral coincide.
then the function λ t in (9) has the following expression The general case of this result is stated in Appendix A as Theorem A2. On their turn, these cross volatilities are stochastic functions that can be estimated using the Fourier estimation method-see Mancino et al. (2017). In particular, consistent estimators of the stochastic functions A t , B t and C t have been studied in Mancino and Sanfelici (2020).
The next subsection illustrates the methodology by means of two examples, the first one is the Black-Scholes model, the second is a level dependent volatility model, the Constant Elasticity of Variance (CEV) model.

Two Examples
We illustrate the methodology first on the Blac-k-Scholes model, then on the CEV model in the case of the European-type options.

The Black-Scholes Model
Under the risk neutral probability, let S t follow a Black-Scholes model Note that this is a linear SDE, thus the SDE for ζ(t) coincides with that of S(t), with different initial data dζ(t) = rζ(t)dt + σζ(t)dW(t), ζ(0) = 1 and s 0 ζ(t) = S(t). Using (7), as in Fournié et al. (1999) we obtain the Delta: Remark 1. Note that, for the Black-Scholes model s 0 , constant-which implies that λ = 0. This can also be derived by using the expression (9).

The CEV Model
Consider now the CEV model for S(t), that is, under the risk neutral probability, The associated SDE of the variation process is The non linearity of the diffusion coefficient-i.e., σS δ -drastically changes the relation between S(t) and ζ(t). The solution of the linear SDE (13) is The computation of the Delta using (7), as in Fournié et al. (1999), gives: which entails the computation of two stochastic integrals.
On the other side, by using the λ-method, through the formula (8) we get: which entails the computation of a single stochastic integral. The computation of λ for the CEV model can be easily obtained from (9) and is equal to Remark 2. Note that even for a plain vanilla call option, there is no explicit expression for the non arbitrage price (and thus the Delta) assuming a CEV model. For the pricing and hedging of options, the available formulae are expressed by a series of incomplete Gamma functions-see, for instance, Davydov and Linetsky (2001) .

Simulation Study
The aim of this section is twofold: first we show the computational advantage of computing the weight function π by a single stochastic integral as in (15) versus a double stochastic integral as in (14); second, we show the computational efficiency of the method when compared to the standard finite difference Monte Carlo method.

Computation of the Weight
As a case study, we consider the CEV model. We compare two alternative ways to compute the Malliavin weight, namely by the processes and Equation (18) exploits the knowledge of the price-volatility feedback rate λ, and requires the computation of a Riemann integral, while the term (17) involves the computation of a stochastic integral.
The simplest numerical approximation to these stochastic processes is given by the Euler-Maruyama method. It is well known that this method has strong order of convergence 0.5 and weak order of convergence 1. Although the Delta values in (14) and (15) are expressed as expected values, it is important to quantify the ability of the method to compute strong (i.e., pathwise) solution on average, because this reflects on the weak convergence as well.
We discretize the interval [0, 1] into M = 400 subintervals of length dt = 1/M with nodes t j = j dt. We generate a discretized Brownian path which is used to compute a reference ("true") trajectory for the processes φ and ψ. Then, we consider other decomposition of the interval [0, 1] into L subintervals in such a way that the stepsize ∆t = 1/L is an integer multiple R ≥ 1 of the increment dt used for computing the Brownian path. This ensures that the set of nodes {t j } contains the points {τ j } at which the convergence of the Euler-Maruyama method is tested. Our choice is to double the discretization step each time by taking R = 2 p , with p = 1, 2, 3, 4. The Brownian motion increments on the reference path are aggregated in order to build the increments W(τ j ) − W(τ j−1 ) on the sparse grids that are then used to implement the Euler method. The parameter values are K = 100, σ = 0.2, δ = 0.8, T = 1, s 0 = 100 and r = 0.0.
The upper panels of Figure 1 show the effects of the discretization on the asset price S(t) and lambda paths λ(t). The blue lines are the reference trajectories on the finer grid, while the red dotted lines are the discrete approximations obtained for p = 4. In the lower panels of Figure 1, we plot the corresponding approximated trajectories of the functions φ and ψ over the discrete grids {τ (p) j } for p = 1, 2, 3, 4. We notice that the function φ(t) is much affected by the discretization over sparse grids, while ψ(t) is not. Nevertheless, the approximated functions match the reference ones computed on the fine grid more closely as ∆t decreases. Therefore, convergence seems to take place. Finally, we remark that the stochastic function ψ(t) seems to be much regular. Indeed, this is in complete agreement with Equation (8), which proves the differentiability with respect to time of the rescaled variation, while the differentiability of function φ(t) is not so visually evident. In order to test the strong convergence of the two different methods, we compute e strong ∆t = E[|ẑ j − z(τ j )|], where z is the reference rescaled variation computed by (17) or by (18) on the fine grid andẑ j its discrete approximation. A method is said to have strong order of convergence equal to γ if there exists a constant C such that for any fixed τ j and ∆t sufficiently small. It is well known that the Euler method of integration is of order γ = 1/2 for Itô stochastic integrals, while it is of order γ = 1 for ordinary differential equations. This marks a departure between the stochastic and the deterministic settings and can support and motivate our approach of computing the rescaled variation by the Riemann integral (18) instead of the Itô integral (17). In our numerical tests, we focus on the error at the endpoint τ = 1 and we test whether the bound e strong ∆t ≤ C∆t γ holds for some γ. We compute N = 1000 different Brownian paths over [0, 1] with M = 400 and then we apply the Euler-Maruyama method with different stepsizes ∆t = 2 p dt, with p = 1, 2, 3, 4. If the condition (19) holds with approximate equality, then we get log e strong ∆t ≈ log C + γ log ∆t, namely plotting the approximated errors e strong ∆t against ∆t on a log-log scale would produce a line with slope γ.
In Figure 2, the blue asterisk line plots our simulated results. For reference, a dashed line of slope γ is plotted as well. In both cases the error e strong ∆t reduces as ∆t becomes smaller but in the case of the lambda approach the slope γ is 1, suggesting a strong order of convergence of 1, while in the other case the convergence rate is γ = 1/2. We can analyze this further by testing a power law relation e strong ∆t = C∆t γ or, equivalently, log e strong ∆t = log C + γ log ∆t. A least squares fit for log C and γ is computed producing the values 1.2945 and 0.4849 for γ in the two cases, respectively. Now we come to examine the convergence of the approximation to the weight π computed in the two alternative forms and As before, we compute the mean of the error E|Ẑ ∆t − Z| of the approximated integralẐ ∆t as ∆t becomes small. In the upper panels of Figure 3 we can see that again the error reduces with ∆t but, while Φ exhibits a strong order of convergence equal to 1/2, the approximation of Ψ still retains a strong order of convergence equal to 1. More precisely, the least squares fit for γ produces the values 1.2025 and 0.6196 for the strong convergence regression. This different behavior of the two methods has important consequences even on the weak order of convergence of the integral. This is a less demanding alternative to strong convergence which measures the rate of decay of the error of the means |E[Ẑ ∆t ] − E[Z]|. This error is plotted against ∆t on a log-log scale in the lower panels of Figure 3. As expected, the approximation of Ψ exhibits a weak order of convergence γ = 1 typical of the Euler-Maruyama method, while in the case of Φ the convergence rate is not much clearer. The least squares fit for γ produces the values 1.0802 and 0.3307 for the weak convergence regression with the two different methods, respectively. However, the least squares residual for Φ is 1.4119 and does not support a clear order of weak convergence for a double stochastic integral approach.

Efficiency of the Method
In order to asses the efficiency of our method, we compute the ∆ in the case of a Black-Scholes and of a CEV underlying diffusion for the following payoffs: using the λ-method, the method proposed by Fournié et al. (1999) and the Monte Carlo finite difference method. The results are displayed in Tables 1-4. "Analytic" denotes the closed form solution.
First, let us consider the Black-Sholes framework. In this case, the λ-method and the method of Fournié et al. (1999) are the same. However we implement the formula (11) for the ∆ in two different ways, either by computing the expectation in the right hand side as in Fournié et al. (1999) or by simulating trajectories of the Brownian motion to compute the integral in the left hand side (λ-method).  Fournié et al. (1999) as well as the Monte Carlo finite difference method for small number of iterations N, i.e., the λ-method exhibits a faster convergence speed. The perturbation of the initial condition in the Monte Carlo approach is set to h = 0.01. The Monte Carlo finite difference method performs quite poorly for Digital options, while the λ-method gives very good results.
When increasing the number of iterations, the results provided by the three methods are more comparable and almost equivalent as shown in Tables 2 and 3.  (15), where the function λ is given by (16), while the Fournié et al. (1999) method is based on the computation of the expectation in (14). The Monte Carlo finite difference method is implemented with common random numbers. The initial perturbation is set to h = 0.01 and in the case of digital options a second order different quotient is used to approximate the delta. The Analytic delta value is computed in terms of the noncentral chi-square distribution-see Hull (2018).
As can be seen from Table 4, the results provided by the λ-method and the method of Fournié et al. (1999) are almost undistinguishable, while the Monte Carlo finite difference method performs slightly better in the European case but is completely unreliable for binary options with a 20% relative error on the Delta value.  Figure 4 shows the performance of the λ-method versus the Monte Carlo finite difference for increasing values of the number of simulations N used to compute the expectation. In the case of the European option the Monte Carlo method provides slightly less biased results. This is a well understood phenomenon and a possible way to fix this problem is to localize the integration by parts around the strike price K, as suggested by Fournié et al. (1999). In contrast, for the Digital option case the λ-method is very efficient, while the Monte Carlo finite difference methods gives very poor results.

Greek Estimation from Observed Underlying Prices. An Empirical Case Study
In order to illustrate the model-free nature of our approach in the context of the Malliavin weight method, we consider a case study based on 5-minute observations (open prices) of the EUR/USD exchange rate on the period from 29 February 2016 to 3 June 2016. We propose a nonparametric estimation procedure for computing the Delta of exchange rate options. Prices of derivative securities depend crucially on the form of the volatility of the underlying process, therefore we work in the framework of local volatility models but leave the volatility function unrestricted and estimate it nonparametrically from the observed prices.
Our sample contains 20,136 observations. Table 5 describes the main features of our data set. Intra-day returns may be contaminated by transaction costs, bid-and-ask bounce effects, etc., leading to biases in the variance measures. A preliminary analysis through the autocorrelation function plot of the 5-min returns is conducted in Figure 5 to detect the presence of market microstructure noise. The degree of first order autocorrelation of log-returns seems very weak, meaning that microstructure effects are not relevant in this sample. However, the autocorrelations decay slowly. To test for stationarity, we performed an augmented Dickey-Fuller nonstationarity test. The nonstationarity hypothesis is rejected at the 90% level. The test is known to have low power, so even a slight rejection means that stationarity of the series is likely.  We assume the risk neutral dynamics for the spot exchange rates S is Then, the logarithm X t = log S t follows a time-homogeneous Itô diffusion process represented by the following stochastic differential equation with initial condition X 0 = x 0 , where W t is a standard real Brownian motion and the real functioñ σ(x) = σ(exp(x)) is such that a single solution X t of the stochastic differential Equation (23) exists. We remark that we are making no assumptions on the functional form ofσ(·). Our specific problem is to estimate the diffusion termσ(x) when we observe a discrete realization of the process X t via n + 1 observationsX 0 , ...,X n at times t 0 = 0 < t 1 < ... < t n = T in the interval [0, T] and then to price derivatives in the context of the Malliavin weight method. FIRST STEP: We estimate the diffusion termσ(x) by means of classical Nadaraya-Watson type estimators of the kindσ whereσ 2 i is a consistent estimator of the spot volatility at time t i .
Following Kenmoe and Sanfelici (2014), our choice forσ 2 i is based on the Fourier-Malliavin estimator of the spot volatility. The spot volatility at time t i is estimated at the opening of each day using the whole high-frequency exchange rate time serieŝ where the k-th Fourier coefficient of the variance process is obtained by Bohr-convolution Figure 6 plots the histogram of the distribution of the EUR/USD exchange rate and the nonparametric estimate of the diffusion function of the exchange rate model obtained by the estimator (24). The function σ 2 (x) is plotted on a uniform grid x 0 < x 1 < ... < x J in the interval [0.08, 0.15]. The diffusion function does not look flat or linear, rather it shows the shape of a "smile" with higher volatilities for low exchange rates.

SECOND STEP:
The estimated volatility functionσ(·) is then used to generate N realizations of the process (22) and to price a given contingent claim in a Monte Carlo simulation. We generate 100 trajectories by dividing the interval [0, T] into n = 1000 subintervals and by discretizing the estimated diffusion process by first order Euler method.
THIRD STEP: For each generated trajectory, we estimate the spot functions A(t), B(t) and C(t) using the Fourier estimation metho-see Mancino et al. (2017)-and we compute the feedback rate function λ(τ) by (A6) in Remark A1: We remark that each time the Fourier cross-volatility estimation procedure is reiterated the resolution at which the estimated function can be observed becomes lower and lower. To maintain a good resolution, we can exploit the fact that, in the trajectories simulation process, the spot volatility A(t) = σ 2 (exp(X t )) is actually known. Therefore, we can easily compute its discrete Fourier coefficients by From c k (A n ) and c k (dX n ), the functions B(t), C(t) and λ(t) can be estimated by the Fourier methodology according to Theorem A2. 4 Furthermore, we highlight that any other approach for computing the Malliavin weight is not possible in this context because the computation of the variation process from (6) requires the knowledge of the functional form of the volatility function to perform differentiation.
FOURTH STEP: We compute via Monte Carlo method, by averaging the argument of the expectation over the simulated sample paths. As a benchmark for our computations, we solved the associated generalized Black-Scholes equation

∂V ∂t
with the final condition V(S, T) = f (S) by (backward Euler) finite difference method. Following Tavella and Randall (2000), we impose linearity conditions on the boundary of a truncated computational domain. The option delta values are also compared to those obtained by Monte Carlo finite difference and by the Black-Scholes formula with volatility parameter set equal to σ =σ(log(S)).
In Table 6, we show the results obtained for the pricing of an European Put option with strike K = 1.1 and expiry time T = 1 year. For simplicity, the risk free interest rate is taken to be r = 0. The perturbation of the initial datum in the Monte Carlo finite difference method is set to 0.001. In particular, as discussed in Mancino and Sanfelici (2008), the cutting frequencies for the functions B(t) and C(t) can be taken as N B = N/2 and N C = N B /2, respectively. Finally, the feedback rate function λ can be estimated at N C distinct points in the interval [0, T] by taking M C = N C /2 in the Fourier expansion for B(t). Being N = n/2, this procedure allows for a reconstruction of the feedback rate function at 125 nodes in the interval [0, T].
We notice that the lambda method outperforms the Monte Carlo finite difference method for at-the-money options, due to the lack of regularity of the payoff function f at S = K. However, the Monte Carlo finite difference method works better for out-of-the-money and deep out-of-the-money options and has in general a much smaller variability.

Conclusions
We have proposed a possible way to compute the hedging Delta using the Malliavin weight method of Fournié et al. (1999) in a local volatility framework. The advantage of our λ-method are twofold: on the one hand the advantage is computational, because it allows us to compute the Malliavin weight by a single stochastic integral versus a double stochastic integral as in Fournié et al. (1999). This yields a faster convergence rate in a Monte Carlo experiment, especially for discontinuous payoffs. On the other hand, our approach is nonparametric, as we only assume a general local volatility model and we can estimate the volatility and the other processes involved in the Greek formula with quantities that can be nonparametrically obtained from observed prices. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This Section resumes the mathematical computation of the λ function. The proof of the Theorems A1 and A2 can be found in Malliavin and Thalmaier (2006). Theorem A1. Assume that the asset price evolves according to the SDE dS t = a(S t )dW(t) + b(S t )dt, with initial value s 0 , where a and b are smooth functions. Consider the associated variation process dζ t = ζ t a (S t )dW t + b (S t )dt , ζ 0 = 1.
Then, the rescaled variation defined by z t = ζ t a(S t ) is differentiable with respect to t and it holds that Theorem A2. Denoting by , the quadratic (co-)variation operation, define the following cross-volatilities: dS t , dS t := A t dt , dA t , dS t := B t dt , dB t , dS t := C t dt , then the function λ t in (A3) has the following expression Remark A1. In particular, if b(S t ) = r − 1 2 a 2 (S t ) and substitute in (A3), then Then, the formula (A4) becomes