Pierre Del Moral

Bruno Rémillard

Sylvain Rubenthaler

Del Pierre

Bruno R Moral

Sylvain Émillard

Rubenthaler

MONTE CARLO APPROXIMATIONS OF AMERICAN OPTIONS THAT PRESERVE MONOTONICITY AND CONVEXITY

Keywords: American option, American option, simulation, exercise region, Snell envelope

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Evaluation of American options is a central problem in financial engineering. Many ways of tackling the problem have been proposed so far, mostly for options on a single asset.

One of the first method was proposed by [START_REF] Brennan | The valuation of American put options[END_REF], and is based on numerical solutions of partial differential equations. Then tree-based methods were first introduced by [START_REF] Cox | Option pricing: A simplified approach[END_REF] for American options. Although Monte Carlo methods were proposed by [START_REF] Boyle | Options: A Monte Carlo approach[END_REF] for European options, it seems that the first simulationbased method is the one of [START_REF] Tilley | Valuing American options in a path simulation model[END_REF]. After that paper, Monte Carlo methods increased in popularity and many papers improved Tilley's results, e.g. [START_REF] Carriere | Valuation of the early-exercise price for options using simulations and nonparametric regression[END_REF], [START_REF] Broadie | Pricing American-style securities using simulation[END_REF], [START_REF] Longstaff | Valuing American options by simulation: A simple least-square approach[END_REF], [START_REF] Broadie | A stochastic mesh method for pricing highdimensional American options[END_REF].

However, even if in some interesting cases, e.g. American call-on-max or American puton-max options of dividend paying assets, where the underlying assets follows geometric Brownian motions, the value of an American option possesses convexity or monotonicity properties, the Monte Carlo methods described above do not preserve these properties. For a related work, see [START_REF] Laprise | Pricing American-style derivatives with European call options[END_REF].

In what follows, a Monte Carlo approach is proposed so that convexity and monotonicity are preserved. It is shown that the algorithm is quite precise for American options on a small number of assets. When the number of assets is large then the proposed method can be used to implement primal-dual methods, as in [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF].

In Section 2, one states the definition of the Snell envelope, together with a review of existing methods for pricing American options, in particular the methods of [START_REF] Carriere | Valuation of the early-exercise price for options using simulations and nonparametric regression[END_REF], [START_REF] Longstaff | Valuing American options by simulation: A simple least-square approach[END_REF] and [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF]. Properties of options prices in terms of assets dynamics and payoff functions are then studied in Section 3.1, and the algorithm for pricing an American option is described in Section 3.2. Finally, implementation issues are presented in Section 4, and a brief discussion of the results is provided in Section 5.

A brief review of algorithms for valuation of American options

As it is usually assumed, one wants to calculate the value of a Bermudan option, that is an option that can be exercised at given fixed periods, instead of a real American option, which can be exercised anytime. For simplicity, assume that the exercise periods are 0, 1 . . . n.

Suppose that, under the risk neutral measure, the value of the assets, including possibly stochastic volatility factors, is modeled by a discrete time Markov chain (X k ) k≥0 , and that the actualized payoff Z k at period k, given by Z k = f k (X k ), is integrable for any k ∈ {0, . . . , n}. For example, if the interest rate is constant and given by r, then f k (x) = f k x (1) , . . . , x (d) = e -rk max max 1≤j≤d x (j) -K, 0 is the actualized payoff of a so-called call-on-max option.

In what follows, one first states the solution of the American option problem, in terms of the Snell envelope in the Probability literature. In Operations Research, such a problem is called a dynamic programming problem. Then one discusses some deterministic and stochastic algorithms for calculating the Snell envelope or the exercise region.

2.1. Snell envelope. For a given filtration F, the value U 0 (x) at period 0 of the "American" option with X 0 = x, is given by

U 0 (x) = sup τ ∈T 0,n E(Z τ |X 0 = x),
where T k,n stands for the set of all F stopping times with values in {k, . . . , n}, 0 ≤ k ≤ n.

The sequence of functions

U k (x) = sup τ ∈T k,n E(Z τ |X k = x), k ∈ {0, . . . , n},
is called the Snell envelope.

It is well-known, e.g. [START_REF] Neveu | Discrete-parameter Martingales[END_REF], that the Snell envelope can be calculated in a recursive way. In fact, U n ≡ V n ≡ f n , and

V k (x) = E{U k+1 (X k+1 )|X k = x}, U k (x) = f k (x) ∨ V k (x) = max {f k (x), V k (x)} , k = n -1, . . . , 0.
Note that U k also has the representation

U k (x) = E(Z τ * k |X k = x),
where

τ * k = min{j ≥ k ; f j (X j ) = U j (X j )}, k = 0, 1, . . . , n. Also τ * k = min{j ≥ k ; f j (X j ) ≥ V j (X j )}, k = 0, 1, . . . , n.
An equivalent approach to solving the problem of American options is to look at the exercise region E, defined as

E = {(k, x) ; f k (x) = U k (x)} = {(k, x) ; f k (x) ≥ V k (x)}.
In fact, knowing E suffices to evaluate the option, since in that case, one can obtain the optimal stopping times τ * k viz.

τ * k = min{j ≥ k ; (j, X j ) ∈ E}, k = 0, 1, . . . , n.
In some situations, the boundary of the exercise region is quite simple, e.g. [START_REF] Broadie | The valuation of American options on multiple assets[END_REF]. For example, for the American call option in a Black-Scholes setting,

E = n k=0 {k} × [e k , ∞),
for some e 0 , . . . , e n , since (k, x) ∈ E implies (k, x ) ∈ E, for any x ≥ x. Therefore it suffices to know the e k 's to evaluate the option.

The boundary of the exercise region for the American put option in a Black-Scholes setting is also quite simple, and is given by

E = n k=0 {k} × [0, e k ],
for some e 0 , . . . , e n , since (k, x) ∈ E implies (k, x ) ∈ E, for any x ≤ x. However, in general, the exercise region is more complex. For more details, see [START_REF] Broadie | The valuation of American options on multiple assets[END_REF]. et al. (1979), and more recently the finite Markov chain methods, e.g. [START_REF] Duan | American option pricing under GARCH by a Markov chain approximation[END_REF], [START_REF] Duan | Numerical pricing of contingent claims on multiple assets and/or factor -A low-discrepancy Markov chain approach derivatives in the GARCH model[END_REF][START_REF] Bally | A quantization tree method for pricing and hedging multidimensional American options[END_REF]. In the latter, a powerful optimal method for choosing the approximation, known as quantization, is discussed. See also [START_REF] Kargin | Lattice option pricing by multidimensional interpolation[END_REF] for a similar approach that can be easily adapted to Bermudan and American options. Broadly speaking, the stochastic method proposed by [START_REF] Broadie | A stochastic mesh method for pricing highdimensional American options[END_REF] can be seen as a random version of the above algorithms. In all cases, theses algorithms are designed to approximate the value of U 0 (x) for a fixed x, not for all x.

