POLYNOMIAL CHAOS EXPANSION APPROACH TO MALLIAVIN CALCULUS ANALYSIS OF BOND OPTIONS SENSITIVITY

The present paper provides a sensitivity analysis of bond options exploiting the probabilistic properties of Malliavin Calculus and the computational benefits of the Polynomial Chaos Expansion. The purpose is to use the integration by parts formula of Malliavin Calculus in order to obtain the Greeks, of a given financial derivative, in a form which is suitable for numerical simulation. In particular, such computations will be performed usign the so called Polynomial Chaos Expansion technique, then comparing obtained results versus those retrieved using both the usual Monte-Carlo approach and the analytical formula. AMS Subject Classification: 33C90, 41A10, 42A61, 45D05, 60H10, 60H30, 65C30, 91G30, 91B16, 91B60, 91G20


Introduction
The Malliavin stochastic Calculus, is a powerful probabilistic tool first developed by Paul Malliavin in the 70's.It is particularly useful to study regularity bounds for the density of solutions of stochastic differential equations.In our case we shall make use of the so called Malliavin derivative which is, roughly speaking, a derivative with respect to a Brownian Motion.The use of Malliavin derivative will greatly simplifies the calculation of the sensitivity parameters characterizing financial derivatives, the main advantage being that it can be applied to the expectation of the contingent claim under the equivalent martingale measure, producing the same result obtained through the well known Malliavin weighted scheme, see [15], that is Greek = E[Xπ], where X is a random variable and π the so called Malliavin weight, which can be random as well.In this way the Greek is nothing else that the expectation of a random variable.It is worth to mention that such an approach can be extended to rather complex financial models, also allowing to obtain insights concerning (stochastic) control issues that can be found in several fields of applications, spanning from engineering, to neurobiology, etc., see, e.g.[4,20] and references therein.The main advantage of these new approach has been shown by Farai in [14] and Benhamou in [5] for cases where the distribution of the random variable inside the expectation is not known or the payoff of the derivative is not differentiable.We will not present such cases but a standard one in order to make possible the use of the second tool, the PCE.
The Polynomial Chaos Expansion (PCE) techniques, allow us to express a squared integrable random variable as linear combination of Hermite or Legendre polynomials.In this way the simulation of the random variable itself become easier and faster, especially concerning the estimate of its mean and variance, see, e.g.[13].This approach will provide a way to estimate the Greeks and will be compared to the Monte-Carlo approach using as benchmark the analytical value.

The Malliavin Calculus
We start giving a brief review about the Malliavin Calculus introduced by Nualart in [21].Such approach allows us to define a new kind of derivative, namely the Malliavin derivative, which could be used to efficiently compute the so called Greeks of a given financial model.We recall that a financial derivative, at time T > 0 pays off Φ(u(T, Π)), with Φ the pay-off function, u(T, Π) the derivative value at time T and Π a set of parameters.Then the derivative price at time t = 0 is u(0, Π) and the Greek is defined as its derivative with respect to one of the parameters α ∈ Π, i.e.

