On the Malliavin approach to Monte Carlo ... - Semantic Scholar

Report 2 Downloads 45 Views
On the Malliavin approach to Monte Carlo approximation of conditional expectations Bruno Bouchard

Nizar Touzi

Universit´e Paris VI

Centre de Recherche

and CREST

en Economie et Statistique

Paris, France

Paris, France

[email protected]

[email protected]

January 2002

Abstract Given a multi-dimensional Markov diffusion X, the Malliavin integration by parts formula provides a family of representations of the conditional expectation E[g(X2 )|X1 ]. The different representations are determined by some localizing functions. We discuss the problem of variance reduction within this family. We characterize an exponential function as the unique integratedvariance minimizer among the class of separable localizing functions. For general localizing functions, we provide a PDE characterization of the optimal solution, if it exists. This allows to draw the following observation : the separable exponential function does not minimize the integrated variance, except for the trivial one-dimensional case. We provide an application to a portfolio allocation problem, by use of the dynamic programming principle.

Key words: Monte Carlo, Malliavin calculus, calculus of variations. MSC 1991 subject classifications: Primary 60H07 65C05 ; secondary 49-00.

1

1

Introduction

Let X be a Markov diffusion process. Given n simulated paths of X, the purpose of this paper is to provide a Monte Carlo estimation of the conditional expectation r(x) := E [g(X2 )|X1 = x] ,

(1.1)

i.e. the regression function of g(X2 ) on X1 . In order to handle the singularity due to the conditioning, one can use Kernel methods developed in the statistics literature, see e.g. [2]. However, the asymptotic properties of the Kernel estimators depend on the bandwidth of the Kernel function as well as the dimension of the state variable X. Therefore, using these methods in a Monte Carlo technique does not induce the √ n rate of convergence. Malliavin integration by parts formula has been suggested recently in [6], [5], [9] √ and [8] in order to recover the n rate of convergence. Let us first discuss the case of the density estimator considered by [8]. The starting point is the expression of the density p of a smooth real-valued random variable G (see [11] Proposition 2.1.1) : "

DG p(x) = E 1{G>x} δ |DG|2

!#

,

(1.2)

where δ is the Skorohod integration operator, and D is the Malliavin derivative operator. Writing formally the density as p(x) = E[εx (G)], where εx is the Dirac measure at point x, the above expression is easily understood as a consequence of an integration by parts formula (integrating up the Dirac). A remarkable feature of this expression is that it suggests a Monte Carlo estimation technique which does not require the use of Kernel methods in order to approximate the Dirac measure. The same observation prevails for the case of the regression function r(x) which can be written formally in : r(x) =

E [g(X2 )εx (X1 )] . E [εx (X1 )]

Using the Malliavin integration by parts formula, [5] suggest an alternative representation of the regression function r(x) in the spirit of (1.2). This idea is further developed in [9] when the process X is a multi-dimensional correlated Brownian motion. An important observation is that, while (1.2) suggests a Monte Carlo estimator √ with n rate of convergence, it also provides the price to pay for this gain in efficiency : the right-hand side of (1.2) involves the Skorohod integral of the normalized 2

Malliavin derivative |DG|−2 DG; in practice this requires an approximation of the continuous-time process DG and its Skorohod integral. In this paper, we provide a family of such alternative representations in the vectorvalued case. As in [5], we introduce localizing functions ϕ(x) in order to catch the idea that the relevant information, for the computation of r(x), is located in the neighborhood of x. The practical relevance of such localizing functions is highlighted in [9]. The main contribution of this paper is the discussion of the variance reduction issue related to the family of localizing functions. We first restrict the family to the Q class of separable functions ϕ(x) = i ϕi (xi ). We prove existence and uniqueness of a solution to the problem of minimization of the integrated variance in this class. i i The solution is of the exponential form ϕi (xi ) = e−η x where the η i ’s are positive parameters characterized as the unique solution of a system of non-linear equations. In the one dimenstional case, this result has been obtained heuristycally by [8]. We also study the problem of minimizing the integrated variance within the larger class of all localizing functions. Assuming existence, we obtain a partial differential equation with appropriate boundary conditions. This is a characterization of the solution, if it exists, since the minimization problem is strictly convex in this setting. In the general multi-dimensional case, the existence of a solution to this equation is left is still an open problem. However, an interesting observation is that separable localizing functions do not solve this equation. The paper is organized as follows. Section 2 introduces the main notations together with some preliminary results. Section 3 contains the proof of the family of alternative representations of the conditional expectation. The variance reduction issues are discussed in Sections 4 and 6. Numerical experiments are provided in Section 7. Finally, Section 8 provides an application of this technique to a popular stochastic control problem in finance, namely find the optimal portfolio allocation in order to maximize expected utility from terminal wealth.

2

Preliminaries

Let (Ω, F, P ) be a complete probability space equipped with a d-dimensional standard Brownian motion W = (W 1 , . . . , W d ). Since we are interested in the computation of the regression function (1.1), we shall restrict the time interval to T := [0, 2].

3