Classes of algorithms for valuation of

(ii) Snell's method: In these algorithms, the Snell envelope U is calculated. The main representatives are those based on partial differential equations, e.g. [START_REF] Brennan | The valuation of American put options[END_REF], and Carriere's algorithm, e.g. [START_REF] Carriere | Valuation of the early-exercise price for options using simulations and nonparametric regression[END_REF], which is a simulation-based method. The latter will be described in more details below.

(iii) Exercise region method: For these algorithms, E is estimated. The first algorithm in this direction appeared in [START_REF] Tilley | Valuing American options in a path simulation model[END_REF] for options based on a single asset. It was later improved by [START_REF] Longstaff | Valuing American options by simulation: A simple least-square approach[END_REF]. Almost immediately, several papers appeared that validated the Longstaff-Schwartz method, e.g. [START_REF] Clément | An analysis of a least squares regression method for American option pricing[END_REF]. Other exercise region methods include the "primal-dual" methods of [START_REF] Rogers | Monte Carlo valuation of American options[END_REF] and [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF]. These methods are also discussed in more details below. For other methods for a single asset based on the exercise boundary, see e.g. [START_REF] Broadie | Nonparametric estimation of American options' exercise boundaries and call prices[END_REF].

The (random) algorithms of [START_REF] Carriere | Valuation of the early-exercise price for options using simulations and nonparametric regression[END_REF] and [START_REF] Longstaff | Valuing American options by simulation: A simple least-square approach[END_REF], all have in common the following features. They both rely on the simulation of a large number of paths, and conditional expectations given F k are approximated by functions in some given

classes A k , k = 0, . . . , n -1.
These algorithms, together with the dual approach algorithm, e.g. [START_REF] Rogers | Monte Carlo valuation of American options[END_REF] and [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], are detailed next.

Carriere algorithm.

Start by simulating N paths, denoted by

S 1 = (X 1,0 , . . . , X 1,n ), . . . , S N = (X N,0 , . . . , X N,n ).
Then, based on approximating classes A 0 , . . . , A n-1 , estimate V by V and U by Û using local regression methods. More precisely, set Vn = Ûn = f n , and then, for all n -1, . . . , 0, let Vk and Ûk be defined by (2.1) Vk = arg min

g∈A k N i=1 Ûk+1 (X i,k+1 ) -g(X i,k ) 2 , Ûk = f k ∨ Vk ,
provided the solution of the minimization problem (2.1) exists.

In [START_REF] Carriere | Valuation of the early-exercise price for options using simulations and nonparametric regression[END_REF], A k is the set of all q-splines.

= k=1 k × E k and E k = {x ; f k (x) ≥ Ṽk (x)}, with Ṽn = f n , τn ≡ n,
and for all k = n -1, . . . , 0,

(2.2) Ṽk = arg min

g∈A k N i=1 f τk+1 (X i,τ k+1 ) -g(X i,k ) 2 ,
where for any path S = (X 0 , . . . , X n ),

τk (S) = k if f k (X k ) ≥ Ṽk (X k ), τk+1 (S) if f k (X k ) < Ṽk (X k ),
provided the solution of the minimization problems (2.2) exist. Then, the estimation of

U 0 (X 0 ) is given by (2.3) Û0 (X 0 ) = 1 N N i=1 f τ0 (X i,τ 0 ).
As before, the same conclusions as those in Remark 2.1 apply.

Remark 2.2. The Carriere and Longstaff-Schwartz algorithms may seem equivalent, but they are not. In the latter one, emphasis is on stopping times. Instead of estimating U k by

f k ∨ Ṽk , as in Carriere's algorithm, U k (X i,k ) is estimated by (2.4) Ûk (X i,k ) = 1 N N i=1 f τk (X i,τ k ), k ∈ {0, . . . , n}.
So, in the Longstaff-Schwartz algorithm, U k is estimated at sample points, not all points.

2.5. Primal-dual approach. As seen before, to approximate the value of an American option, one can estimate directly E, as proposed in Section 2.4, or approximate first U or V , and then approximate E, as proposed in Section 2.3 or Section 3.2.

In [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], the authors proposes a two-stage approach: one for obtaining a lower bound for U 0 (X 0 ) and the other one for an upper bound. Suppose that Ê is an estimation of E, e.g. Ê is obtained by the methods discussed previously.

2.5.1. Algorithm for the lower bound.

A1: For each i = 1, . . . , N 0 , simulate a path S i = (X i,0 , . . . , X i,n ) starting from X 0 ; A2: Define the entrance time of S i in E as

τ E,i = min j ≥ 0 ; (j, X i,j ) ∈ E ; A3: Calculate (2.5) LB = 1 N 0 N 0 i=1 f τ E,i X i,τ E,i .
Then LB is a pointwise estimation of a lower bound for U 0 (X 0 ). One could also calculate the standard deviation σLB associated with the data

LB i = f τ E,i X i,τ E,i , 1 ≤ i ≤ N 0 .
One can assume that (0, X 0 ) ∈ E, for otherwise, Û0 (X 0 ) = f 0 (X 0 ).

Remark 2.3. As emphasized by [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], the last procedure yields a lower bound for U 0 (X 0 ). In fact, even if one can approximate U and V , to estimate U 0 (X 0 ) more precisely, it is recommended to calculate (2.5) in addition, since errors in U or E often transform in smaller errors for τ E , as many simulations showed. Unfortunately, so far there is no mathematical justification of that property.

For the upper bound, [START_REF] Rogers | Monte Carlo valuation of American options[END_REF] showed that for any martingale M ,

U 0 (X 0 ) ≤ M 0 + E max 0≤k≤n {f k (X k ) -M k } ,
with equality for the (unique) martingale M associated with the Doob-Meyer decomposition of the supermartingale U k (X k ).

Remark 2.4. In fact, the starting point of the martingale is not important. For simplicity, one could restrict ourselves to martingales starting from M 0 = 0. In that case, the equality

holds for M k = M k -M 0 ,
where M is the (unique) martingale M associated with the Doob-Meyer decomposition of the supermartingale U k (X k ).

Based on E, [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF] suggests a Monte Carlo algorithm to construct a martingale M , close to the optimal Doob-Meyer martingale. For all k ∈ {0, . . . , n}, let l k be the indicator function of the set (k,

X k ) ∈ E , let τ E,k = min{j ≥ k; (j, X j ) ∈ E} and define L k = E Z τ E,k |F k . Further set set M 0 = L 0 = Û0 (X 0 ), and for any 1 ≤ k ≤ n, set M k = M k-1 + L k -L k-1 -l k-1 E (L k -L k-1 |F k-1 ) .
Note that the last equation can be written in the much simpler form