The Wiener Space
Let us consider a filtered probability space (Ω, F, F = {F t } {0≤t≤T } , P ) and a Brownian motion W = {W t } t∈[0,T ] with T ∈ R + , T < ∞, defined on it, assume in addition that F = {F t } {0≤t≤T } is the filtration generated by the Brownian motion W t .We recall that, almost surely for every ω ∈ Ω, W (ω Our aim is to treat the universe space Ω as a vector space, so exploiting the continuity property of the Brownian motion paths, for a.e.ω ∈ Ω, we identify With t varying between 0 and T .According to [14, Sec.2.1], in force of (1), we can consider the space Ω as the space of continuous functions from [0, T ] to R, that is Ω = C([0, T ], R).The advantage deriving from this identification is the possibility of defining a concept of derivative with respect to ω ∈ Ω in a similar way to the usual vector space definitions.
Let us consider the Hilbert space ( Let S n be the space of all polynomials of degree n with domain the set W (H).
In addition, let S be space of all random variables F such that where

The Malliavin derivative
Taking into account the definitions given in Sec.2.1, let us consider the random variable defined as in eq. ( 3), namely where h ∈ H is a deterministic function with H = L 2 (0, T ) from now on.Consider a function γ ∈ Ω.It will act as direction function in the derivation definition and this is possible since we have identified the universe space Ω with the space of continuous functions C([0, T ]).We will define γ as follows Definition 1.Let F : Ω → R be a random variable of the form as in eq. ( 4) and let γ ∈ Ω be defined as in eq. ( 5).Then the directional derivative D γ F of F at the point ω ∈ Ω, in the direction γ ∈ Ω is given by if the limit exists.
Let us define the space which is called the Cameron-Martin space.In addition, note that it is not possible to proceed in the Then, by Def.(1) Since dγ(t) = g(t)dt from (5), we get Where of course •, • H is the inner product on H, which here is assumed to be L 2 ([0, T ]).Finally, recalling that h ∈ H is a deterministic function, taking γ(t) = t and thus g(t) = 1, we can write From now on, the derivative in the direction γ(t) = t = id L 2 ([0,T ]) will be denoted by D. In general, the notation to represent the directional derivative of F in the direction γ(t) = t 0 g(s)ds is Definition 2. The derivative with ∂f ∂x i the derivative with respect to the i-th component, h i ∈ H and D is called Malliavin derivative on S.
So far we have derived in the Malliavin sense only functions belonging to the space S of the form (3). Let us now introduce the following norm:

The Skorohod integral
We introduce now the divergence operator, which is defined as the adjoint of the derivative operator.In particular, when the Hilbert space H is an L 2 space of the form L 2 (T, B, µ) with µ a σ-finite measure, we interpret the divergence operator as a stochastic integral and we call it the Skorohod integral which, in the brownian motion case is a generalization of Ito Integral, see [21,Sec. 1.3].The Malliavin derivative D is a closed linear operator defined on D 1,2 and Note that D is an unbounded operator with domain dense in L 2 (Ω) and its domain and range belong to different Hilbert spaces.From theory we have that the adjoint of unbounded operators have the same domain and range.In the following a definition of the adjoint operator δ of D is given.
where c is some constant depending on u.If u belongs to Dom(δ), then is the element of L 2 (Ω) such that the integration by parts formula holds: Moreover, δ is linear, densely defined and we call it Skorohod integral.
An important property of the Skorohod integral that makes its computation easier is the fact that its domain Dom(δ) contains all the adapted stochastic processes which belong to L 2 (Ω × [0, T ]).For these processes the Skorohod integral coincides with the Ito stochastic integral as shown in the following proposition, for more details see, e.g.[14,Sec 2.3].

An integration by parts formula
In this article the Mallivin derivative is used in the in order to simplify the calculation of the Greeks.The next proposition (from [21, Sec.6.2]) plays a key role in the derivation of quantities of the form E[f (F )G] with F, G random variables and f a continuously differentiable function.
Proposition 6.Let F , G be two random variables such that F ∈ D 1,2 .Consider an H-valued random variable u such that D u F =< DF, u > H = 0 a.s. and Gu(D u F ) −1 ∈ Dom(δ).Then, for any continuously differentiable function f with bounded derivative we have Where δ(u) is the Skorohod integral of u and D u F is the Malliavin derivative in the direction u.
Proof.Let us first apply the chain rule (Prop.3), so we get Since by hypotesis DF, u H = 0 we have Consequently, we have Now, since Gu( DF, u H ) ∈ Dom(δ), an application of eq. ( 18) gives which proves the Proposition.

Generalized Polynomial Chaos
In what follows we shall consider the Generalized Polynomial Chaos (gPC), which will be then exploited to deal with the PCE in the case of normal and uniform random variables.As gPC identifies a set of polynomials whose unknown is a given random variable, denoted by ξ.It follows that the elements of gPC can be seen as functionals constituting an orthonormal basis of a Hilbert space (L 2 (Ω, Σ, P ) in our case).However, not every random variable ξ defines a gPC basis, but as long as the basis is defined, it is possible to represent the elements of the Hilbert space as linear combinations of the gPC components, so we have an expansion, called Polynomial Chaos Expansion.In this way initial element of interest can be approximated by truncating the expansion at a desired order, which we will call N .To clarify the notation, from now on we indicate the elements of the gPC basis as {Ψ n (x)} n∈N the family of Hermite or Legendre polynomials.We would like to underline that, from a more theoretical point of view, other (approximations) approaches can also be pursued, as, e.g., exploiting techniques outlined in [19], or studying related asymptotics for the involved financial quantities, see, e.g., [2,10] and references therein.The latter approach having the advantage to treat directly small perturbations arising in financial markets, particularly with respect to the consideration of limiting cases for the parameters, as, e.g., w.r.t. the strike, the time to maturity, etc., involved in the dynamic of the product we are interested in.

One Dimensional Polynomial Chaos Expansion
Let (Ω, Σ, P ) be a probability space, with Ω an abstract set of elementary events, Σ is a σ-algebra of subsets of Ω, and P is a probability measure on Σ. Definition 7. The space L 2 (Ω, Σ, P ) is the Hilbert space of scalar realvalued random variables X : Ω → R defined on (Ω, Σ, P ) such that where dP (ω) = f (x)dx and f (x) is the probability law of X and the weight function which defines the scalar product of the Hilbert space.It follows that L 2 (Ω, Σ, P ) is a Hilbert space endowed with the scalar product with associated norm In order to simplify the notation we set ||X|| 2 P := ||X|| 2 L 2 (Ω,Σ,P ) and the convergence, in this norm is always referred as mean square convergence or strong convergence.
Let us state the two main elements to define the PCE scheme.The first is the random variable of interest that we want to analyse.From now on we will call it Y and it belongs to L 2 (Ω, Σ, P ).The other important element is the class of basic random variables ξ, which belongs as well to the space L 2 (Ω, Σ, P ).In particular, they are used as entries of the polynomials that compose the gPC basis.These basic functions ξ : Ω → D, with D ⊂ R, must satisfy the following properties: • ξ has finite moments of all orders.
The two types of basic random variables used from now on are the uniform U (−1, 1) and the normal N (0, 1).
A distinction between the basic random variable that rules the decomposition and the generic elements Y that one wishes to decompose needs now to be done.
Let us denote by σ(ξ) the σ-algebra generated by the basic random variable ξ.We have σ(ξ) ⊂ Σ.Having in mind to express the random variable Y in terms of ξ using a polynomial decomposition, it needs at least to be measurable with respect to the σ-algebra σ(ξ) generated by ξ.
By the Doob-Dynkin Lemma, see e.g.[22], Lemma 1.2.1, we have that Y is σ(ξ)-measurable, if there exists a Borel measurable function g : R → R, such that Y = g(ξ).Therefore we restrict our analysis to compute the decomposition with respect to the space L 2 (Ω, σ(ξ), P ).
Given D the range of the basic random variable ξ and B(D) the Borel set generated by it.We have that ξ induces, via its cumulative density function, a measure dF ξ on (D, B(D)), endowed with a Borel σ-algebra on D. We have now two possible situations: , while dF ξ = ω(x)dx is nothing but the density function • if ξ is standard Gaussian random variable, namely ξ ∼ N (0, 1), then The classes of Legendre and Hermite polynomials {Ψ n (x)} n∈N are orthogonal with respect to the weight function P (ω) defined above, which represents the probability density function p of the random variable ξ ∼ U (−1, 1), respectively of the random variable ξ ∼ N (0, 1).In addition, the family {Ψ n (x)} n∈N is orthogonal with respect to the measure dF ξ , in the Hilbert space L 2 (D, B(D), dF ξ ).
Considering ξ still as a standard normal or uniform random variable, the sequence {Ψ n (x)} n∈N is an orthogonal system of functional polynomials in Hilbert space L 2 (Ω, σ(ξ), P ) since inherits the orthogonality property of {Ψ n (x)} n∈N in L 2 (D, B(D), dF ξ ).There is also the fact that ψ 0 is considered as a degenerate random variable, also referred as an almost surely constant random variable, i.e.
Definition 8. Given a basic random variable ξ, the associated orthogonal system {Ψ n (x)} n∈N is called generalized polynomial chaos (gPC) basis for the space L 2 (Ω, σ(ξ), P ).Moreover if Y ∈ L 2 (Ω, σ(ξ), P ) its PCE reads as follows where, for all i ∈ N, it holds Theorem 9. Let ξ be a basic random variable with associated gPC basis {Ψ n (ξ)} n∈N .If Y ∈ L 2 (Ω, σ(ξ), P ), then the truncated series, at level N ∈ N + of the PCE expansion defined by eq. ( 25) reads as follows where the coefficients are defined as in equation (26).The sequence of random variable {Y (N ) } n∈N converges, in mean square sense, to Y as N → +∞.
Since the set of orthogonal polynomials {Ψ n (x)} n∈N is a maximal system in L 2 (D, B(D), dF ξ ), see, e.g., [13] and references therein, then, for every real, ǫ > 0 there exists a N ∈ N + such that

Properties of gPC basis
In the same general framework defined so far, by the fact that the basic random variable ξ induces the measure dF ξ on (D, B(D)), some properties of Legendre and Hermite polynomials are transferred to the gPC basis in the space L 2 (Ω, σ(ξ), P ).In particular we have • Property 1: which follows by straightforward calculations, namely, if i = 0 we have otherwise, by orthogonality with respect to the Ψ 0 polynomial, we obtain • Property 2: for each i ∈ N, let us define then {η i Ψ i (ξ)} n∈N is an orthonormal and maximal system in L 2 (Ω, σ(ξ), P ) and, applying the Parceval identity, we have where d i is the i-th Fourier coefficient, i.e.
Moreover, we have that where c i is the i-th coefficient of the PCE decomposition, see (26), therefore and we can derive a precise formula for the second raw moment of Y = g(ξ).In addition, the estimation of the mean square error, taking into account the N -th approximation level, see eq. ( 27), equals P which follows from the basic properties of orthogonal projection in Hilbert space, see e.g.[3,Sec. 14].Eventually, using the Parceval identity and the relation between d i and c i coefficients, see eq. ( 30), we have defined as the mean square error of the truncated expansion of degree N.
• Property 3 Given Y ∈ L 2 (Ω, σ(ξ), P ), and using the definition of the coefficients c i , along with the fact that Ψ 0 = 1, we have Moreover, we have

Delta of Bond Options in the two-Additive-Factor Gaussian Model
In what follows we will study the Greeks computed on the price of a European option written on a zero coupon bond whose yield curve is described by a two-Additive-Factor Gaussian Model, as defined in [7,Sec. 4.2], where the authors exploit such a model for interest rates.

The Short-Rate Dynamics
Let us consider the following dynamic under the risk neutral measure Q: where the processes {x(t) : t ≥ 0} and {y(t

Pricing a European Option on a Zero-Coupon Bond
At time t, the price of a European call option with maturity T and strike K, written on a zero-coupon bond with unit face value and maturity L is given by, see, e.g.[7,Sec. 4.2].Where the three aforementioned times t, T and L are such that: In order to explicitly compute expectation appearing in eq. ( 37) we need to change probability measure to an equivalent martingale measure and, in our case, this measure is the so called forward measure, denoted by Q T , for any fixed maturity time T , see [8, Sec.2.4] w.r.t. the Vasicek model.To have a little insight of what this measure is, observe that the Radon-Nykodim derivative used to pass to the risk-neutral measure Q, reads as follow The following lemma shows how the dynamic of the interest rate changes under the new measure Q T , see [7,Sec. 4.2] for more details.
Lemma 1.The processes x and y, defined in eq. ( 35), under the forward measure Q T , evolve according to where W T 1 and W T 2 are two correlated Brownian motions under Moreover, the explicit solutions of equations 39 are, for s ≤ t ≤ T , given by where So that, under Q T , the distribution of r conditional on F s is normal with known mean and variance.In addition, the conditional expectations of x and y are: Theorem 10 (Analytical formula for Option price on a zero-coupon bond).The price at time t of a European call option expiring at T , with strike K, and written on a zero-coupon bond with maturity L is: where while N is the cumulative distribution function of a standard Gaussian random variable. Proof.See kP (0, T ) Σ(0, T, L)P (0, L)kP (0, T ) having that   ln P (0,L) The formula (47) is the one that we will use later to compute the exact value of the Delta.

Malliavin Delta
To apply the Malliavin weighted scheme we need a formula for P (T, L) depending on the initial bond value P (0, L), which is exactly the result stated in [7, Sec.4.2, Thm.4.2.2], which claims that Proposition 11.Maintaining the same assumption made in Sec.3.2.3, the price of a zero-coupon bond at time T , with unit face value and maturity L, with 0 ≤ T ≤ L, is moreover under the forward measure Q T , the logarithm of P (T, L) conditional on F t is normally distributed with mean and variance Finally, under Q T , the price of the option is where Φ is the payoff function.
Proof.For the proof, see [7], Appendix C of Chapter 4.
Proposition 12.The Delta for a European Option with expiry time T, written on a zero-coupon bond with unitary face value and maturing at time L (with 0 ≤ T ≤ L) is: Proof.The idea is to apply the integration by parts formula stated in eq. ( 20), in order to avoid the differentiation of the payoff and get a convenient formula for simulation.
Omitting the conditioning on F 0 , we start from differentiating the formula (51) with t = 0, w.r.t.P (0, L): T ] using integration by parts formula (20) with Since the argument of δ is deterministic, from Prop.( 5) And thus, the thesis.

Numerical simulation
The purpose of this section is to investigate the accuracy of the PCE simulation over the classical Monte-Carlo, both of them applied on the Malliavin formula The results in terms of accuracy of the Delta estimation are displayed in Fig.
( Observing the results, we have performed three experiments for each approach in order to obtain three different levels of accuracy and easily compare PCE and Monte Carlo.For the first two levels (10 −2 and 10 −3 ) we can observe that the computational time cost of the PCE approach is bigger than the Monte Carlo.If the accuracy required increases to 10 −4 , the Monte Carlo method becomes, as expected, slower than PCE and we need 33 554 432 number of MC simulations and a PCE degree of 35.Note that the cost of the PCE includes the time spent in computing all the coefficients.Therefore, although to reach an order of accuracy of 10 −4 a big number of MC simulations and PCE degree are needed, the computational time required is negligible, making both methods reliable and precise.

Conclusion
This study produced a good result in terms of combined use of the two techniques of PCE and Mallivin Calculus.In particular, we have shown the main advantage of the Malliavin Calculus, which is the possibility to derive the expression for the Greeks of a financial derivative without making the derivative of the payoff appear in the final formula.In this specific case, this approach was not necessary, but it can become useful whenever an analytical formula is not available (because the probability law of the process is unknown) or in those cases in which the payoff is not differentiable.Thus, an interesting development would be that of using these tools in more complex financial derivatives like European options with stochastic volatility or Exotic that have more complicated payoffs, like the Asian or Digital.Moreover it remains rather not analyzed the case of financial models aiming at capture the dynamics of structured derivatives characterized by jumps and abrupt changes in their sample paths.Concerning the latters, insights can be gained looking at the invariant measures associated to their theoretical counterpart, see, e.g., [1] and references therein, or [9] for a more financially oriented discussion.
A further line of research is related to the fine analysis of the possible symmetries, see, e.g., [11,12,16] and references therein, characterizing the particular problem into considerations.In fact, using the F eynman − Kac formula one can reduce the original stochastic problem into a deterministic (PDE) one, then exploit possible Lie symmetries to simplify the framework.Going back to the stochastic setting, one can newly use the PCE techniques, but with a considerable low degree for the involved orthogonal polynomials.
coming with the scalar product a, b H = T 0 a(t)b(t)dt and define W (H) := { T 0 h(t)dW t |h ∈ H}.

[ 7 ,
Appendix C ofChapter 4].So far the formula needed for the computation of the Delta has been stated, therefore taking into account eq.(45), and since its Delta is obtained by deriving at time t = 0, with respect the initial value of the underlying P (0, L), then for 0 = t ≤ T ≤ L, we have ∂ZBC(0, T, L, K) ∂P (0, L)

Figure 1 :
Figure 1: Accuracy (absolute error) and execution times (in seconds) for different configurations of Monte Carlo and PCE approaches.
Let D 1,2 be the closure of S in the norm || • || 1,2 .Can be shown that the derivative D on S can be extended to a closed operator on D 1,2 .An important result is the chain rule for the Malliavin derivative.
In our numerical experiments, we have considered different values of N , in order to see how the computational accuracy varies in dependence of such a parameter.Obtained results have been compared with the ones obtained by simple Monte-Carlo simulations with different sample sizes.

Table 1 :
[6]nd in Table(1).The hardware used is a CPU Intel Core i5, 2.20 GHz x 2 and 8 GB of RAM.The values used in the experiments are taken from a calibration performed in[6], on the British market: Accuracy (absolute error) and execution times (in seconds) for different configurations of Monte Carlo and PCE approaches.