We denote by IF := {Ft , t ∈ T} the P −completion of the filtration generated by W . Throughout this paper, we consider a Markov process X defined as the unique strong solution of the stochastic differential equation : dXt = b(Xt )dt + σ(Xt )dWt ,

(2.1)

together with an initial condition X0 . Here, b, σ and σ −1 are Cb∞ vector and matrixvalued functions. Under the above conditions, for all t ∈ T, Xt belongs to the Sobolev spaces IDk,p (p, k ≥ 1) of k− times Malliavin differentiable random variables satisfying : ||Xt ||IDk,p



:= E(|Xt |p ) +

k X

1

p



E ||D

j

Xt ||pL2 (Tj )



j=1



0, σ 22 > 0. An admissible strategy is a U -valued predicable process (U is some compact subset of IR). We denote by U the set of such processes. Given a strategy ν ∈ U, the corresponding wealth process Y 0,y,ν is the solution of dYt = νt

dXt1 , Xt1

with the initial condition Y00,y,ν = y. Since the contingent claim can not be perfectly hedged, we intend to find a partial hedging strategy defined so as to maximize the expected exponential utility 

G

−a(YT0,y,ν −G)

v (0, y, x) := sup E −e



,

ν∈U

where a > 0 is the risk aversion parameter and T > 0 is a given time horizon. Comparison with the maximal expected utility in case where the option is not delivered 

0

−a(YT0,y,ν )

v (0, y, x) = sup E −e



ν∈U

leads to the well-known definition of the indifference utility valuation rule (see e.g. [7] ) n o p(G, x) := inf π ∈ IR : v G (0, π, x) ≥ v 0 (0, 0, x) . Observing that v G (t, y, x) = e−ay v G (t, 0, x) , we see that

1 v 0 (0, 0, x) p(G, x) = − ln G a v (0, 0, x) 27

(8.1) !

.

Hence, the computation of p is reduced to the computation of the involved value functions. The value function v G satisfies the dynamic programming equation 

h

t,y,ν v G (t, y, x) := sup E v G t + ∆, Yt+∆ , Xt+∆ ν∈U



| Xt = x

i

∆ > 0.

Using again the identity (8.1), the above equation may be written in t,0,ν

h

i

v G (t, 0, x) := sup E e−aYt+∆ v G (t + ∆, 0, Xt+∆ ) | Xt = x ; ν∈U

∆ > 0,

so that it is sufficient to compute the value function for y = 0. We also observe that v is independent of x1 . Omitting this variable in the definition of v G , we simply write : t,0,ν

h



2 v G (t, 0, x2 ) := sup E e−aYt+∆ v G t + ∆, 0, Xt+∆ ν∈U



i

| Xt2 = x2 ;

∆>0.

(8.2)

Hence, in the context of the particular model studied in this section, the number of state variables is reduced to one, which considerably simplifies the computations. The reason for considering such simplifications is that we are mainly concerned by the dynamic application of the conditional expectation estimation presented in the first sections of this paper.

8.2. The case G = 0. When there is no contingent claim to be delivered, the value function v 0 does not depend on the state variable X 2 . It follows from (8.2) that the optimal control process is constant in time. The Monte Carlo estimation procedure is then considerably simplified, as it is sufficient to perform it on a single time step.

8.3. Discretization in time for v G . Set tk := kT /n, n ∈ IN and k = 0, . . . , n. The above dynamic programming equation leads to : v

G

(tk , 0, Xt2k )



t ,0,ν

−aYt k

:= ess sup E e

k+1

v

G

ν∈U



tk+1 , 0, Xt2k+1



|

Xt2k



; 0≤k ≤n−1.