M k = M k-1 + L k -E(L k |F k-1 ) (2.6) = M k-1 + E {f τ k (X τ k )|F k } -E {f τ k (X τ k )|F k-1 } . (2.7)
Since the conditional expectations must be estimated, they suggest a procedure to approximate M . However, based on (2.6)-(2.7), the following modification seems more natural.

For each simulated path S i , i ∈ {1, . . . , N 1 }, set M i,0 = 0, and repeat the following steps, for each k = 1, . . . , n:

B1:

-

For each j = 1, . . . , N 2 , simulate a new sub-path S i,j = {X i,j, } n =k , starting from X i,k ; -Calculate Li,k = 1 N 2 N 2 j=1 f τ E,k X i,j,τ E,k ; B2: -For each j = 1, . . . , N 3 , simulate a new sub-path S i,j = {X i,j, } n =k-1 , starting from X i,k-1 ; -Calculate Ê(L i,k |F k-1 ) = 1 N 3 N 3 j=1 f τ E,k X i,j,τ E,k ; B3: Set M i,k = M i,k-1 + Li,k -Ê(L i,k |F k-1 ).
Finally, the upper bound UB for U 0 (X 0 ) is then approximated by

(2.8) UB = 1 N 1 N 1 i=1 max 0≤k≤d {f k (X i,k ) -M i,k } .
One could also calculate the standard deviation σUB associated with the data

UB i = max 0≤k≤d {f k (X i,k ) -M i,k }, 1 ≤ i ≤ N 1 .
Because the lower and upper bounds estimations are conditionally independent given E, a 95% percent confidence interval for U 0 is given by (2.9)

LB -1.9545 σLB √ N 0 , UB + 1.9545 σUB √ N 1 .
Remark 2.5. In [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], it is argued that the upper bound constructed by their Monte Carlo methods will always be greater that the Monte Carlo lower bound, i.e., the value of (2.5) should be smaller than the value of (2.8). If true, their reasoning should apply as well to the modified algorithm described above. However their proof is clearly incorrect1 . As shown in Tables 12, it is possible that the (pointwise) Monte Carlo upper bound is smaller that the (pointwise) Monte Carlo lower bound.

Approximation of the Snell envelope

From now on, assume that the Markov chain (X k ) k≥0 takes values in a convex subset X of [0, ∞) d . Note that this setting excludes finite state space Markov chains, which can be treated much more easily. As customary, when s, t ∈ R d , s ≤ t means that s (j) ≤ t (j) for all j = 1, . . . , d. For x ∈ X, let x = max 1≤i≤d x (i) . Finally, recall that a real valued function g defined on X is non decreasing if for any x, x ∈ X, x ≤ x implies g(x) ≤ g(x ); conversely, g is non increasing if -g is non decreasing.

Throughout the rest of the section, one has to make some hypotheses. The first assumption is related to the law of the Markov chain.

Assumption 1. For all 1 ≤ k ≤ n, X k has the following representation

(3.1) X k = π k (X k-1 , Y k ), Y k ∈ Y,
where

Y k has law μ k , is independent of F k-1 , and x → π k (x, y) is continuous on X, for any fixed y ∈ Y.
Since Monte Carlo simulations play an essential role in the proposed methodology, the following assumption is needed too.

Assumption 2. For each k ∈ {1, . . . , n}, f k (X k ) is integrable.
That assumption insures that the Snell envelope is well defined.

Assumption 3. For all 0 ≤ j ≤ k ≤ n, the mappings x → E {f k (X k )|X j = x} and
x → f k (x) are continuous.

In the next section, desirable properties like convexity and monotonicity are studied for American options. The main result is that when the payoff has nice properties, then they are inherited by U and V as well. Next, in Section 3.2, the proposed algorithm is described. of simplexes with disjoint non empty interiors, so that K = m j=1 S j . The set of vertices of the partition P is denoted by V(P).

Note that K is then the convex hull generated by V(P).

The algorithm is based on Monte Carlo simulations, combined with a sequence of approximations on compact sets R 0 , . . . , R n-1 , determined by partitions P 0 , . . . , P n-1 .

The idea behind the algorithm is quite simple. Given approximations Ṽk , . . . , Ṽn of V k , . . . , V n , one first get Vk-1 by estimating V k-1 at every vertices x ∈ V(P k-1 ), using Monte Carlo simulations, and then, one uses a linear interpolation of Vk-1 , to define Ṽk-1 at any point x ∈ R k-1 . Then Ṽk-1 may be extended to all of X using projections.

More precisely, one may proceed through the following steps.

3.2.1. Algorithm.

C1: Set Ṽn = Ũn = fn , where fn = f n on R n and for any x ∈ R c n , fn (x) = fn (x ), where

x is the unique point in R n so that d(x, R n ) = d(x, x ); C2: For k = n, . . . , 1, generate ξ 1,k , . . . , ξ N k ,k according to μ k ; C3: For every x ∈ V(P k-1 ), define Vk-1 (x) = 1 N k N k i=1 Ũk {π k (x, ξ i,k )} ,
where Ũk = max fk , Ṽk , and where fk = f k on R k and for any x ∈ R c k , fk (x) = fk (x ), where x is the unique point in R k so that d(x, R n k) = d(x, x ); C4: Interpolate linearly Vk-1 over R k-1 , as in Definition C.1, and call it Ṽk-1 .

For the precise meaning of linear interpolation and ways to implement it, see Section C. Remark 3.1. First, contrary to most Monte Carlo algorithms, one does not simulate trajectories of the price process S. This is the key to preserve monotonicity and convexity, as illustrated in Figure 1.Next, in the case where μ 1 = . . . , μ n = μ and N 1 = . . . = N n = N , one could take the same set of random points ξ 1 , . . . , ξ N from distribution μ, thus reducing calculations. Last but not least, note that since f k is not interpolated, which would induce extraneous errors, our algorithm cannot be interpreted as generating a tree-based algorithm.

As a by-product of the method, one obtains an estimation E of the exercise region E. Such approximations are illustrated in Figures 4 and6. The following results show that the algorithm produces good uniform approximations of the Snell envelope U . Theorem 3.5. Suppose that for all k ∈ {0, . . . , n}, f k , V k and U k are continuous. Let R 0 be a given compact convex subset of X. Let > 0 be given. Then one can find compact convex sets R 1 , . . . , R n-1 ⊂ X, partitions P 0 , . . . P n-1 generating respectively R 0 , . . . , R n-1 , such that almost surely, one can find integers N 10 , . . . , N n0 , so that for the simple interpolation method,

