Tristan Guillaume

Discretely monitored lookback and barrier options : a semianalytical approach Revised version, 2006

Keywords: option, lookback option, barrier option, discrete monitoring, numerical integration, dimension. JEL classification : G13

All the explicit formulae for the valuation of lookback and barrier options available in the financial literature assume continuous monitoring of the underlying asset. In practice, however, monitoring is always discrete, and the gap between continuously and discretely monitored option values can be very large. In this paper, we provide explicit formulae for discretely monitored lookback and barrier options. They allow for non-constant volatility, interest rate, dividend rate and barrier parameters that vary as step functions of time. They can deal with any number and spacing of monitoring dates. They are not restricted to particular payoffs or strike price specifications. We also provide a simple rule for the numerical integration of these high-dimensional formulae, as well as an efficient interpolation method.

Introduction

Lookback and barrier options are among the most heavily traded OTC derivatives, particularly in the foreign exchange markets. They are also embedded in a lot of popular structured products in equity and interest rate markets. Lookback options allow investors to « buy at the lowest » and « sell at the highest », but at the cost of a substantially increased premium. Barrier options owe their success to their low price and their strong leverage effect, as well as to the precision and the flexibility with which they can adapt to the needs or views of market participants.

In a stylised Black-Scholes framework, closed form formulae can be obtained for the value of these contracts. They assume continuous monitoring of the underlying asset. In the real markets, however, monitoring is always discrete, for practical but also financial reasons : discretely monitored lookback options are more affordable, thus contributing to eliminate the main obstacle to their commercial success ; discretely monitored knock-out barrier options bear a diminished risk of ending worthless ; and discretely monitored knock-in barrier options benefit from enhanced leverage. Not to mention, in the trader's perspective, the case for facilitated hedging.