Restricting the control process to be constant on [tk , tk+1 [, we obtain the following approximation G



(tk , 0, Xt2k )



:= ess sup E e ν∈U

t ,0,ν

−aYt k

k+1

G





tk+1 , 0, Xt2k+1



|

Xt2k



; 0 ≤ k ≤ n − 1 .(8.3)

We now can appeal to the conditional expectation representation studied in this paper. Using a slight adaptation of Corollary 3.1 (see (3.3)), it is easily checked 28

that, for each ν ∈ U , 

t ,0,ν

−aYt k

E e

G



k+1

tk+1 , 0, Xt2k+1 t ,0,ν



=



−aYt k

E 1x2 <Xt2 e

G



k+1

k







|

Xt2k

=x

tk+1 , 0, Xt2k+1



2

 

∆k

ϕ(Xt2k



2

−x )





(8.4)

E 1x2 <Xt2 ∆k ϕ(Xt2k − x2 ) k

where ϕ(x) = e−ηx for some η > 0 and tk

Z

∆k (·) =

0

·hk dWt2 −

Z

tk+1

¯ k dW 2 ·h t

tk

¯ k = [(tk+1 − tk )σ 22 X 2 ]−1 . with hk = [tk σ 22 Xt2k ]−1 and h tk We can then use a backward induction procedure to compute the value function as in [3] (see also [9] for the case of optimal stopping problems). 8.4. Monte Carlo approximation. We start by simulating N paths, X i , of the process X. Given the value function for k = n (i.e tn = T ), 2,i

vˆG (tn , 0, Xt2,i ) = −e−a(0−g(Xtn )) n we then use the other simulated paths, X j , j 6= i, in order to approximate (8.4) at the point (tn−1 , Xt2,i ) by n−1 P

tn−1 ,0,ν,j

1n j6=i

X

2,i 2,j <X tn−1 tn−1

P

j6=i

where

t ,0,ν,j Ytnn−1



Xt1,j n

Xt1,j

o e−aYtn

1n

X

2,i 2,j <X tn−1 tn−1

!

−1



vˆG (tn ,0,Xt2,j ∆j ϕ(Xt2,j n ) k

o ∆j

k



ϕ(Xt2,j

n−1

−Xt2,i

n−1

n−1

)

−Xt2,i

n−1

)





and ∆jk is the Skorohod integral operator corre-

n−1

sponding to the j-th path, i.e. : 



∆jk ϕ(Xt2,j − Xt2,i ) n−1 n−1 = ϕ(Xt2,j n−1



Xt2,i ) n−1

Wt2,j /tn−1 −(Wt2,j −Wt2,j )/(tn −tn−1 )+σ 22 n n−1

n−1 σ 22 Xt2,j n−1



!

(see Remark 3.4). The optimization over the parameter ν is achieved by a simple Newton-Raphson algorithm. This produces of vˆG (tn−1 , 0, Xt2,i ). Iterating this pron−1 2,i cedure backward, we can estimate, for each k and i, vˆ(tk , 0, Xtk ) together with the corresponding optimal control νˆ(tk , 0, Xt2,i ). k

29

8.5. Numerical experiments. We consider the contingent claim defined by the payoff function g(y) = 5 ∗ min{K, y}. We fix x = (1, 1)∗ , µ = (0.1, 0.1)∗ , σ 11 = 0.15, U = [0, 40], T = 1, n = 10, and we perform the computations for different values of K, σ 22 , σ 21 and a. Each estimation of the value functions v G , and the induced price pG , is based on 8192 simulated paths. For each experiment, we compute 200 estimators. The average and the standard deviation in pourcentage of the average (in brackets) are collected in the following table. K = 1.2, σ 21 = 0.1, σ 22 = 0.1 pG vG a = 0.25 5.20 [0.46%] −2.94 [0.55%] a=1 5.26 [1.54%] −154.04 [14.36%] K = 1.2, σ 21 = 0.05, σ 22 = 0.2 pG vG a = 0.25 5.26 [0.50%] −2.98 [0.56%] a=1 5.40 [0.31%] −177.91 [1.66%] K = ∞, σ 21 = 0.05, σ 22 = 0.2 pG vG a = 0.25 5.59 [0.52%] −3.23 [0.61%] a=1 6.29 [1.30%] −433.25 [8.91%] K = ∞, σ 21 = 0.15, σ 22 = 0.005 pG vG a = 0.25 5.09 [0.22%] −2.85 [0.39%] As expected, the price is increasing with the risk aversion parameter a and with K. It is decreasing with the ”relative correlation” between S 1 and G. For K = ∞, G = 5 ∗ X 2 which is not ”far” from 5 ∗ X 1 for σ 21 = σ 11 >> σ 22 . For a low risk aversion parameter we should then obtain a price of order of 5 ∗ X01 = 5 which is nearly what we get for a = 0.25.

References [1] V. Bally and G. Pag`es (2001). A quantization algorithm for solving multidimensional optimal stopping problems, preprint. [2] D. Bosq (1998). Non-parametric Statistics for Stochastic Processes, Springer Verlag, New York. 30

[3] B. Bouchard (2000). Contrˆole Stochastique Appliqu´e `a la Finance, PhD thesis, Universit´e Paris Dauphine. [4] J. F. Carri`ere (1996). Valuation of the Early-Exercise Price for Options using Simulations and Nonparametric Regression, Insurance : mathematics and Economics, 19, 19-30. [5] E. Fourni´e, J.-M. Lasry, J. Lebuchoux and P.-L. Lions (2001). Applications of Malliavin calculus to Monte Carlo methods in finance II, Finance and Stochastics 5, 201-236. [6] E. Fourni´e, J.-M. Lasry, J. Lebuchoux, P.-L. Lions and N. Touzi (1999). Applications of Malliavin calculus to Monte Carlo methods in finance, Finance and Stochastics 3, 391-412. [7] S. Hodges, A. Neuberger (1989). Optimal Replication of Contingent Claims under Transaction Costs. Review of Futures Markets 8, 222-239. [8] A. Kohatsu-Higa and R. Pettersson (2001). Variance reduction methods for simulation of densities on Wiener space, preprint. [9] P.-L. Lions and H. Regnier (2001). Calcul du prix et des sensibilit´es d’une option am´ericaine par une m´ethode de Monte Carlo, preprint. [10] F. A. Longstaff and R. S. Schwartz (2001). Valuing American Options By Simulation : A simple Least-Square Approach, Review of Financial Studies 14, 113-147. [11] D. Nualart (1995). The Malliavin Calculus and Related Topics Springer Verlag, Berlin.

31