sup x∈R k |U k (x) -Ũk (x)| = U k -Ũk R k ≤ V k -Ṽk R k < , whenever N 1 ≥ N 10 , . . . , N n ≥ N n0 .
The next result shows that using convex interpolations, as in Definition C.2, one can preserve the convexity of the Snell envelope. The proofs of these results are given is Appendix B.

Implementation issues

The first example of application is the classical Black-Scholes-Merton setting. Examples of calculations are provided for the case of an American call-on-max option on one and two assets. The second example is the N-GARCH model, as studied in [START_REF] Duan | American option pricing under GARCH by a Markov chain approximation[END_REF] and [START_REF] Duan | Approximating American option prices in the GARCH framework[END_REF]. Examples of calculations are given for the American put option. In each numerical example, estimated lower and upper bounds are calculated, using the modified Andersen-Broadie algorithm described in Section 2.5. 4.1. Geometric Brownian motion. Suppose that under the risk neutral measure Q, the model satisfies (4.1)

dX i (t) = (r -δ i )X i (t)dt + σ i X i (t)dW i (t),
where W 1 , . . . , W d are dependent Brownian motions with correlation matric ρ. Setting

α i = r -δ i -σ 2 i /2 Δ and β i = σ i Δ 1/2 , then the conditional law of X t given X t-1 = x is
log-normal, with mean vector log(x) + α, and covariance matrix C, where log(x) stands for the vector with components log(x i ), 1 ≤ i ≤ d, and

C ij = β i β j ρ ij , 1 ≤ i, j ≤ d. It follows that Assumption 1 is satisfied with Y k ∼ N (α, C) and π (j) k (x, y) = x j e y j , j = 1, . . . , d.
Furthermore, for each fixed y = (y 1 , . . . , y d ), π k (•, y) is continuous, non decreasing, convex and concave. Since x → max 1≤j≤d x j is non decreasing and convex and x → min 1≤j≤d x j is non decreasing and concave, it follows that for the call-on-max, U k and V k are convex and non decreasing, while for the put-on-min, U k and V k are convex and non increasing.

4.1.1. Numerical illustration for the American call option. For this example, S t satisfies (4.1) with the same parameters as in [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], i.e., r = 0.05, δ = 0.1, σ = 0.2, T = 3 and K = 100. They did their calculations assuming that S 0 = 100, that is the option is at-the-money. For this reason, one could take

R 0 = • • • = R n = R = [23, 230],
so that P (S T ∈ R|S 0 = 100) ≥ 0.997. The Monte Carlo algorithm described in Section 3.2 was implemented using M = 1000 equally spaced vertices in R. Also, one used

N 10 = • • • = N n0 = N = 10000 new trajectories each time.
As argued previously, for each k = 0, . . . , n, U k and V k are convex and non-decreasing.

Since d = 1, one can take take the convex interpolation described in Example 1 of Appendix C. It follows from Corollary 3.6 that Ũk and Ṽk are also convex and non decreasing.

Table 1 reports the results of the simulations for n = 2 and n = 10 exercise periods. In each case, the primal-dual algorithm was implemented with N 0 = 2000000, N 1 = 1500 , N 2 = N 3 = 10000. The corresponding estimated values of Ũ0 are displayed in Figure 2. As argued previously, for each k = 0, . . . , n, U k and V k are convex and non-decreasing. However, due to computation time constraints, the interpolation method used was only (locally) convex in each of the 199 2 sub-rectangles, using the the interpolation described in Example 2 of Appendix C.1. Therefore the estimated functions Ũk and Ṽk are not necessarily convex over R. However, they are non decreasing, using Corollary 3.6.

Table 2 reports the results of the simulations for n = 9 exercise periods, as in [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], except that they only consider initial values S 0 = (90, 90), S 0 = (100, 100) and S 0 = (110, 110). In each case, the primal-dual algorithm was implemented with the same number of iterations, i.e., N 0 = 2000000, N 1 = 1500 , N 2 = N 3 = 10000.

The results are quite similar. However, one obtains as a bonus the estimation of U k over R = [25,230] 2 . This is illustrated in Figure 3, where the values along the diagonal are also displayed. One can remark that around the diagonal, the Û0 is not convex, due to the fact that the convex interpolation was not used.

It is quite interesting to plot the estimated exercise region at time 0 given by the set of all points s = (s 1 , s 2 ) ∈ [25, 230] 2 such that Ṽ0 (s) ≤ f 0 (s) = max{0, max(s 1 , s 2 ) -100)}. This is done in Figure 4. It is quite surprising to see that along the diagonal (s, s), the option is not exercised unless the s > 220. The values S 0 = (160, 160) and S 0 = (200,200) were included in Table 2 to show that the option is not exercised at time 0 for these values, even if they are quite large compared to the strike price K = 100. Looking at the confidence intervals, one sees that 60 < 69.4746 and 100 < 107.3452.

S k = S k-1 e r-δ-h k /2+h 1/2 k Y k h k = β 0 + β 1 h k-1 + β 2 h k-1 (Y k-1 -θ) 2 ,
where Y k ∼ N (0, 1). Then {S k } k≥0 or {(S k , h k )} k≥0 are not Markov processes. However, it is easy to check that {X k = (S k , h k+1 )} k≥0 is a Markov chain. In fact, Assumption 1 holds

with π k (s, h, y) = se r-δ-h/2+h 1/2 y , β 0 + β 1 h + β 2 h(y -θ) 2 .
It follows that for any given y ∈ R, π k (s, h, y) is monotone in s but it is not monotone in h, nor it is convex or concave. It is well-known, e.g., [START_REF] Duan | American option pricing under GARCH by a Markov chain approximation[END_REF], that under the condition

β 1 +β 2 (1+θ 2 ) < 1, X k has a stationary distribution and lim k→∞ E(H k ) → h = β 0 {1 -β 1 -β 2 (1 + θ 2 )} -1 .
In [START_REF] Duan | American option pricing under GARCH by a Markov chain approximation[END_REF] and [START_REF] Duan | Approximating American option prices in the GARCH framework[END_REF], the authors proposed two kinds of approximations for the American put option on S k . They calculated the option prices for Table 3. Estimated values of the American put options at time 0 for n = 10 and n = 30 exercise periods, with initial asset value S 0 = K = 50, three initial volatilities h 1 and strike price K = 50.

Conclusion

Monotonicity and convexity are important properties shared by many American option prices, especially when it comes to define the exercise region. In this paper, it was shown that these properties can be preserved by a Monte Carlo algorithm which is easy to implement and yields the whole Snell envelope. Also a modified approach inspired by the work of [START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF] can be used to improve the precision of the estimation. Numerical results showed that the proposed methodology provides accurate results