As the pricing bias caused by the assumption of continuous monitoring can be very large, alternative approaches are needed to tackle discrete monitoring. Monte Carlo simulation, which is rather slow and inaccurate when it comes to valuing path-dependent contracts with continuous monitoring, is a better choice when extrema are discretely monitored. Nevertheless, the use of advanced variance reduction techniques is uneasy. In particular, closed form lookback and barrier option formulae that might be used as control variates are available only for simple contract specifications; and even when they do exist, they can hardly be used if the number of fixing dates is moderate or low. Likewise, conditional Monte Carlo cannot be easily implemented to reduce the number of random number samplings at each simulation. Lattice-based methods, using either binomial trees [START_REF] Cheuk | Currency lookback options and observation frequency : a binomial approach[END_REF] or, better, trinomial trees [START_REF] Cheuk | Complex barrier options[END_REF][START_REF] Ahn | Pricing discrete barrier options with an adaptive mesh[END_REF]) can be accurate and fast. However, they entail stability and convergence issues, especially when parameters and barrier levels are time-varying, and they do not easily cope with complex payoffs. [START_REF] Broadie | A continuity correction for discrete barrier options[END_REF]1999) propose an approximation formula that involves a correction to the continuous-monitoring closed form formulae. This approach is simple to use and efficient. But, as mentioned earlier, continuous-monitoring closed form formulae are scarce. Furthermore, the accuracy of this approach deteriorates as the number of fixing dates decreases or as the fixing dates become more and more unevenly spaced (a limitation shared by the extrapolation method developed by [START_REF] Levy | Discrete by nature[END_REF]. Sullivan (2000) introduces an efficient method combining Gaussian quadrature and Chebyshev polynomial approximation, but his analysis covers only the most simple barrier option payoffs and assumes all parameters are constant. Finally, Öhgren (2001) develops an interesting approach based on a result known as the Spitzer identity but, as the author points out himself, it is too restrictive inasmuch as it can be used only if the option is in-the-money and for specific forms of lookback payoff.

In this paper, we provide explicit formulae for discretely monitored European lookback and barrier options that aim to avoid the above limitations :

-their accuracy is not contingent on the spacing of monitoring dates -they are not restricted to particular payoffs or strike price specifications -they allow for non-constant volatility, interest rate, dividend rate and barrier parameters that vary as step functions of time, without loss of convergence or stability Besides, those formulae can show us how the various parameters affect the solution and bring out the interrelationships among them by mere differentiation. However, they have to be numerically integrated, in dimensions that may look daunting as the number of monitoring dates increases. We address this issue by providing simple numerical integration schemes. Furthermore, we show that, as long as monitoring dates are evenly spaced, option values with an arbitrarily large number of monitoring dates can be easily interpolated by means of cubic and quintic splines within our analytical framework, which drastically cuts computational time.

This paper is organized as follows. Section 2 derives explicit formulae for discretely monitored European lookback options. Section 3 derives explicit analytical representations for discretely monitored European barrier options in the form of multiple integrals. Section 4 dicusses the numerical implementation of the analytical results provided in Section 2 and Section 3.

Valuation formulae for discretely monitored lookback options

The option life is divided into n intervals, n being the number of fixing dates. We start by showing how to value a fixed-strike discrete lookback option with three fixing dates.

The method will then be easily extended to a greater number of fixing dates.

The option life starts at time 0 t and ends at time n t . Thus, our three fixing dates are : Following the risk-neutral valuation approach, at the contract inception 0 t , the undiscounted value of a fixed-strike discrete lookback call option with three fixing dates, is given by :

1 2 3
      1 2 3 0 0 0 1 2 3 , , , Q t t t C S K t t t t E S S S K S              (1)
where Q E is the expectation operator under the equivalent martingale measure and K is the strike price. .

Let

1 , the conditional expectation in (1) can be expanded as :

    2 1 1 1 1 1 1 1 1 2 1 3 /2 0 , , t t t t t t r t W Q S K S S S S E S e d s s               1          2 2 1 1 1 1 2 2 2 2 1 1 2 1 2 1 2 2 1 2 3 /2 /2 0 , , t t t t t t t t r t r t t W W W Q S K S S S S E S e d s d s s s                      1 (2)               2 2 2 1 1 1 1 2 2 2 2 1 3 3 3 3 2 1 2 3 1 2 1 3 2 3 3 1 3 1 /2 / 2 /2 0 , , t t t t t t t t t t r t r t t r t t W W W W W Q S K S S S S S e E d s d s d s s s s                                1       1 1 2 1 3 2 2 1 2 3 3 3 1 3 1 , , ,
                    1 1 1       1 1 1 1 1 2 1 3 0 , r t t t t t t S e Q X k Q X X X X d                1 1 1 2 2 2 1 2 2 1 2 3 0 , r t r t t t t t t t S e Q X k X X Q X X d d          (3)           1 1 1 2 2 2 1 3 3 3 2 3 3 1 3 2 0 , , r t r t t r t t t t t t t S e Q X k X X X X d d d                       1 1 2 1 3 2 2 1 2 3 3 3 1 3 2 , ,
, ,

t t t t t t t t t t t t t t t Q X k Q X X X X Q X k X X Q X X K Q X k X X X X                        
To get to (3), we have used the property of independence of Brownian increments on nonoverlapping time intervals, so that :

      1 1 2 1 3 1 1 2 1 3 , ,
We have also applied Girsanov's theorem by defining a new measure Q  such that :

    1 2 1 1 exp 2 i i t n n i t t i i i i F dQ W W t t dQ s s                      (6)
Let us show how to obtain such probabilities under the Q  measure. A change of drift from :

   2 1 1 / 2 n n i i i i i i r t t m d s        (7) to :    2 1 1 / 2 n n i i i i i i r t t m d s         (8)
provides the corresponding probabilities under the Q   measure.

First,

    1 2 1 3 2 1 3 1 , 0, 0 t t t t t t t t Q X X X X Q X X X X       
is the joint cumulative distribution function of two correlated increments of a geometric Brownian motion with parameters that vary as step functions of time . Both increments are normally distributed :

        2 1 2 2 2 1 2 2 1 t t X X t t t t m s       (9)             3 1 2 2 2 2 1 3 3 2 2 2 1 3 3 2 t t X X t t t t t t t t m m s s           (10)
The pair  

2 1 3 1 , t t t t X X X X   is bivariate normal. The correlation coefficient between   2 1 t t X X  and   3 1 t t X X  is :             2 1 3 1 3 1 2 2 2 1 2 2 2 2 1 3 3 2 2 1 cov , var var t t t t t t X X X X t t t t t t X X X X s r s s           (11)
Thus :

            1 2 1 3 2 2 1 2 2 1 3 3 2 2 2 2 2 2 1 2 2 1 3 3 2 , , t t t t t t t t t t Q X X X X N t t t t t t m m m r s s s                         (12)
where

  2 .,. N r 
refers to the bivariate standard normal cumulative distribution function.

Next, the probability  

2 2 1 , t t t Q X k X X   can be decomposed as :       2 2 1 2 1 1 2 1 , , , t t t t t Q Q t t t X k X k X k X X Q X k X X E E                         1 1       1 2 1 2 1 , 0 t t t t t Q X k X k Q X k Q X X        (13)
The correlation between

1 t X and 2 t X is equal to   2 1 1 2 2 1 1 2 2 1 t t t t s r s s   
, so we have:

  2 2 1 , t t t Q X k X X          2 2 1 1 1 2 2 1 1 1 1 1 2 2 2 1 1 1 1 2 2 1 1 1 2 2 1 , k t t t t t k t k t N N N t t t t t t t m m m m m r s s s s s                                        (14)
Next, we deal with

  3 3 1 3 2 , , t t t t t Q X k X X X X    . The events     3 3 1 , 0 t t t X k X X    and   3 2 0 t t X X 
 are all correlated, so there is no way to get round a triple integral. The correlation between

3 t X and   3 1 t t X X  is :         2 2 3 3 2 2 2 1 1 2 2 2 3 3 2 2 2 1 1 1 t t t t t t t t t s s r s s s         (15) while the correlation between   3 2 t t X X  and   3 1 t t X X  is :       2 3 3 2 2 2 2 3 3 2 2 2 1 t t t t t t s r s s      (16)
Denoting by

  3 1 2 .,.,. , N r r 
the special form of the trivariate normal cumulative distribution function defined in Appendix 1, we obtain :

  3 3 1 3 2 , , t t t t t Q X k X X X X                      3 3 2 2 2 1 1 1 3 3 2 2 2 1 2 2 2 2 2 3 3 2 2 2 1 1 1 3 3 2 2 2 1 3 3 3 2 1 2 3 3 2 , , ,
k
                                      (17)
Retracing our steps, the value of a fixed-strike lookback discrete option with three fixing dates 1 2 , t t and 3 t within the option life  

0 3 t t  is given by : PROPOSITION 1   0 0 1 2 3 , , , V S K t t t t                   3 1
P P e q q              (18)
where :

                  1 2 2 3 2 11 12 2 1 2 2 3 2 3 , k f f f f P N P N f f f f q m q m q m s s s s s                                            1 1 2 1 21 2 1 1 2 1 2 , k f k f f P N f f f q m q m s s s s                               1 2 3 22 1 2 3 k f f f N N P N f f f q m q m q m s s s                                                      1 3 2 3 3 2 3 3 3 3 1 3 2 3 3 1 3 2 3 , , , k f f f f f P N f f f f f q m q m q m s s s s s s s                    1 1 q         if the option is a call if the option is a put    0 1 1 exp n n i i i i i F S r t t d                     1 1 exp n n i i i i r t t e                           2 2 1 1 n m n i i i i i i i i m f t t f t t s s s s                   1 1 n m n i i i i i i i i m f t t f t t m m m m                   1 1 n m n i i i i i i i i m f t t f t t m m m m              
11 21 , P P   and 3 P  are the same as 11 21 , P P and 3 P respectively, except for the drift coefficients, given by : 2 / 2

i i i i r m d s     , instead of : 2 / 2 i i i i r m d s    .
A formula for a fixed-strike discrete lookback call option with four fixing dates can be derived by following the same steps as with three fixing dates. The option life starts at time 0 t and ends at time 4 t , so that our four fixing dates are : 1 2 3 4 , , , t t t t , with :

0 1 2 3 4 t t t t t     .
Performing the necessary calculations, it can be shown that the value of such an option is given by : PROPOSITION 2

  0 0 1 2 3 4 , , , , V S K t t t t t                 4 1
P e q q                 
31 32 4

P P P   (19)
where :

                  1 2 2 3 2 4 11 12 3 1 2 2 3 2 4 , , k f f f f P N P N f f f f q m q m q m q m s s s s                                                  2 2 3 1 1 2 1 21 2 2 3 2 4 1 1 2 1 2 , , f f k f k f f P N f f f f f s s q m q m s s s s s s                                               1 2 3 3 4 3 22 2 1 2 3 3 4 3 4 , k f f f f f N N P N f f f f f q m q m q m q m s s s s s s                                                          1 3 2 3 3 2 3 3 31 3 1 3 2 3 3 1 3 2 3 , , , k f f f f f P N f f f f f q m q m q m s s s s s s s                                       4 1 4 2 4 3 4 32 4 4 4 1 4 2 4 3 4 , , , f k f f f P N P N f f f f q m q m q m q m s s s s                                         4 2 4 3 4 4 4 1 4 2 4 3 4 , , f f f f f f f f q m s s s s s s s           11 21 31
, , P P P    and 4 P  are the same as 11 21 31 , , P P P and 4 P respectively, except for the drift coefficients, given by : 2 / 2

i i i i r m d s     , instead of : 2 / 2 i i i i r m d s   
. The other notations are defined in Proposition 1.

Appendix 2 provides a compact expression for the value of a fixed-strike discrete lookback option with n fixing dates, for any n such that n   and 4 n  .

Knowing a formula for the value of a fixed-strike discrete lookback option, we can easily deduce the value of a floating-strike discrete lookback option. Indeed, the well known parity relation between fixed-strike and floating-strike lookback options in a continuous framework is readily transposed into a discrete framework. If we denote by  

0 0 1 2 , , ,..., n C S K t t t t  
the value of a discrete lookback call and by

  0 0 1 2 , , ,..., n P S K t t t t  
the value of a discrete lookback put, we have the following parity relations :

        0 0 0 0 1 2 1 2 0 1 1 1 1 , , ,..., , , ,..., exp exp Fixed floating n n n n i i i i i i i i C S K t t t t P S K t t t t S t t K r t t d                                          (20)         0 0 0 0 1 2 1 2 0 1 1 1 1
, , ,..., , , ,..., exp exp

Fixed floating n n n n i i i i i i i i P S K t t t t C S K t t t t S t t K r t t d                                          (21)

Analytical valuation of discretely monitored barrier options

The same model assumptions and notations as in section 2 are used here. In particular, in

each interval   1 , i i t t 
, the underlying asset follows a geometric brownian motion with a given set of piecewise constant parameters. We first seek an analytical formula for the value of a European-style up-and-out barrier option, with two different barrier levels 1 H and 2 H discretely monitored at times 1 t and 2 t respectively ; 0 t is the option contract inception, 3 t is the contract expiry and we have : 0 1 2 3 t t t t    . Following the riskneutral valuation approach, the value of such an up-and-out call option is given by :

        1 1 2 2 1 3 3 2 3 1 2 1 2 0 , t t r t r t t r t t Q t S H S H e E S K S                  1       1 1 2 2 1 3 3 2 1 2 3 0 1 2 , , t t t t t t t t S e Q X h X h X k d d d                 1 1 2 2 1 3 3 2 1 2 3 1 2 , , r t r t t r t t t t t e K Q X h X h X k          (22) with   0 ln / i i h H S 
; the other notations are the same as in section 1.

By conditioning and applying the Markov property of Brownian motion, simple calculations yield :

  1 2 3 1 2 , , t t t Q X h X h X k      3 2 1 2 1 2 2 , 1 t t Q Q t t X h X h E E X k X h                    (23)     2 2 2 2 1 1 1 1 2 2 2 1 1 1 3 3 2 3 3 2 1 1 2 2 1 2 2 y x t
                                                    (24)         2 1 1 2 2 1 1 1 2 2 1 3 3 2 1 1 1 1 1 2 2 3 3
, ,

h
m m m m m m s s s                   (25)
where :

           2 2 1 1 2 2 2 2 2 3 3 1 2 2 3 2 2 1 / 2 1 / 3 / 2 1 2 2 3 , , 2
p                                      (26)
The function

  1 2 1 1 2 1 , ,..., , n n t t t t n Q X h X h X h X k         1 2 1 1 2 1 1 2 1 ... , ,..., , ... n n n n n n x x x x dx dx dx dx                  (27)
where :

            2 1 1 2 2 1 1 1 1 1 2 1 1 2 2 1 1 1 2 2 1 1 1 2 1 1 1 1 1 2 2
                                                   2 2 2 1 1 2 2 1 1 2 2 1 1 2 1 1 2 2 1 1 1 2 2 1 1 1 2 1 ... 2 2 1 / 2 1 / 2 1 / / 2 1 2 2 1 1 , ,..., , 2 1 /
p                                                                        (28) 1 1 q         if the option is a call if the option is a put > <         if the option is a call if the option is a put
The value of an up-and-in option is easily deduced by subtracting the value of the corresponding up-and-out option from that of a plain vanilla option. As for the values of down-and-out or double knock-out options, they can be written in a similar manner simply by modifying the integration bounds.

Numerical implementation

The formulae provided in Section 2 and Section 3 raise the question of multidimensional integration. Indeed, as the number of monitoring dates increases, so does the dimension of the integrals involved in the computation of the multivariate distribution functions that appear in these formulae.

In low dimensions, this problem can be analytically tackled. The following results provide two exact integration rules that reduce the trivariate integral to a univariate one and the quadrivariate one to a bivariate one. Let

                             (
                                                    (30)
The numerical integrations implied by ( 29) and ( 30) are very easy to perform. A level of at least 6 10  accuracy, which is more than enough for option pricing, can be achieved with a mere 16-point Gauss-Legendre rule (and a lower bound of -8.5 in the integral), which is extremely efficient.

However, as dimension rises, such analytical simplifications seem more difficult to carry out. To adress the issue of multidimensional integration in high dimension, one can turn to Monte Carlo integration, the convergence rate of which is independent of dimension.

To compute the m  dimensional integral

      2 2 1 1 2 1 1 1 2 1 1 1 , 1 ,..., 1 j j j j j j j m m m m m w y w w y w w y r r r r           (31)
and then test the relevant conditions for each deviate in each sample :

  1 1 2 2 , , ..., , 1,..., j j j m m w x w x w x j n    
By the law of large numbers, repeating this procedure will achieve convergence to the exact value of

  1 2 1 2 1 , ,..., , ,..., m m m N x x x r r r  
as n tends to infinity. The convergence rate is in the order of 1/ n . Besides, this method is easily combined with classical variance reduction techniques such as antithetic variates, or stratified and Latin hypercube sampling.

Let us now implement the formula for a discretely monitored fixed strike lookback call option. For this first example, we use the same volatility, riskless rate and dividend rate parameters during the entire option life. This will enable us to draw a comparison with the exact value of a fixed strike lookback call with continuous monitoring, which admits a well-known explicit formula. We assume volatility is equal to 32%, while the riskless rate is 5%, the dividend rate is 1.5% and the option expiry is one year. The spot value of the underlying asset at the contract's inception is 100 and the option is at-the-money. Under those assumptions, the exact value of a plain vanilla call is 14.074 (in a Black-Scholes framework, and to four decimal places), while the exact value of a fixed strike lookback call with continuous monitoring is equal to 28.634. Let us now examine how a finite number of fixing dates affects the option value compared with the case in which there is no fixing date at all before expiry (i.e., the case of a plain vanilla call), which is a lower bound, and the case in which there is an infinite number of fixing dates (i.e., the case of a fixed strike lookback call under the assumption of continuous monitoring), which is an upper bound. If we set four evenly spaced fixing dates during the option life (that is, the first fixing occurs at t = 0.25), the option value increases from to 14.074 to 19.732, applying the analytical formula given by Proposition 2, and following the implementation rule given by eq. ( 30). This value is obtained in less than one second 1 , which is much faster than the 26 seconds required to obtain an approximate value of 19.761 after performing 2,000,000 Monte Carlo simulations using antithetic variates and stratified sampling. This is a large increase in value (roughly 40%), which suggests that the bulk of the price transformation is concentrated in the first few fixing dates. This observation is confirmed by computing the option value with eight evenly spaced fixing dates (first fixing at t  0.125), which yields 22.045. Thus, setting only eight evenly spaced fixing dates accounts for 54.5% of the total increase in value that results from the transformation of a plain vanilla contract into a contract with an infinite number of fixing dates. This is obviously linked to the fact that our fixing dates uniformly fill the option life. If we set all eight fixing dates during, say, the sixth month of the option life, the option value diminishes to 18.7 (and to 17.2 if we set all eight fixing dates during the first month of the option life).

Let us now turn to higher monitoring frequencies of the underlying asset. Monte Carlo integration will tend slowly to exact option values when the number of fixing dates becomes large. This efficiency issue is more acute for discretely monitored lookback options than for discretely monitored barrier options. Indeed, the formulae for discretely monitored lookback options are more complicated and they involve the computation of many more high dimensional integrals. Implementing the formula for a one-year fixed strike lookback call with daily monitoring or simply weekly monitoring would be slow and quite cumbersome to code. Fortunately, all that is needed to obtain a fast approximation for the value for our fixed strike lookback call with an arbitrarily large number of fixing dates is to compute its value with four and eight fixing dates, which is easy and accurate using the analytical formulae provided in this paper, and to know the value of the two limiting cases (plain vanilla call and fixed strike lookback call with continuous monitoring), for which there are explicit formulae available. Then, all fixed strike lookback call values with a number of fixing dates greater than eight can be interpolated by means of quintic splines.

Let us denote by   s x the logarithm of the value of a fixed strike lookback call with x e fixing dates (with a given set of spot, strike, volatility, riskless rate, dividend rate and expiry parameters). Let us approximate continuous monitoring by 12500 fixing dates in a year (which corresponds to 50 monitoring times per business day 2 ). Resuming the inputs of our first example above, we have :

  0 2.644351543 s  ,   1.386294361 2.982647 s  ,   2.079441542 3.09285898 s  ,   9.433483923 3.35461126 s 
If we construct a quintic spline out of these data, we obtain 3 , for 2.079441542 x  (that is, for more than eight fixing dates), the following polynomial of degree 5 :

2 Adding more fixing dates to the approximation has insignificant effet on the spline value 3 All the spline polynomials reported in this section were obtained using the spline function in the Maple software .

  2 3 4 5 s x =2.631426028+0.3324445005x-0.659113161x +0.006854474588x -0.0003633055742x +0.000007702468719x (32)
To obtain the value of a fixed strike lookback call with daily monitoring, for example, it suffices to take the exponential of The average discrepancy between interpolated values and Monte Carlo simulation values is equal to 0.36%. Note that it is extremely time-consuming to obtain a reliable Monte Carlo simulation estimator of the fixed strike lookback call with 250 fixing dates, not only because of the number of pseudo random numbers to be drawn, but also because of the very large number of conditions to be tested.

Our second example involving discretely monitored lookback options provides a basic illustration of the advantage of using a model that allows for different parameters over various sub-intervals of the option life. We have observed that setting all fixing dates in the sixth month of the option life substantially decreases the option value in a framework in which volatility remains constant during the entire option life. Now suppose that volatility jumps to 0.75% in the sixth month. Applying the formula provided in section 2, our option value rises from 18.7 (case where volatility is constantly 32%) to 24.1, while the option value when all fixing dates are uniformly spread throughout the option life is 26.5 with this new term structure of volatility. This stylised example hints at the mispricings we are exposed to if our model cannot allow for volatility, interest rate and dividend rate parameters to vary.

Next, let us move on to discretely monitored barrier options. We focus on an up-and-out put option written on an underlying with spot value 100 S 

, and four different barriers t t  . Those contract specifications are designed so that the distance between the upper barrier and the lower barrier widens as volatility increases. This option better the value of a plain vanilla call minus the value of a double knock-out call with the same specifications. Implementing our semi-analytical approach, we find a value of 1.60 with eight monitoring dates, 3.13 with sixteen monitoring dates, and 3.77 with twenty-four monitoring dates.

Conclusion

In view of the variety of valuation methods dealing with discretely monitored lookback and barrier options, it is not always an easy task to choose the most appropriate one for a given problem. We should like to close this paper with the following recommendations :

If you need to price a contract with complex payoff specifications (such as a many different single or double barriers) or non-standard monitoring specifications (such as unequally spaced monitoring dates), use : If you want to price a contract with simple payoff and monitoring specifications :

(i) Continuity adjustment is the fastest and most straightforward method available, but its accuracy depends on the type of option being valued as well as on the number of monitoring dates (ii) Sequential quadrature is extremely efficient too, but its implementation is lengthy and rather complicated (iii) Robust lattice methods and explicit formulae combined with spline interpolation appear to us to be a good trade-off between accuracy and efficiency If you need to price an American-style contract with early exercise features, use robust lattice methods.

Appendix 1

In general, if

  1 ,..., n X X X 
is a vector of n joint standardized normal random variables with a symmetric, positive-definite   n n  matrix of variances and covariances  , then the density of X is given by :

      1 2 1 /2 ,..., 2 T X X n n e f x x Det p      (36)
where T X is the transpose of X ,

 

Det  is the determinant of  and 1   is the inverse of  (for a proof, see e.g., Tong, 1990).

However, using these general multinormal expressions becomes analytically cumbersome and computationally inefficient as the dimension of the integral rises. Actually, simplified expressions can be used when dealing with the finite-dimensional distributions of geometric Brownian motion GBM. Let   1 2 , ,..., n t t t be n successive points in a time interval. Then, it is easily shown [START_REF] Guillaume | Window Double Barrier Options[END_REF] that :

  , 1 2 1 2 ,..., n n t t t P X x X x X x               2 1 2 1 2 1 1 1 2 1 ... ... 1 2 , 1 2 1 , 2 2
           (38)
where :

              1 2 2 3 11 12 1 1 2 2 3 , ,..., n k f f f P N P N f f f q m q m q m s s s                                       2 2 2 3 2 1 2 2 3 2 4 2 , ,..., n n n n f f f f f f f f q m s s s s s s s                               1 1 2 1 21 2 1 1 2 1 2 , k f k f f P N f f f q m q m s s s s                           1 2 1 2 k f f N N f f q m q m s s                                           
f f f f f f f f f f P N f f q m q m q m s s s s s s s s s                                                         1 3 2 3 3 2 3 3 31 3 1 3 2 3 3 1 3 2 3 , , , k f f f f f P N f f f f f q m q m q m s s s s s s s                                          
f f f f f f f f f f P N f f q m q m q m s s s s s s s s s                                                                 1
k f f f f f f f f f f P N f f f f q m q m q m q m s s s s s s s s s s                                                         
f f f f f f f f f f P N f f q m s q m q m s s s s s s s s                                   …                                 1 1 2 1 2 1 1 1 2 1 2 1 1 1 1 1 2 1 3 1 1 1 1 1 2 1 2 1 , , ..., , , ,...
k f f f f f f P N f f f f f f f f q m q m q m s s s q m s s s s s s s                                                           1 2 n n n f P N f q m s                                          1 2
k f f f f f f f f P N f f f f f f q m q m q m q m s s s s s s s s s s                                    
The mn P  terms are the same as the mn P terms, except for the drift coefficients, which are :

2 / 2 i i i i r m d s     , instead of : 2 / 2 i i i i r m d s   
. The other notations are defined in Proposition 1.

  (i) Monte Carlo simulation if you need to incorporate several stochastic factors into the equation modelling the motion of the underlying asset (ii) Analytical formulae if you can model the motion of the underlying asset as a continuous process with piecewise constant parameters (in particular, if you want to use the market information contained in volatility surfaces and in the term structure of interest rates)

  is the value of the underlying asset at time t . In each interval  

	, , t t t with : 0 t = contract inception 1 t    = expiry. Any spacing can be chosen 2 3 t t
	between the fixing dates.
																	i t	1 , i t
	dS	t		i r		i d	t S dt		i s	t S dW t	,	t		i t	1 , i t	
																	i t	1 , i t
														i t	1 , i t
																	, , F P 		, with
	t F			s	 W s t , s 		being the natural filtration of t W .

t S  , it is assumed that the underlying asset follows a geometric brownian motion under the equivalent martingale measure denoted by Q :     where : i r is the risk-free interest rate in   1 , i i t t  i d is the dividend rate paid out by the underlying asset in    i s is the volatility parameter in    t W is Brownian motion defined on a probability space 

  ,

	1	t S		t K S S S S t t t  				1	t S		K	1	t S	t S S S t t  	(4)
	and :
	 1	2 t S		, K S	2 t		1 t S S , t 2		t S	3			 1	2 t S		2 t K S ,	1 t S 	  1	2 t S	t S 	3	

  

				4	,	4 5	,...,	4	n	4	,	4 5	,...,
				4		4 5		4	n	4 5	4 6
	32	n	3
			4	n	1
			4	n

  ,

		n		n	n	n
		n		n	n	n
	n	n
		n	n	n	n
		n	n	n	n	n

All reported values were obtained on a

2.4 Ghz -clock single computer

Those splines are of degree 3, instead of 5 in the lookback option case first examined. This is because the gradual loss in the knock-out option value as the number of monitoring dates increases is smoother than the increase in the lookback option value as the number of fixing dates increases. The following tables, where all monitoring dates are uniformly spread throughout the option life, report a few numerical results : protects investors when volatility is high, by reducing the likelihood of knocking-out, while it diminishes protection when volatility is low, thus reducing the overall cost of the option. This ability to adapt to the fluctuations of volatility is one of the main attractions of « corridor » options, so that it would be particularly inappropriate to value them with a model that would set all parameters constant during the entire option life. Table 5 presents a few numerical results, obtained by implementing the analytical representation provided in Section 3 along with the Monte Carlo integration scheme explained at the beginning of this Section. We could have constructed a polynomial interpolation again. But, as this option admits no known explicit solution under the assumption of continuous monitoring, we would have had to begin by numerically integrating the formula with discrete monitoring for a very large number of monitoring dates, making it an approximation to the continuous monitoring limit, which would have been very slow. In practice, this will be a worthy time investment only if a lot of option values with discrete monitoring are subsequently required.

Eventually, we illustrate the case of a double knock-in call option, with an up-barrier set at 150 and a down-barrier set at 60, and all monitoring dates set within the last quarter of the option life, during which volatility jumps to 45% after being equal to 15% during the first three quarters (riskless rate = 5%, dividend rate = 1.5%, spot value = 100, strike price = 100, expiry = one year). As mentioned in section 3, this option value must be equal to

Appendix 2

The no-arbitrage value of a fixed-strike discrete lookback option with n fixing dates, for any n such that n   and 4 n  , is given by : Sullivan, M.A. (2000), Pricing discretely monitored barrier options, Journal of computational finance, 3 (4), [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52]


Frédéric Gueydan

Jacques Précigout

Modes of continental rifting as a function of ductile strain localization in the lithospheric mantle

Analogue and numerical models have shown that the strength of the lithospheric mantle controls the mode of lithosphere deformation. In extension, the presence or absence of a high strength brittle mantle respectively leads to localized or distributed rifting. However, first order geophysical data question the existence of such a brittle mantle. Here we use 2-D finite-element large strain modelling to quantify the impact of a ductile localizing mantleinstead of brittlein triggering continental rifting. As a novelty, the mantle rheology considers the effect of grain boundary sliding during strain-induced grain size reduction, which may promote a significant strength drop and subsequent strain localization at low mantle temperature (< 700-800°C). Our results reveal that such ductile localizing mantle implies varying modes of continental rifting that mainly depend on both the amount of weakening in the ductile mantle and the strength of the lower ductile crust. A medium to strong lower crust implies coupling between the upper crust and ductile localizing mantle, yielding to narrow continental rifting. In contrast, a weak lower crust implies decoupling between the upper crust and ductile localizing mantle, giving rise to a switch from distributed faulting at incipient strain to localized faulting at large strain. Ductile strain localization in the lithospheric mantle is therefore sufficient to trigger continental rifting, although a critical amount of weakening is required. Such ductile localizing mantle provides a relevant geological and mechanical alternative to the brittle mantle. It moreover provides a wider variety of modes of upper crustal faulting that are commonly observed in nature.

Introduction

Thermal and rheological layering of the continental lithosphere exerts a direct control on the mode of lithosphere deformation. While the crustal strength dominates the lithosphere strength for a Moho temperature (T M ) higher than 700°C, the high-strength sub-Moho mantle controls the lithosphere strength for T M lower than 700°C [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. The presence of a high strength uppermost mantle together with a weak deep crust that permits horizontal ductile shearing along the Moho is essential to promote strain localization at lithosphere scale [START_REF] Allemand | Width of continental rifts and rheological layering of the lithosphere[END_REF][START_REF] Brun | Mantle exhumation at passive margins[END_REF][START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF].

In analogue and numerical experiments, the high strength uppermost mantle is modelled as a brittle mantle, such as predicted by the classical strength profile of the lithosphere [START_REF] Brace | Limits on lithospheric stress imposed by laboratory experiments[END_REF][START_REF] Sawyer | Brittle failure in the Upper Mantle during extension of continental lithosphere[END_REF][START_REF] Carter | Flow properties of continental lithosphere[END_REF]). Some of these numerical models promote the development of mantle fault zone through an ad-hoc drop of friction coefficient with increasing strain [START_REF] Huismans | Symmetric and asymmetric lithospheric extension: Relative effects of frictional-plastic and viscous strain softening[END_REF][START_REF] Huismans | Effect of plastic-viscous layering and strain softening on mode selection during lithospheric extension[END_REF][START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. [START_REF] Bos | Frictional-viscous flow of phyllosilicate-bearing fault rock: Microphysical model and implications for crustal strength profiles[END_REF] predict indeed a decrease of analogue fault strength with increasing strain. This weakening process can be applied to upper crustal faulting and could well explain the weakening of major lithospheric faults [START_REF] Holdsworth | Weak Faults-Rotten Cores[END_REF]. However, the applicability of a brittle weakening process to the uppermost mantle remains to be demonstrated. Furthermore, the existence of a brittle mantle is now questionned. The rare occurrence of earthquake within the continental lithospheric mantle [START_REF] Maggi | A reassessment of focal depth distributions in southern Iran, the Tien Shan and northern India: Do earthquakes really occur in the continental mantle?[END_REF], together with the low elastic thickness inferred in deformed region [START_REF] Jackson | Strength of the continental lithosphere: time to abandon the jelly sandwich?[END_REF][START_REF] Audet | Dominant role of tectonic inheritance in supercontinent cycles[END_REF], suggest in opposite that the crustal strength controls the lithosphere strength. These features contradict therefore the existence of a high strength brittle uppermost mantle, but they also conflict with pre-requisites of the lithosphere dynamics that involve a high strength and localizing mantle to explain both lithosphere-scale strain localization and the large values of Elastic thickness beneath cratons [START_REF] Burov | The long-term strength of continental lithosphere: "jelly sandwich" or "crème brûlée[END_REF].

As an alternative, an entirely ductile mantle with an initial high strength that decreases during deformation can reconcile these opposite views of lithosphere strength. Indeed, [START_REF] Frederiksen | Numerical modelling of strain localisation during extension of the continental lithosphere[END_REF] have shown that mantle shear zones can be generated in response to strain-induced weakening of the uppermost lithospheric mantle, providing a simple explanation for lithosphere-scale strain localization. Based on structural and microstructural observations in the Ronda peridotite and using available olivine flow laws [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF], we have shown that viscous deformation dominated by dislocation-accommodated grain boundary sliding (disGBS) can induce a significant weakening during dynamic recrystallization at low mantle temperature [START_REF] Précigout | Strain localisation in the subcontinental mantle --a ductile alternative to the brittle mantle[END_REF]. 1-D aggregate-scale numerical modelling of the disGBS-induced weakening moreover shows that this process is efficient to strongly reduce the mantle strength for low mantle temperature [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. These results exemplify the role of ductile strain localization within the ductile mantle and allow us to define a ductile localizing mantle.

Nonetheless, the impact of such a ductile localizing mantle on lithosphere-scale deformation remains to be constrained. This is the main purpose of the present study through 2-D numerical analyses focusing on 1/ the role of the amount of weakening in the ductile mantle and 2/ the impact of ductile crust rheology for the mode of continental rifting.

Continental lithosphere rheology

Mantle rheology

We account here for a composite olivine rheology that combines four deformation mechanisms with their respective flow laws: dislocation creep (r), diffusion creep (d), disGBS (g) and low-temperature plasticity (e) [START_REF] Goetze | The mechanisms of creep in olivine[END_REF][START_REF] Drury | Dynamic recrystallization and Strain Softening of Olivine Aggregates in the Laboratory and the Lithosphere[END_REF][START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF][START_REF] Précigout | Strain localisation in the subcontinental mantle --a ductile alternative to the brittle mantle[END_REF][START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. Each mechanism contributes to the overall strain rate ( e ) of an olivine aggregate as follows:

e = e r +e d +e g +e e (1) With Dislocation creep e r = A r .exp(

-Q r RT ).t n r (2)
Diffusion creep

e d = A d .exp( -Q d RT ).t n d .d -m d (3) disGBS e g = A g .exp( -Q g RT ).t n g .d -m g (4)
and Low-temperature plasticity (or exponential creep)

e e = A e .exp[-(Q e / RT).(1-t / t p ) n e ] (5)
where, e , , d, m, T, R, A, Q, and n are, respectively, the strain rate (in s -1 ), the equivalent shear stress (in MPa), the grain size (in µm), the grain size exponent, the temperature (in K), the gas constant (J/°), the pre-exponential constant (MPa n /s), the activation energy/enthalpy (in J, the effect of pressure is disregarded here) and the stress exponent. The subscripts r, d, g and e stand for respectively dislocation creep, diffusion creep, disGBS and exponential creep.

p  is the Peierls stress defined for low-temperature plasticity (Goetze, 1978).

The slopes for the iso-temperature curves (Fig. 1). This degree directly depends on the grain size exponent (m) of the olivine flow law (eqs. 3 and 4).

During strain-induced grain size reduction of olivine, the contribution of GSS mechanisms to the deformation of mantle rocks is enhanced at the expense of GSI mechanisms. As a consequence, olivine aggregate suffers a stress drop/weakening that may trigger strain localization depending on both the weakening rate and temperature [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. If only grain size reduction by dynamic recrystallization is involved, the weakening induced by GSS processes occurs until equilibrium of grain size. The nature of this equilibrium is still subjected to some hypotheses, which imply either a balance at the boundary between the stress-grain size fields of GSS and GSI mechanisms (dashed thick line in Fig. 1; balance hypothesis; [START_REF] De Bresser | On dynamic recrystallization during solid state flow: Effects of stress and temperature[END_REF]2001) or an experimentally constrained paleo-piezometer (dashed thin line in Fig. 1; Piezometer hypothesis; Twiss, 1997;[START_REF] Karato | Dynamic recrystalization of olivine single crystals during hightemperature creep[END_REF][START_REF] Van Der Wal | Relationships between dynamically recrystallized grain size and deformation conditions in experimentally deformed olivine rocks[END_REF]. We will consider these two hypotheses in this study (Fig. 1).

The rate of dynamic grain size reduction is intimately related to the diffusioncontrolled migration processes of intra-crystalline dislocations and is therefore quite complex to quantify. Several theoretical and lab-based grain size evolution laws have been proposed

and discussed so far (e.g. [START_REF] Montesi | Grain size evolution and the rheology of ductile shear zones: from laboratory experiments to postseismic creep[END_REF], but this level of discussion is beyond the scope of this paper. We thus chose to follow [START_REF] Braun | A simple parametrization of strain localization in the ductile regime due to grain size reduction: A case study for olivine[END_REF] and [START_REF] Kameyama | The interaction of viscous heating with grain-size dependent rheology in the formation of localized slip zones[END_REF] using the following phenomenological grain size evolution law to simulate dynamic grain size reduction:

d = -e .(d -d ¥ ) / e T , ( 6
)
where d is the changing rate of grain size, e the overall strain rate, d the mean grain size of an olivine aggregate, d  the recrystallized grain size and T  the characteristic strain required to approach the recrystallized grain size [START_REF] Braun | A simple parametrization of strain localization in the ductile regime due to grain size reduction: A case study for olivine[END_REF]Montési and Hirth, 2003). The value of T  has been set to 40% in order to be consistent with our previous modelling [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. The value of T  defines the time/strain necessary to reach strain localization, but does not change the overall pattern of deformation. However, future studies that compare numerical modelling of shear zone formation and large strain experiments should be done in order to better constrain this parameter. As an example, Hansen et al. (2012) experimentally obtained T  of 100%. Finally, the recrystallized grain size will be defined by either the recrystallized piezometer (Piezometer hypothesis) or the boundary between diffusion creep and disGBS/dislocation creep (Balance hypothesis, Fig. 1).

Studied rheological parameters.

For a given continental geotherm, we wish studying the role of 1/ the ductile localizing mantle, and 2/ the lower crustal rheology in defining the modes of continental rifting. For that purpose, we will document in our 2-D numerical models the role of the selected disGBS flow law, the recrystallized grain size hypothesis and the ductile crust rheology.

The drop of stress/strength (e.g. weakening) during grain size reduction will depend on the selected disGBS flow law, which mostly differs in terms of grain size exponent (m, eq. 3).

So far, two different rheologies have been proposed to account for the dislocationaccommodated Grain Boundary Sliding [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF][START_REF] Hansen | Grain boundary sliding in San Carlos olivine: Flow law parameters and crystallographic-preferred orientation[END_REF]. In ). For H&K rheology, the amount of weakening increases for decreasing temperature, promoting an ideal potential for strain localization at low mantle temperature [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. In contrast, a weakening is still predicted for the H rheology, but the slop of iso-temperature curves is significantly lower, potentially inhibiting the potential of strain localization. So far, only the rheology of H&K has been tested in terms of potential of strain localization, and never at lithosphere scale [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. We thus consider here the effect of these different rheologies on the mode of lithosphere extension.

The selected hypothesis for the final grain size (Balance vs. Piezometer) leads much smaller differences in amount of weakening (Fig. 1). The balance hypothesis implies more weakening than the piezometer hypothesis. The 2-D numerical models presented below aim at quantifying the impact of such variations in the amount of weakening at lithosphere scale.

The ductile crust strength controls the amount of coupling/decoupling between the crust and the mantle [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. In order to study this impact in the present study, we have selected two end-members: a medium/strong ductile crust [START_REF] Hirth | An evaluation of quartzite flow laws based on comparisons between experimentally and naturally deformed rocks[END_REF]) and a weak ductile crust [START_REF] Kirby | Rheology of the lithosphere: selected topics[END_REF]. In table 2, we provide the varying parameters used for different set of 2-D numerical models.

Numerical set-up

Selected continental geotherm

To address the effect of the crust/mantle rheologies on lithospheric deformation, we have selected here one continental geotherm typical of the thermal conditions that can lead to lithospheric necking. In previous studies [START_REF] Buck | Modes of Continental Lithospheric Extension[END_REF][START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF], it has been shown that a cold and stiff lithosphere (Moho temperature T M lower than 400°C) induces a coupling between crust and mantle deformation, and hence no major necking occurs. A hot and weak lithosphere promotes wide rifting (T M > 700°C) without major necking either. In contrast, for T M values between 400 and 700°C, the high strength mantle, which is mostly brittle in previous studies, controls the lithosphere deformation; the deep crust is also sufficiently weak to produce deep crustal décollement. In these conditions, mantle shear zones develop and trigger lithosphere necking. However, several geophysical features highly question the presence of such a brittle mantle (see introduction). We thus consider here the impact of a ductile mantle on the mode of continental rifting assuming a Moho temperature of T M = 500°C and considering the occurrence of GBS-related rheology in the uppermost mantle. The strength of the lithospheric mantle will be therefore changing with time, depending on changes in grain size.

Model set-up and boundary conditions

The model set-up is close to that used in previous modelling. Numerical results are presented for a Moho temperature of 500°C. The initial grain size is 3000 µm in the overall lithospheric mantle, which implies the dominance of dislocation creep for temperatures higher than ~650°C and the dominance of exponential creep for lower temperatures at the onset of deformation (Fig. 1A). Grain size reduction during lithosphere deformation will change the dominant deformation mechanism in the lithospheric mantle. Such a change may imply strain localization at lithosphere scale, as it will be now discussed in 2-D.

Results

The reference model (model number 1, Table 2) is first presented in figure 3. It considers a strong ductile crust [START_REF] Hirth | An evaluation of quartzite flow laws based on comparisons between experimentally and naturally deformed rocks[END_REF], the disGBS mantle flow law from [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF] and the balance hypothesis for the recrystallized grain size (Table 2).

Models 2 to 4 will be next presented in order to document the role of the amount of weakening in the mantle (Models 2 and 3) and the role of the ductile crust rheology (model 4).

Strain, grain size and deformation mechanisms during rifting (reference model)

In figure 3, we show the finite strain through the entire lithosphere, and the grain size with the percentage of disGBS within the lithospheric mantle (Table 2). Three steps of extension are shown: after 60, 100 and 140 km of horizontal displacement, e.g. after 6, 10 and 14 Ma of extension. After 60 km of extensional displacement, strain becomes heterogeneous in the lithospheric mantle, with a region of higher strain in the model centre. DisGBS becomes the dominant deformation mechanism to accommodate grain size reduction in the model centre at depths between 30 and 40 km (between 50 and 75% of disGBS, Fig. 3). This implies a more efficient grain size reduction, associated with significant weakening (Fig. 1A), which results from the higher strain in the model centre. The onset of upper crustal faulting occurs in a narrow region above the high strain zone in the mantle. This suggests a strong mechanical coupling between the upper crust and mantle, likely permitted by high stresses in the ductile crust [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF].

These features are more pronounced with an increase in the amount of horizontal extension. Both the coupling between crust and mantle, and the mantle strain localization promote necking of the upper part of the lithosphere, including the crust and uppermost mantle. A set of conjugated mantle shear zones develops after 140 km of extension with very large strain and very small grain sizes. In shear zone centre, disGBS and diffusion creep equally contribute to the deformation because of fully recrystallized rock (50% of disGBS, Fig. 3). These results exemplify the role of ductile strain localization in triggering continental rifting. They show also the role of a medium/strong ductile crust in coupling upper crustal deformation and mantle deformation, as required to trigger narrow rifting [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. Our findings therefore highlight this ductile localizing mantle as a relevant mechanical alternative to the brittle mantle to account for strain localization in the lithospheric mantle, and hence, it should be now considered for future modelling of the continental lithosphere deformation.

Role of the amount of mantle weakening

In Figure 4, we show the finite strain of the whole lithosphere after 140 km of horizontal extension for 1/ the reference model (Fig. 4A The model 2 shows a very limited amount of heterogeneous deformation in the lithospheric mantle (Fig. 4B), which can be explained as follows. The dataset of Hansen et al.

(2011) implies a grain size exponent of 0.7 instead of 2 for the reference model [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF]. As a consequence, grain size reduction will lead to differences in the amount of stress drop/weakening. As an example, a division of the grain size by a factor of 2 leads to a division of stress by 4 for the reference model out of 1.62 for the Model 2. These strong differences in the amount of weakening during mantle grain size reduction (also shown in Fig. 1) thus imply strong differences for strain localization in the ductile mantle. The very low value of weakening for Model 2 explains the absence of strain localization at lithosphere scale. We show therefore that a minimum amount of weakening in the ductile mantle is required to trigger mantle strain localization and subsequent continental rifting. Hence, our results indicate that defining precisely the disGBS flow law is crucial to account for the mode of lithosphere extension. In contrast, the selected hypothesis for the equilibrium of the recrystallized grain size (Fig. 4C) has a very limited impact on the pattern of lithosphere deformation, mostly because it implies limited differences in the amount of mantle weakening (Fig. 1, and discussion in the Rheology section).

Role of the ductile crust rheology

In figure 5, we show the mode of lithosphere extension considering a weaker lower crust than in the reference model [START_REF] Kirby | Rheology of the lithosphere: selected topics[END_REF]. We show here three steps of extension because the pattern of lithospheric deformation drastically changes.

After 60 km of extensional displacement, distributed rifting occurs in the upper crust.

While the rift spacing is controlled by the crustal thickness, the rift width is controlled by the thickness of the upper brittle crust, as previously showed in [START_REF] Allemand | Width of continental rifts and rheological layering of the lithosphere[END_REF] and [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. The weak deep crust sustains shearing parallel to the Moho that remains flat, and the lithospheric mantle is weakly deformed. These features characterize a strong decoupling between the crust and mantle (Fig. 5). At this level of extension, the olivine grain size is of the order of 1500 µm with lower values in the middle of the mantle because of little increase in mantle deformation in this region. Such an olivine grain size implies that GSI mechanisms mostly accommodate mantle deformation: exponential creep below 650°C and dislocation creep above 650°C. As a consequence, the contribution of disGBS (less than 30 %, Fig. 5) is not sufficient to promote weakening and strain distributes over the whole lithosphere.

After 100 km of extensional displacement, the disGBS becomes the dominant mechanism to accommodate mantle deformation in the model centre at depths between 30 and 40 km (between 50 and 75% of disGBS, Fig. 5). This implies a more efficient grain size reduction, associated with significant weakening (Fig. 1). As a consequence, strain becomes heterogeneous in the mantle with a region of higher strain in the model centre, where the disGBS is more efficient. The increase in mantle deformation in the model centre is also accommodated in the deep crust by horizontal shearing. The onset of heterogeneous mantle deformation and deep crustal shearing both slightly modify the faulting pattern in the upper crust, where two main rifts accumulate strain in the model centre.

After 140 km of extensional displacement, intense strain localization occurred in the model centre and formed two conjugated lithosphere-scale mantle shear zones, where the strain has reached values as large as 10. The deep crustal décollements are also more efficient with sense of shearing towards the mantle shear zones. In the core of these mantle shear zone just below the Moho, dynamic recrystallization is complete. DisGBS and diffusion creep again equally contribute to the deformation of this fully recrystallized zone (50% of disGBS, Fig. 5). The intense strain in the mantle together with the deep crustal décollement drastically changes the pattern of upper crustal faulting. Only two upper crustal faults now accommodate strain, defining a narrow rift that leads to lithosphere necking.

These results exemplify the role of the lower crust in defining the level of coupling between the crust and mantle, and hence the modes of upper crust faulting [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. A strong deep crust (Model 1) implies mechanical coupling between the crust and mantle, which leads to narrow rifting. In contrast, a weak deep crust (Model 4) leads to a switch from distributed to localized faulting during progressive strain localization in the lithospheric mantle.

Implications for the mode of continental rifting

The above results show that the ductile localizing mantle allows a wide variety of modes of continental rifting. Such a difference of mode was not possible to predict with a single geotherm considering a brittle mantle only, which implies localized fracturing at early stage of rifting [START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. [START_REF] Schueller | Brittle-ductile coupling: Role of ductile viscosity on brittle fracturing[END_REF][START_REF] Schueller | Mechanics of the transition from localized to distributed fracturing in layered brittle-ductile systems[END_REF]

Conclusion

In this paper, through 2-D numerical modelling of continental extension, we have investigated the role of 1/ a grain-size-sensitive ductile rheology (disGBS) as an alternative to the brittle rheology for the uppermost mantle, and the impact of 2/ strong vs. weak rheologies of the lower crust, both on the modes of continental rifting. We summarize here our main conclusions:

1/ A ductile localizing mantle appears to be relevant to model continental lithosphere rifting.

It moreover allows predicting a wide variety of modes of continental rifting, depending on the mantle rheology involved.

2/ A weakening induced by disGBS does not necessarily promote strain localization. A critical amount of weakening during grain size reduction is thus required to trigger lithosphere-scale strain localization. Future large strain deformation of olivine aggregates should be performed in order to better identify disGBS flow laws and the related weakening.

3/ The strength of the ductile crust defines the degree of mechanical coupling between the upper crustal faulting and lithospheric mantle deformation. A weak ductile crust triggers distributed faulting at early stages of rifting, and then permits a progressive evolution from distributed to localized faulting. In contrast, a medium/strong crust promotes a mechanical coupling between crust and mantle that leads to lithosphere-scale necking.

Accounting for such a ductile localizing mantle in future modelling is therefore a challenging issue in order to better constrain the large diversity of types of continental rifting and/or non-volcanic passive margins. Indeed, several geological studies show that faulting during continental break-up in the absence of large magma generation is marked by three successive steps: 1/ distributed faulting that forms a large number of small basins; 2/ progressive localization of deformation within a narrow region that accumulates more strain;

and then 3/ the formation of a single narrow rift (see for example [START_REF] Cowie | Spatio-temporal evolution of strain accumulation derived from multi-scale observations of Late Jurassic rifting in the northern North Sea: A critical test of models for lithospheric extension[END_REF][START_REF] Whitmarsh | Evolution of magma-poor continental margins from rifting to seafloor spreading[END_REF]. However, the precise comparison between a natural case study and our models are beyond the scope of this paper, which remains a theoretical study focusing on the impact of ductile lithosphere-scale strain localization for the mechanics of continental rifting. , 1998) or experimentally constrained paleo-piezometer (thin dashed line; [START_REF] Van Der Wal | Relationships between dynamically recrystallized grain size and deformation conditions in experimentally deformed olivine rocks[END_REF]. Values of the ductile flow laws are given in Table 1. Table 1: Rheological parameters used for the different sets of numerical models. Rheological parameters of [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF] have been corrected following discussion in [START_REF] Hansen | Grain boundary sliding in San Carlos olivine: Flow law parameters and crystallographic-preferred orientation[END_REF].

Figure Captions

  figure 1A and B, we show the deformation maps respectively considering Hirth and Kohlstedt

  Readers may thus referto[START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF] for a detailed description of the brittle rheology (Von Mises approximation). The initial model width is 300 km. A boundary extensional velocity of V=1cm/a is imposed at the right vertical boundary of the model, while the left vertical boundary is only allowed to move vertically. The isostatic equilibrium at the base of the lithosphere is achieved by the presence of an underlying asthenospheric layer whose viscosity is set to 10 21 Pa.s. A small geometrical perturbation for both the crustal brittleductile transition and the Mohoa deflection with amplitude of 100 mis introduced at the model centre to initiate heterogeneous deformation.

  same as Fig.3), 2/ a model that considers the rheology of[START_REF] Hansen | Grain boundary sliding in San Carlos olivine: Flow law parameters and crystallographic-preferred orientation[END_REF] (Model 2; Fig.4B), and 3/ a model that consider the piezometer hypothesis[START_REF] Van Der Wal | Relationships between dynamically recrystallized grain size and deformation conditions in experimentally deformed olivine rocks[END_REF] rather than the balance hypothesis[START_REF] De Bresser | Grain size reduction by dynamic recrystallization: can it result in major rheological weakening?[END_REF] (Model 3; Fig.4C).

  have indeed shown through simple 2-D numerical experiments that the viscosity of ductile layers in brittle-ductile systems controls the numbers of faults in brittle layers, defining the brittle-ductile coupling.The observed variety of upper crustal fracturing reflects the role of brittle-ductile coupling at the scale of the lithosphere. The increase in ductile viscosity decreases the fault velocity (e.g. viscous friction) and implies the nucleation of new faults to accommodate the displacement rate. As a consequence, a high viscosity promotes distributed faulting while a low viscosity implies localized fracturing. In the present 2-D models, the two end-members configurations occur. 1/ When crust and mantle are decoupled (with a low deep crust strength), the homogeneous deformation of the high-strength ductile mantle controls the deformation of the lithosphere at the early stage of extension, and hence, distributed faulting affects the upper crust (high brittle/ductile coupling). During progressive strain localization in the upper lithospheric mantle, the strength of the mantle decreases and implies a decrease of the degree of brittle/ductile coupling, yielding localized faulting. 2/ When crust and mantle are mechanically coupled (medium/strong crust), necking of the upper part of the lithosphere dominates and inhibits brittle/ductile coupling.

Figure 1 :

 1 Figure 1: Olivine deformation maps (shear stress vs. mean grain size) at a constant strain rate (10 -15 s -1 ) showing the four deformation mechanisms that compete to control the mantle

Figure 2 :

 2 Figure 2: Set-up of numerical modelling showing the initial model geometry, the numerical

Figure 3 :

 3 Figure 3: Model results for the reference model (Model 1): patterns of finite strain, grain size

Figure 4 :

 4 Figure 4: Patterns of finite strain after 140 km of boundary extensional displacement, with

Figure 5 :

 5 Figure 5: Model results for Model 4: patterns of finite strain, grain size (only in the mantle)

Table 2 :

 2 Selected parameters for the four presented numerical results. The reference model is model number 1.

Acknowledgments

We wish to thank the Editor Rob Govers, Dr. Lars Hansen and an anonymous reviewer that provide very detailed and constructive comments on a previous version of our manuscript.