for American option prices when assets follows a Markov chain with continuous state space.

In addition, the Monte Carlo algorithm can also be applied for finding optimal hedging strategies, as in [START_REF] Hocquard | Optimal hedging strategies with an application to hedge fund replication[END_REF] and [START_REF] Papageorgiou | Replicating the properties of hedge fund returns[END_REF].

g(φ(x, y)) is continuous and non decreasing, and (ii) follows from the Bounded Convergence Theorem.

To prove (iii), assume that g is convex and non decreasing, and that φ(•, y) is convex, for any fixed y ∈ E. Then, for any λ ∈ [0, 1], for any y ∈ E, and for any x, x ∈ X,

φ(λx + (1 -λ)x , y) ≤ λφ(x, y) + (1 -λ)φ(x , y), so g(φ(λx + (1 -λ)x , y)) ≤ g(λφ(x, y) + (1 -λ)φ(x , y)) ≤ λg(φ(x, y)) + (1 -λ)g(φ(x , y)).
Integrating with respect to y, one obtains Mg(λx

+ (1 -λ)x ) ≤ λM g(x) + (1 -λ)Mg(x ),
proving that Mg is convex. Finally, if g is convex and non increasing, and φ(•, y) is concave for any fixed y, then, for any λ ∈ [0, 1], and for any y, φ(λx

+ (1 -λ)x , y) ≥ λφ(x, y) + (1 - λ)φ(x , y), so g(φ(λx + (1 -λ)x , y)) ≤ g(λφ(x, y) + (1 -λ)φ(x , y)) ≤ λg(φ(x, y)) + (1 -λ)g(φ(x , y)).
Integrating with respect to y, one obtains Mg(λx

+ (1 -λ)x ) ≤ λM g(x) + (1 -λ)Mg(x ).
This completes the proof.

The next result is easy to prove.

Proposition A.2. For any real numbers x, y, z, w,

|x ∨ z -y ∨ w| ≤ |x -y| ∨ |z -w|.
In particular, |x ∨ zx ∨ w| ≤ |z -w|.

Appendix B. Proofs of the main results

B.1. Proof of Proposition 3.2. First, for any k ∈ {0, . . . , n}, is also continuous.

U k (x) = sup τ ∈T k,n E {f τ (X τ )|X k = x} , so it follows that V k (x) ≤ U k (x) ≤ H k (x) = n j=k E {f j (X j )|X k = x}. Next,
So suppose that V k is continuous. Then so is U k . Suppose that (x j ) j≥1 is a sequence converging to x ∈ X. Then, U k (π k (x j , y)) ≤ H k (π k (x j , y)), and for all y ∈ Y,

U k (π k (x j , y)) → U k (π k (x, y)), H k (π k (x j , y)) → H k (π k (x, y))
as j → ∞. Moreover, by hypothesis,

H k (π k (x j , y))μ k (dy) = H k-1 (x j ) -f k-1 (x j ) j→∞ -→ H k-1 (x) -f k-1 (x) = H k (π k (x, y))μ k (dy).
Thus, by Fatou's Lemma, V k-1 (x j ) → V k-1 (x) as j → ∞. Therefore V k-1 and U k-1 are continuous. Hence the result.

B.2. Proof of Theorem 3.5.

Proof. Recall that for all k ∈ {0, . . . , n}, the f k , V k and U k are continuous. Set

M k g(x) = Y g(π k (x, y))μ k (dy) and Vk-1 (x) = Mk Ũk (x) = 1 N k N k i=1 Ũk {π k (x, ξ i,k )} .
The approximation algorithm can be summarized as follows:

Ũn = f n and ⎧ ⎨ ⎩ Vk-1 = Mk Ũk on V k-1 = V(P k-1 ), Ṽk-1 = I P k-1 Vk-1 Ũk-1 = f k-1 ∨ Ṽk-1 . k = n, . . . , 1,
where I P g denotes the linear interpolation of g over the compact set R with associated partition P of R as in Definition C.1. For simplicity set I k-1 = I P k-1 and δ k-1 = mesh(P k-1 ).

First, for any compact subset R of X, Ũk -

U k R ≤ Ṽk -V k R , 0 ≤ k ≤ n, by Proposition
A.2. Therefore, to prove the result, it suffices to show that for any k = 1, . . . , n, given a compact set R k-1 and > 0, one can find δ k-1 > 0, a partition P k-1 with mesh(P k-1 ) < δ k-1 , another compact set R k and an integer N k0 so that Ṽk-1

-V k-1 R k-1 ≤ Ṽk -V k R k + .
To begin, one has, for any k = 1, . . . , n,

V k-1 -Ṽk-1 = (V k-1 -I k-1 V k-1 ) + I k-1 V k-1 -I k-1 Vk-1 . Since I k-1 V k-1 -I k-1 Vk-1 R k-1 = sup x∈V k-1 V k-1 (x) -Vk-1 (x) = V k-1 -Vk-1 V k-1 = M k U k -Mk Ũk V k-1
and

M k U k -Mk Ũk V k-1 ≤ M k U k -Mk U k V k-1 + Mk U k -Mk Ũk V k-1 ≤ M k U k -Mk U k V k-1 + V k -Ṽk R k + Mk (U k 1 R c k ) V k-1 + Ũk R k Mk 1 R c k V k-1 , it follows that Ṽk-1 -V k-1 R k-1 ≤ ω(V k-1 , R k-1 , δ k-1 ) + M k U k -Mk U k V k-1 + V k -Ṽk R k + Mk (U k 1 R c k ) V k-1 + Ũk R k Mk 1 R c k V k-1
, using the above inequalities together with Remark C.1.

The aim is to show that if N 1 , . . . , N n are large enough, and if R 1 , . . . , R n are large enough, and if δ 0 , δ n-1 are small enough, then the terms ω(

V k-1 , R k-1 , δ k-1 ), M k U k -Mk U k V k-1 , Mk (U k 1 R c k ) V k-1 and Ũk R k Mk 1 R c k V k-1
can be arbitrarily small, for all 1 ≤ k ≤ n.

The first term is easy to handle since V k-1 is continuous and R k-1 is compact. To handle the second term, just use the strong law of large numbers, since V k-1 is finite. The last two terms can be made arbitrarily small since, by the strong law of large numbers,

Mk (U k 1 R c k ) V k-1 and Mk 1 R c k V k-1 converge almost surely to M k (U k 1 R c k ) V k-1 and to M k 1 R c k V k-1
, which both can be as small as one wants, choosing R k large enough. over P is the (unique) function g defined in the following way:

• If x ∈ S ⊂ R, where S ∈ P is a simplex with vertices x 1 , . . . , x d+1 , g(x) = d+1 i=1 λ i g(x i ),
where the barycenters {λ 1 , . . . , λ d+1 } are the unique solution of

x = d+1 i=1 λ i x i , d+1 i=1 λ i = 1, λ i ∈ [0, 1], i = 1, . . . d + 1.
• If x ∈ R, let x R be the (unique) closest point to x that belongs to R, and set g(x) = g(x R ).

Uniqueness follows from the convexity of R and the strict convexity of the Euclidean norm.

Remark C.1. Note that since each x i is extreme in S, the unique solution of

x i = d+1 j=1 λ j x j , d+1 j=1 λ j = 1, λ j ∈ [0, 1], j = 1, . . . d + 1,
is λ i = 1 and λ j = 0 for all j = i, yielding g(x i ) = g(x i ) for all 1 ≤ i ≤ m. Moreover, g is affine on each simplex, justifying the term "linear interpolation".

= {a = x 0 < x 1 < • • • < x m = b}, i.e the simplexes are [x i-1 , x i ], i = 1, . . . , m. Set Δ i = g(x i )-g(x i-1 ) x i -x i-1 , 1 ≤ i ≤ m.
Then the linear interpolation of g over R is given by

g(x) = ⎧ ⎨ ⎩ g(a), x ≤ a, g(x i ) + (x -x i )Δ i+1 , x ∈ [x i , x i+1 ], i = 0, . . . , m -1, g(b)
x ≥ b.

If g is monotone, the slopes Δ i all have the same sign, so g has the same monotonicity. If g is convex, the slopes Δ i are non decreasing, so g is also convex.

Definition C.2. Given a convex function g and a partition P of R, a convex linear interpolation of g over P is a function g defined on X such that (CLI1) g(x) is convex on R;

(CLI2) g(x) = g(x), for all x ∈ V(P);

(CLI3) For any x ∈ R, there exists non negative numbers λ z , z ∈ V(P), such that The following result proves that convex linear interpolation exists and it can be obtained via the simplex algorithm.

Lemma C.3. Suppose C is the convex hull of {x 1 , . . . , x n }. Suppose also that g 1 , . . . , g n are real numbers. For any x ∈ C, set

Λ x = λ ∈ [0, 1] n ; n i=1 λ i x i = x, n i=1 λ i = 1 . Define g(x) = inf λ∈Λx n i=1 λ i g i .
Then g is convex over C and g(x i ) ≤ g i . If in addition, g is a convex function over C such that g(x i ) = g i , then g is a linear convex interpolation of g. In particular, g(x) ≤ g(x) for all x ∈ C, and g(x i ) = g i , for all 1 ≤ i ≤ n.

Suppose that the points x 1 , . . . , x n form a grid. Then if g is monotone, then g has the same monotonicity, i.e., if g is non-decreasing in x j the so is g, for any 1 ≤ j ≤ d.

Proof. First, it is obvious that f (x i ) ≤ g i for any 1 ≤ i ≤ n. Next, suppose that a, b ∈ C.

Take λ ∈ Λ a and λ ∈ Λ b . Let t ∈ [0, 1] be given and set

x = ta + (1 -t)b. Then tλ + (1 -t)λ ∈ Λ x , so g(x) ≤ n i=1 tλ i + (1 -t)λ i g i = t n i=1 λ i g i + (1 -t) n i=1 λ i g i .
Taking the infimum over all λ ∈ Λ a and all λ ∈ Λ b in the last inequality yields g

(x) ≤ tg(a) + (1 -t)g(b). Hence g is convex over C.
Next, suppose that g is convex over C, and g(x i ) = g i for all i = 1, . . . , n. Let x ∈ C be given. Since λ : Λ x → n i=1 λ i g i is continuous and Λ x is compact, there exists λ ∈ Λ x so that g

(x) = n i=1 λ i g i . Therefore g(x) ≤ n i=1 λ i g(x i ) = n i=1 λ i g i = g(x). Hence g(x) ≤ g(x).
In particular,

g i = g(x i ) ≤ g(x i ) ≤ g i , showing that g(x i ) = g i , for all 1 ≤ i ≤ n.
Finally, to prove the monotonicity property of g, note that it suffices to prove it for the interior points of C. For simplicity, the proof is given for a two-dimensional grid but it holds true for any dimension. Let the grid be given by G = (y i , z j ); 1 ≤ i ≤ M 1 , 1 ≤ j ≤ M 2 ) and set g ij = g(y j , z j ). From linear programming, it is known that for any x = (y, z), then g(y, z) = sup a,b,c∈A a + by + cz,

where A = {a, b, c ∈ R; a + by i + cz j ≤ g ij , for all (y i , z j ) ∈ G}.

It is known that there are 3 points on the grid where the constraints are in fact an equality. One of these points, say (y i , z j ) must an interior point of the grid, for otherwise (y, z) could not be an interior point of C. Therefore a + by i + cz j = g ij . One also has g ij + b(y i+1y i ) = a + by i+1 + cz j ≤ g i+1,j , so b ≤ g i+1,j -g ij y i+1 -y i . Similarly, g ijb(y iy i-1 ) = a + by i-1 + cz j ≤ g i-1,j , so b ≥ g ij -g i-1,j y i -y i-1 . It follows that if g is non-decreasing in y, b ≥ 0.

Hence, for > 0 small enough, g(y + , z)g(y, z) = b ≥ 0. The same argument applies in g is non-increasing in y. Hence the result. Example 2. Suppose that one wants to interpolate linearly a convex function f over [0, 1] 2 , knowing only the values at extremes, that is f (0, 0), f (1, 0), f (0, 1) and f (1, 1). One could take the partition P = {C 1 , C 2 }, where C 1 = {(x, y) ∈ [0, 1] 2 ; y ≤ x} and C 2 = {(x, y) ∈ [0, 1] 2 ; x ≤ y}. This would lead to the quick interpolation described above, namely g(x, y) = f 1 (x, y), (x, y) ∈ C 1 f 2 (x, y), (x, y) ∈ C 2 , where f 1 (x, y) = x{f (1, 0)-f (0, 0)}+y{f (1, 1)-f (1, 0)}+f (0, 0) and f 2 (x, y) = x{f (1, 1)f (0, 1)} + y{f (0, 1)f (0, 0)} + f (0, 0).

On the other hand, one could also consider the partition P = {C 3 , C 4 }, where C 3 = {(x, y) ∈ [0, 1] 2 ; x + y ≤ 1} and C 4 = {(x, y) ∈ [0, 1] 2 ; x + y ≥ 1}. That would lead to the following interpolation g of f : g (x, y) = f 3 (x, y), (x, y) ∈ C 3 f 4 (x, y), (x, y) ∈ C 4 , where f 3 (x, y) =

x{f (1, 0)-f (0, 0)}+y{f (0, 1)-f (0, 0)}+f (0, 0), f 4 (x, y) = x{f (1, 1)-f (0, 1)}+y{f (1, 1)f (1, 0)} -Δ + f (0, 0), and Δ = f (1, 1)f (0, 1)f (1, 0) + f (0, 0).

So which partition should be chosen? If convexity is at stake, the answer relies only on the sign of Δ. In fact, it it easy to check that g is convex on [0, 1] 2 if and only if Δ ≤ 0, while g is convex on [0, 1] 2 if and only if Δ ≥ 0. It is also easy to check that Lemma C.3 applied with the four extreme points (0, 0), (1, 0), (0, 1) and (1, 1) of [0, 1] 2 yields g when Δ ≤ 0 (resp. g when Δ ≥ 0).

  American options. Algorithms for evaluating American options can be divided in three classes. (i) Tree-based methods: For these algorithms, values of U are estimated at a finite number of points, possibly random. In most cases, the underlying process X is replaced by a finite state Markov chain Xk k≥0 . The main representatives are the algorithms based on the binomial and multinomial methods, initiated by Cox

Figure 1 .

 1 Figure 1. Full Monte Carlo method vs suggested Monte Carlo method with N=100

Corollary 3. 6 .

 6 Suppose in addition that for each k, f k is monotone and π k (•, y) is non decreasing for any fixed y ∈ Y. If the quick interpolation method is used, then Theorem 3.5 holds true and the functions Ṽk and Ũk have the same monotonicity as the U k 's and V k 's when restricted to R k . If a convex interpolation on a grid is used instead of a linear interpolation, then the conclusion of Theorem 3.5 holds true and the functions Ṽk and Ũk have the same monotonicity as the U k 's and V k 's when restricted to R k . If in addition f k is convex and non decreasing and π k (•, y) is convex and non decreasing for any fixed y, then the U k 's and V k 's are convex when restricted to R k . The same holds true if f k is convex and non increasing and π k (•, y) is concave and non decreasing, for any fixed y.

  with the same parameters as in[START_REF] Andersen | Primal-dual simulation algorithm for pricing multidimensional American options[END_REF], i.e., r = 0.05, δ = 0.1, σ = 0.2, ρ = I 2 , T = 3 and K = 100. They did their calculations assuming that S , that is the option is at-the-money.

Figure 2 .

 2 Figure 2. Estimated value of U 0 (s) for an American call option for s ∈ [25, 230], with n = 2 exercise periods (panel a) and n = 10 exercise periods (panel b).

Figure 3 .Figure 4 .

 34 Figure 3. Estimated value of a call-on-max American option at period 0 for s ∈ [25, 230] 2 (panel a) and along its diagonal (panel b), with n = 9 exercise periods.

S 0 =Figure 5 .

 05 Figure 5. Estimated value of U 0 for an American put option based on a N-GARCH model, with 10 exercise periods (panel a) and 30 exercise periods (panel b).

Figure 6 .

 6 Figure 6. Estimated exercice region E at period 0 for the American put option based on a N-GARCH model, for 10 exercise periods (panel a) and for 30 exercise periods (panel b).

B. 3 .

 3 Proof of Corollary 3.6. It follows from that Vk-1 = Mk Ũk is continuous, monotone and convex, if Ũk is, by Proposition A.1. As shown in Proposition C.4 and Lemma C.3, the quick interpolation and the convex interpolation on a grid both preserve all monotonicity. In addition, the convex interpolation preserves convexity, by Lemma C.3. Finally, Theorem 3.5 yields the almost sure uniform convergence. Appendix C. Linear interpolations Definition C.1. Given a function g and a partition P of R, the linear interpolation of g

Finally, g is

  continuous and bounded on X andsup x∈R |g(x)g(x)| ≤ ω(g, R, mesh(P)),where mesh(P) = max S∈P sup x,z∈S xz and ω(g, R, δ) is the modulus of continuity of g over R, i.e., ω(g, R, δ) = sup x,z∈R, x-z ≤δ |g(x)g(z)|.Example 1. Suppose d = 1. Then the linear interpolation g of a monotone (respectively convex) function g on R = [a, b] is monotone (respectively convex). To see that, let P be the partition given by P

  For any x ∈ R, g(x) = g(y), where y is the unique point satisfying d(x, R) = d(x, y).

C. 1 .

 1 Quick linear interpolation on rectangles. Suppose that g is defined on R = [0, 1] d , and g is known at all extreme points of R. If one wants to linearly interpolate g, as in Definition C.1, a convenient choice for the partition P is the set of all d! simplexes defined, for any permutation π of {1, . . . , d}, byS π = {y ∈ [0, 1] d ; y π 1 ≤ • • • , y π d }. Note that the vertices of S = {y ∈ [0, 1] d ; y 1 ≤ • • • , y d } are 0 and the vectors u 1 , . . . , u d are such that the A = (u 1 | • • • |u d ) is defined by A ij = 1 if 1 ≤ j ≤ i ≤ d, and A ij = 0 otherwise. Any y ∈ S can be uniquely written as y = d j=1 λ j u j , with λ 1 = y 1 , λ k = y ky k-1 , 2 ≤ k ≤ d,andλ d+1 = 1y d .To define g at x = (x 1 , . . . , x d ) ∈ R, proceed in the following way:• Let π be a permutation such that x π 1 ≤ x π 2 ≤ • • • ≤ x π d and let P be the associated permutation matrix. Set y = P x;• Set λ 1 = y 1 , λ k = y ky k-1 , 2 ≤ k ≤ d, and λ d+1 = 1y d ; • Set ĝ(x) = d j=1 λ j g P u j + λ d+1 g(0). Remark C.2. If R = d j=1 [a j , b j ], then set x j =x j -a j b j -a j , 1 ≤ j ≤ d and apply the last procedure. The extreme points necessary for the interpolation are then given by a+DP u j , 1 ≤ j ≤ d and a, where D is the diagonal matrix with diagonal b-a. The quick interpolation on a grid is just the repeated application of the quick interpolation on the partition of rectangles generating the grid.

Remark 2.1. In

  general, if A k consists of all possible linear combinations of a finite number d k of basis functions, then Vk can be found using regression, and the computations can be quite fast. Unfortunately, if d k is small compared to N , then the approximation Vk of V k can be poor. However if V k is convex or monotone, Vk will not inherit the same properties.

	2.4.

Longstaff-Schwartz algorithm. Using

approximating classes A 0 , . . ., A n-1 , estimate E by E, where E

  When the monotonicity properties do not hold for π k , one needs another assumption to get the continuity of the Snell envelope.The next proposition deals with convexity of the Snell envelope.

	Proposition 3.2. Under Assumptions 2 and 3, V k and U k are continuous for all k ∈
	{0, . . . , n}.

3.1. Properties of U and V . The first two propositions are about monotonicity and continuity. Proposition 3.1. Assume that for every k ∈ {0, . . . , n}, f k is non decreasing (respectively non increasing) and π k (•, y) is non decreasing for any y ∈ Y. Then, for every k ∈ {0, . . . , n}, U k and V k are non decreasing (respectively non increasing). If in addition, f k is continuous for any k ∈ {0, . . . , n}, then U k and V k are also continuous for all k ∈ {0, . . . , n}. Proposition 3.3. Assume that for every k ∈ {0, . . . , n}, f k is convex and non decreasing, and π k (•, y) is convex and non decreasing for any y ∈ Y. Then U k and V k are convex and non decreasing, for all k ∈ {0, . . . , n}. If f k is convex and non increasing, and if π k (•, y) is concave and non decreasing for any y ∈ Y and any k ∈ {0, . . . , n}, then U k and V k are also convex and non increasing, for all k ∈ {0, . . . , n}.

The proofs of Propositions 3.1 and 3.3 follow readily from Proposition A.1, since convexity, continuity and monotonicity are preserved by taking the max, that is, if g and h are convex and non decreasing (respectively non increasing), then max(g, h) = g ∨ h is also convex and non decreasing (respectively non increasing). Clearly, if g and h are continuous, then g ∨ h is continuous. Proposition 3.2 is proven in Section B.1 3.2. Description of the algorithm and justification. Before describing the algorithm, one has to define what one means by a partition. Definition 3.4. A partition P of a compact convex set K, is any finite set P = {S 1 , . . . , S m }

Table 1 .

 1 Estimated values of the call option with initial value S 0 = 100 and strike price K = 100.

	Method	n = 2	n = 10
	Monte Carlo	7.027575	7.950981
	Lower bound	7.179106 (0.009548) 7.980827 (0.008106)
	Upper bound	7.173063 (0.004365) 7.973237 (0.007579)
	95% Confidence interval	[7.1604, 7.1816]	[7.9650, 7.9881]
	Binomial tree	7.18	7.98
	4.1.2. Numerical illustration for the American call-on-max option on two assets. For this
	example, S t = S	(1) t , S

Table 2 .

 2 Estimated values of the call-on-max option with several initial values S 0 and strike price K = 100.

	S 0	Method	Ũ0 (S 0 , S 0 )
	90	Monte Carlo	8.095579
		Lower bound	8.068165 (0.008518)
		Upper bound	8.070320 (0.008704)
		95% Confidence Interval	[8.0515, 8.0873]
		Binomial tree	8.075
	100	Monte Carlo	13.900166
		Lower bound	13.912179 (0.010687)
		Upper bound	13.904105 (0.009686)
		95% Confidence Interval	[13.8913, 13.9230]
		Binomial tree	13.902
	110	Monte Carlo	21.300620
		Lower bound	21.334454 (0.012411)
		Upper bound	21.338053 (0.011061)
		95% Confidence Interval	[21.3102, 21.3597]
		Binomial tree	21.345
	160	Monte Carlo	69.244785
		Lower bound	69.507220 (0.016665)
		Upper bound	69.519498 (0.011474)
		95% Confidence Interval	[69.4746, 69.5419]
		Binomial tree
	200	Monte Carlo	106.990926
		Lower bound	107.372955 (0.014175)
		Upper bound	107.370157 (0.008543)
		95% Confidence Interval [107.3452, 107.3869]
		Binomial tree

  by hypothesis,V n and H k (π k (•, y)) are continuous for any y ∈ Y and any k ∈ {0, . . . , n}. To show that V k is continuous for any k ∈ {0, . . . , n}, it suffices to show that if V k is continuous, then V k-1

If true, their argument would lead to the following (false) result: if ( k ) n k=1 are martingale differences and if τ is a random time with values in {1, . . . , n}, τ not necessarily a stopping time, then E( τ ) = 0.

Funding in partial support of this work was provided by the Natural Sciences and Engineering Research Council of Canada and the Fonds québécois de la recherche sur la nature et les technologies.

1

Suppose that g ∈ B + (X), that is g : X → [0, ∞) is measurable. Suppose also that φ : X × Y → X is measurable and that for any y ∈ Y , φ(•, y) is continuous. Moreover, assume that for any x ∈ X, g(φ(x, •)) is integrable with respect to a measure ν on Y. Then the function Mg, defined on X by x → Mg(x) = Y g(φ(x, y))ν(dy), belongs to B + (X).

then Mg is also convex and non decreasing.

(iv) If g is convex and non increasing, and if, for any fixed y ∈ Y, φ(•, y), is concave, then Mg is also convex and non increasing.

Proof. Suppose g is non decreasing for any fixed y. Therefore, for any fixed y, x → g(φ(x, y)) is non decreasing. If g is non increasing, then for any fixed y, x → g(φ(x, y)) is non increasing. Integrating with respect to ν yields (i). If, for any fixed y, g is non decreasing and continuous, then x → g(φ(x, y)) is continuous and non decreasing. Hence (ii) follows from the Dominated Convergence Theorem. If g is non increasing and continuous then x → Proposition C.4. Suppose that g is monotone and g represents the quick interpolation of g over a grid. Then g has the same monotonicity as g, i.e., if g is non-decreasing in x j the so is g, for any 1 ≤ j ≤ d.

Proof. It suffices to prove the result on any rectangle forming the partition of the grid.

By translation and rescaling, one can consider that R = [0, 1] d . In that case, note that one can write g(x) = d j=1 y j g(P u j )g(P u j+1 ) , where u d+1 = 0. Next, P u j -P u j+1 = P e j , where e j is the unit vector with (e j ) j = 1 and (e j ) i = 0, i = j. Hence P u j ≥ P u j+1 , so if g is non-decreasing in x j , then b j = g(P u j )g(P u j+1 ) ≥ 0. It follows that if x in an interior point of one of the simplexes, then x + e j belongs to the same simplex if > 0 is small enough. Suppose that the rank of x j is k. Then, the rank of x j + is also k and g(x + e j )g(x) = g(P u k )g(P u k+1 ) ≥ 0, since P e k = e j .

Hence the monotonicity in x j is preserved.

INRIA Bordeaux-Sud Ouest and Université Bordeaux I, HEC Montréal, Université Nice-Sophia-Antipolis E-mail address: bruno.remillard@hec.ca