Sebastian Knoedel
email: sebastian.knoedel@googlemail.com

Martin Hachet
email: martin.hachet@inria.fr

Multi-touch RST in 2D and 3D Spaces: Studying the Impact of Directness on User Performance

Keywords: H, 5, 2 [Information Interfaces and Presentation]: User Interfaces-Input devices and strategies

The RST multi-touch technique allows one to simultaneously control Rotations, Scaling, and Translations from multi-touch gestures. We conducted a user study to better understand the impact of directness on user performance for a RST docking task, for both 2D and 3D visualization conditions. This study showed that direct-touch shortens completion times, but indirect interaction improves efficiency and precision, and this is particularly true for 3D visualizations. The study also showed that users' trajectories are comparable for all conditions (2D/3D and direct/indirect). This tends to show that indirect RST control may be valuable for interactive visualization of 3D content. To illustrate this finding, we present a demo application that allows novice users to arrange 3D objects on a 2D virtual plane in an easy and efficient way.

INTRODUCTION

Recent work on multi-touch surfaces has focused on direct interaction, where the input space and the visualization space are colocated. It is commonly argued that direct-touch favors fast and straightforward interaction. On the other hand, indirect multi-touch interaction has some substantial advantages that have been highlighted in the recent literature (eg. [START_REF] Moscovich | Indirect mappings of multi-touch input using one and two hands[END_REF] [START_REF] Schmidt | A comparison of direct and indirect multi-touch input for large surfaces[END_REF]). In particular, indirect multi-touch design is not guided by co-location concerns, which seems particularly interesting for interactive 3D visualization. The goal of this paper is to understand how directness influences user performance in 2D and 3D spatial contexts. In our approach, we do not attempt to use multi-touch gestures to control the 6 degrees-offreedom (DOFs) that can be implied in 3D interaction. Instead, we focus on the well-known RST technique that allows the simultaneous control of Rotations, Scaling, and Translations from two finger inputs [START_REF] Kurtenbach | The design of a gui paradigm based on tablets, two-hands, and transparency[END_REF] [START_REF] Moscovich | Multi-touch interaction[END_REF]. In particular, we wanted to investigate if users were able to interact on 3D-oriented planes as they do with the standard 2D RST technique operated on the screen-plane. To complement the study we developed a demo application, which is illustrated in Figure 1.

A lot of research was conducted in past years to investigate and evaluate multi-touch interaction techniques in 2D spatial context. For example, Latulipe et al. [START_REF] Latulipe | Bimanual and unimanual image alignment: an evaluation of mouse-based techniques[END_REF] showed that symmetric bimanual interaction outperforms the single mouse interaction during image alignment tasks, when users have to manipulate multiple degrees of freedom simultaneously. Moscovich and Hughes [START_REF] Moscovich | Indirect mappings of multi-touch input using one and two hands[END_REF] analyzed two-touch interaction with one and two hands for position and orientation control. They showed that performance of two hands was better when separating interaction tasks. Schmidt et al. [START_REF] Schmidt | A comparison of direct and indirect multi-touch input for large surfaces[END_REF] compared direct and indirect multi-touch interaction on large surfaces. They analyzed completion times of dragging and resizing operations of 2D content during a docking task experiment in direct and indirect conditions. In contrast to their approach, we consider rotation operations as well and conducted further measurements on efficiency and precision to draw conclusions on users coordination. In addition, we are interested in both 2D and 3D visualization setups.

First attempts to create applications that allow the control of 3D content through multiple touch input were made by Zeleznik et al. [START_REF] Zeleznik | Two pointer input for 3d interaction[END_REF]. In recent years, researchers started to combine the multi-point input technology with interactive 3D graphics. Hancock et al. [START_REF] Hancock | Exploring tangible and direct touch interfaces for manipulating 2d and 3d information on a digital table[END_REF] proposed different approaches using direct and indirect multi-touch interaction techniques for tabletop surfaces to manipulate 2D and 3D content. Martinet et al. have studied the effect of DOF separation in a 3D manipulation task on a direct multi-touch display [START_REF] Martinet | The effect of dof separation in 3d manipulation tasks with multi-touch displays[END_REF]. In their work, Reisman et al. [START_REF] Reisman | A screen-space formulation for 2d and 3d direct manipulation[END_REF] proposed co-location constraints for direct-touch object manipulation. The Cubtile has been designed by de la Rivière et al. [START_REF] De La Rivière | Cubtile: a multi-touch cubic interface[END_REF], where users interact with 3D content from an indirect multi-touch box-shaped device.

In our approach, we focus on several aspects that were not addressed by previous work. In particular, we explore how well the RST technique adapts to the 3D visualization space, by looking closely at user performance over direct and indirect conditions. Our overall hypothesis was that indirect interaction may be valuable when visualizing 3D data.

USER STUDY 2.1 Design

We designed an experiment where the task consisted in aligning one source object with a target object. To complete the task, subjects had to control rotations (1 DOF), scaling (1 DOF), and translations (2 DOF), on either the screen plane (2D) or on a 2D plane oriented in 3D space (3D). The source object to manipulate was a blue rectangle in the 2D scenario and a blue box in the 3D scenario (see Figure 2). In all cases both objects were textured with an arrow indicating one distinct orientation. The source object appeared always in the center of the scene at the same position. Initial size and orientation changed randomly. The grey target object was placed with random position, size and orientation around the source object. During the experiment, the viewing angle was set to 45 • to ensure a good visualization of the 3D data. Four red spheres on the upside corners of the target highlight the tolerance areas. They represent the threshold to successfully complete the task. When aligning one corner of the source object correctly with the corresponding corner of the target object in range of the threshold, the color of the sphere turned from red to green. When correctly aligned, the source object had to remain 0.7 seconds to complete the task. We used a within-subject repeated-measures experiment. The independent variables were the visual representation (2D/3D) and the separation condition (direct/indirect). Note that for the 3D indirect condition, the fingers' movements on the touch sensor were directly mapped to cursor movements on the 2D virtual plane that is oriented in 3D space. Participants performed for each of the four conditions 4 repeated tests, which consisted of 8 randomly shuffled trials. Each subject conducted 128 trials in total. The order of presentation of the four conditions was balanced according to a Latin square. Figure 3 24 subjects (8 women, 16 men) participated in the experiment. The subjects used index fingers of both hands to interact with the interface. They were regular computer users, but had no previous experience using a multi-touch interface. Each subject needed approximately 30 minutes to perform the task and to fill out a short questionnaire.

During our study we used a Stantum SMK resistive multi-touch screen as haptic sensor. The display has a size of 15" with a reso-lution of 1280x800 for direct interaction and a touch repositioning accuracy of less than four pixels or less than 0.5 mm. For indirect interaction we used a Dell Notebook display with equal size and identical resolution. All content was rendered with constant 60 frames per second synchronized with the display refresh rate under all four conditions.

We measured efficiency of the trajectories, accuracy, and completion time. For efficiency, we extended the inefficiency measure introduced by Zhai and Milgram [START_REF] Zhai | Quantifying coordination in multiple dof movement and its application to evaluating 6 dof input devices[END_REF]. Inefficiency is the ratio between the distance traveled by users d u and the optimal distance d o between current position and target; defined as

I = d u -d o d o .
Inefficiency can be computed for translation, rotation and scaling, as well as for 2 dimensional spaces (e.g. I T R for translation-rotation) and 3 dimensional spaces (ie. I RST for rotation-scaling-translation). The final I RST measure is the ratio between the user's curve length l RST and the length of the optimal curve l o , I RST = l RST -l o l o . To measure accuracy, we use a target reentry metric. Completion time is given by the duration starting from the time where the subjects touch the multi-touch surface, to the time when the task is successfully completed.

Quantitative Results

We performed two-way ANOVAs on I RST , target reentry, and completion time, which showed the following effects. Figures 4 illus-trates the obtained means.

The subjects performed the task 27.1% more efficiently in 2D than in 3D (2D:1.322, 3D:1.814, F 1,23 = 34.243, p < .0001). Similarly, they were 22% faster (2D:4125ms, 3D:5290ms, F 1,23 = 36.125, p < .0001). On the other hand, no significant difference allows us to confirm an effect on precision.

If we compare direct and indirect modes, we see that indirect interaction is 14.4% more efficient (1.429) than direct interaction (1.708) (F 1,23 = 6.576, p < .05), and subjects are 21.1% more precise (direct:0.47 indirect:0.38, F 1,23 = 5.860, p < .05). This indicates that indirect interaction does not lead to less efficient/precise interaction. Instead users interact even more efficiently and more precisely in indirect mode. Moreover, significant interaction effects show that the influence of the separation mode is bigger in 3D than in 2D (I RST : F 1,23 = 6.439, p < .05, target reentry: F 1,23 = 13.940, p < .05).

Regarding completion time, direct interaction (4115ms) is 22.4% faster than indirect interaction (5300ms) (F 1,23 = 37.858, p < .0001). An interaction effect shows that the influence of distance on completion time is smaller in 3D than in 2D (F 1,23 = 4.345, p < .05).

One important premise of this paper is to investigate user performance in 3D conditions. Therefore, we analyzed interaction effects on precision and efficiency between direct and indirect manipulation modes in 3D. A t-test on I RST reveals a significant effect (t = 3.607, p < .05) showing that in 3D indirect users are 21.6% more efficient (1.595) than in 3D direct (2.034). A similar significant effect can be observed for target reentry (direct:0.481, indirect:0.296, t = 4.735, p < .0001). Regarding completion time, direct interaction remains significantly faster (4820ms) than indirect interaction (5760ms) (t = -4.345, p < .05). Furthermore, our questionnaire shows that 14 out of 24 subjects prefer to manipulate objects from a distance (Wilcoxon Z Test: -2.521, p < .05), when interacting in 3D.

To evaluate the learning effect in the 3D indirect condition, we conducted one-way ANOVAs on I RST and target reentry. Both show no significant difference between the mean values of each consecutive tests during the study (respectively, p > .500, p > .900). Consequently, we cannot show any performance improvement over time. Regarding completion time, our data shows that the completion of the first 8-trials are more time consuming than the following 8-trials tests. No improvement operates starting from the second test. A similar behavior can be observed for the other conditions (2D/3D, direct/indirect).

Discussion

The results presented above allow us to better understand the influence of direct/indirect interaction under 2D and 3D visualization on user performance. In particular, this study highlights that, although direct touch is faster, indirect interaction is more efficient and more precise for RST docking tasks. This can be explained by the fact that dissociating the visual space and the interaction space releases the co-location constraints. We showed that this effect was bigger in 3D than in 2D. Consequently, we believe that indirect interaction may be more suitable for interactive visualization of 3D data. Moreover, the results of the study showed that, with direct interaction, the speed gain is lower in 3D than in 2D.

In our experiment, we set the orientation of the virtual ground to 45 • as a compromise between a good 3D visualization and a convenient configuration for both direct and indirect 3D interaction modes. We can presume that, with orientations of the ground closer to the direction of the viewing direction, the differences of performance between the 3D direct mode and the 3D indirect mode would have been bigger. Indeed, with direct-touch, small finger movements on the screen can project on large distances in the 3D environment, resulting in a loss of precision.

With 3D indirect, the mapping between the multi-touch surface and the virtual ground seems to work well. The analysis of performance evolution over time did not show learning effect. This might indicate that learning of the technique is very fast and takes place during the first trials, or that improvements in performance with time are not large enough to be of statistical or practical significance. This has been confirmed by the subject's comments we obtained from the post-experiment questionnaire. The subjects rated their experiences on a 5 point Likert scale (1 fully disagree and 5 fully agree). They rated the indirect 3D multi-touch interaction technique as easy to use with an average score of 4.4. Moreover, the subjects stated that they were able to solve the task easily with an average score of 4.0. These scores are almost similar to the scores obtained for the 2D direct approach. 2D direct multi-touch interaction currently knows many successes, in particular because users are able to complete interaction tasks in a very easy and convenient way. The experiment we have conducted tends to show that a similar benefit can be obtained under 3D visualizations.

Coordination

Coordination represents the ability of users to control several DOFs at the same time. Since one key component of multi-touch interaction is the possibility to manage several DOFs simultaneously, it is interesting to understand how coordinated users' movements are. The study of the shapes of the obtained inefficiency curves teach us about coordination. We computed the mean curves in the translation-rotation space I T R , the translation-scaling space I T S , the rotation-scaling space I RS , as well as in the rotation-scalingtranslation space I RST . Figure 2.4 (a) shows that the I T R curves, in all conditions, follow the characteristic of the optimal curve (purple). Consequently, translations and rotations seem to be controlled in a coordinated way, independently of the visualization mode and the separation condition.

In contrast, the characteristic of I T S and I RS mean curves (Figure 2.4 (b)) show that users tend to separate scaling from the other DOFs. In both cases, users perform rotation and translation first, and then scale the object. The same behavior is observed in any conditions. This tends to show that users interact with 3D data ata-distance in a similar way they would do with 2D data displayed on a multi-touch screen.

DEMO APPLICATION

We designed a demo application benefiting from RST multi-touch interaction where users manipulate 3D objects on a 2D plane from a distance. This application, dedicated to the general public, can be seen as an extension in 3D of the standard 2D multi-touch image viewer. Such an application can be valuable in many domains such as entertainment and video games, cultural heritage, education, automotive industries, and so on. In this work, we used a scenario where users can move, turn and resize furniture to arrange their living room by sliding their fingers on the multi-touch surface (see Figure 1). Selection. Since visualization and interaction are decoupled, a remote selection mechanism must be introduced. In our approach, users are able to move a standard cursor on the screen -as known from desktop applications-by touching the tactile sensor with a unique finger. When touching the tactile surface with their second finger, users select the pointed object. The application moves from a cursor mode to a manipulation mode. A blue circle appears below the selected object, and the cursor corresponding to the user's fingers are highlighted by way of yellow disks. The application stays in manipulation mode until a single finger movement is detected. Hence, users are able to release their fingers to manipulate comfortably the selected object. To select a new object, they only need to slide one finger on the multi-touch surface to make the cursor appear.

Global camera control. In addition to selection and manipulation, RST multi-touch interaction can be used for global camera control, similar to [START_REF] Zeleznik | Two pointer input for 3d interaction[END_REF]. After having selected the ground surface, users are able to zoom in and out, drag the world and orbit around it in a similar way they do when manipulating a single object. The previously described RST technique is used to transform the global scene, which results in the modification of the camera viewpoint. Finally, a virtual slider, located on the left side of the multi-touch surface, allows changing the orientation of the ground plane (pitch angle).

First user feedback. We asked colleagues to test this demo application to obtain a initial feedback. The general usability appeared to be good. Users had a very short experience with multi-touch surfaces, and they had never experimented our application before. They found the application very easy to use. They were able to position the 3D objects where they wanted, and they felt comfortable using the selection mechanism. They easily switched between the object manipulation mode and the camera control mode. They reported that they particularly appreciated the possibility to change the view to the world in a fast and easy way. After a short time of use, they were able to interact efficiently with the application. This confirmed the results we had from the user study described above. Of course, such an interface does not allow the control of all the DOFs that could be manipulated in the 3D space (eg. roll and vertical translations). It rather provides a simple interface for manipulation and observation of 3D objects resting on a virtual ground. Such a 3D interaction paradigm may be very efficient for many 3D applications where gravity plays an important role.

A direct touch approach could be used for the same application. It may appear as more "natural", and the selection process can be sped up. On the other hand, our study and our experience with the demo application showed that users performance and general usability of the 3D-extended RST technique is good when used in indirect mode. Consequently, similar applications can benefit from the advantages of indirect-touch interaction (e.g. for stereoscopic and large display visualization) and can still be used easily and efficiently by inexperienced users.

CONCLUSION AND DIRECTIONS TO FUTURE WORK

2D applications based on multi-touch input are now well-known and accepted by the general public. Such applications are changing the way we interact with digital content. On the other hand, interaction with 3D content on multi-touch surfaces is not as mature. In particular, the difference between the 2D structure of the input surface and the 3D structure of the visualized data leads to a more complex interaction process. In this paper, we have shown that an adapted version of the RST technique allows users to benefit from multi-touch input for constrained manipulations of 3D objects. By moving cursors on a virtual ground plane, users can move, rotate, or scale 3D objects in a very convenient way.

The user study showed that indirect interaction favors efficiency and precision for a RST docking task. This is particularly true for 3D visualization, because users can separate their finger movements from the 3D visualization. Consequently, we argue that users can benefit from indirect interaction when dealing with 3D data. We illustrated such a configuration with a demo application, which shows the potential of such an interface.

In this paper, we focused in interaction movements on a virtual ground. It would be interesting to see if the same approach could be adapted to other surfaces. For example, we would like to know if the semantic link between the fingers movements on the sensor and the resulting actions on the target surface still works if the latter has an arbitrary orientation (eg. a wall of a room that does not face the camera). One could also imagine different kind of mappings with curved surfaces. These investigations will be part of our future work.

The 3D-oriented RST technique has a limited range of actions. However, we are convinced that such an interface is very valuable compared to the existing mono-point approaches. It allows novice users to well perceive 3D objects by manipulating them easily. In the future, we want to explore new directions for 3D multi-touch interaction. In particular, it would be interesting to investigate good mappings between the multi-touch input and other DOFs in the 3D environment. For example, one can imagine to control nonhomogeneous scaling, or non-canonical rotations. Many other 3D tasks are still to be studied. Multi-touch brings new challenges to research in the scope of 3D user interfaces. It also opens many possibilities for future 3D applications.

Figure 1 :

 1 Figure 1: Indirect control of the RST technique in a 3D furnitures manipulation scenario.

Figure 2 :

 2 Figure 2: Docking task. The source object (blue) has to be aligned with target object (grey transparent), in 2D (a) and in 3D (b).

Figure 3 :

 3 Figure 3: Experimental conditions.

Figure 4 :

 4 Figure 4: The total mean values for inefficiency, target reentry, and completion time for all four conditions.

Figure 5 :

 5 Figure 5: Mismatch curves for all four conditions and optimal curve (purple).

ACKNOWLEDGMENTS

This work was supported by the ANR project ANR-09-CORD-013 InSTInCT -http://anr-instinct.cap-sciences.net.


Céline Labart
email: celine.labart@univ-savoie.fr

Jérôme Lelong
email: jerome.lelong@imag.fr

A Parallel Algorithm for solving BSDEs -Application to the pricing and hedging of American options

Keywords: backward stochastic differential equations, parallel computing, Monte-Carlo methods, non linear PDE, American options, local volatility model

We present a parallel algorithm for solving backward stochastic differential equations (BSDEs in short) which are very useful theoretic tools to deal with many financial problems ranging from option pricing option to risk management. Our algorithm based on Gobet and Labart (2010) exploits the link between BSDEs and non linear partial differential equations (PDEs in short) and hence enables to solve high dimensional non linear PDEs. In this work, we apply it to the pricing and hedging of American options in high dimensional local volatility models, which remains very computationally demanding. We have tested our algorithm up to dimension 10 on a cluster of 512 CPUs and we obtained linear speedups which proves the scalability of our implementation.

Introduction

Pricing and hedging American options with a large number of underlying assets is a challenging financial issue. On a single processor system, it can require several hours of computation in high dimensions.

Recent advances in parallel computing hardware such as multi-core processors, clusters and GPUs are then of high interest for the finance community. For a couple of years, some research teams have been tackling the parallelization of numerical algorithms for option pricing. [START_REF] Thulasiram | Performance evaluation of parallel algorithms for pricing multidimensional[END_REF] developed a parallel algorithm using MPI for pricing a class of multidimensional financial derivatives using a binomial lattice approach. [START_REF] Huang | Parallel Algorithm for pricing american Asian Options with Multi-Dimensional Assets[END_REF] presented algorithms for pricing American style Asian options using a binomial tree method. Concerning the parallelization of Monte-Carlo methods for pricing multi-dimensional Bermudan/American options, the literature is quite rare. We refer to [START_REF] Toke | Monte Carlo Valuation of Multidimensional American Options Through Grid Computing[END_REF] and to Dung [START_REF] Doan | Parallel pricing algorithms for multimensional bermudan/american options using Monte Carlo methods[END_REF]. Both papers propose a parallelization through grid computing of the [START_REF] Ibáñez | Valuation by simulation of american options through computation of the optimal exercise frontier[END_REF] algorithm, which computes the optimal exercise boundary. A GPGPU approach based on quantization techniques has recently been developed by [START_REF] Pagès | GPGPUs in computational finance: Massive parallel computing for American style options[END_REF].

Our approach is based on solving Backward Stochastic Differential Equations (BSDEs in short). As explained in the seminal paper by El [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF], pricing and hedging European options in local volatility models boil down to solving standard BSDEs. From El Karoui et al. (1997a), we know that the price of an American option is also linked to a particular class of BSDEs called reflected BSDEs. Several sequential algorithms to solve BS-DEs can be found in the literature. [START_REF] Ma | Solving forward backward stochastic differential equations explicitly-a four step scheme[END_REF] adopted a PDE approach, whereas [START_REF] Bouchard | Discrete time approximation and Monte Carlo simulation of backward stochastic differential equations[END_REF] and [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF] used a Monte-Carlo approach based on the dynamic programming equation. The Monte-Carlo approach was also investigated by [START_REF] Bally | Error analysis of the optimal quantization algorithm for obstacle problems[END_REF] and [START_REF] Delarue | A forward-backward stochastic algorithm for quasi-linear PDEs[END_REF] who applied quantization techniques to solve the dynamic programming equation. Our approach is based on the algorithm developed by [START_REF] Gobet | Solving BSDE with adaptive control variate[END_REF] which combines Picard's iterations and an adaptive control variate. It enables to solve standard BSDEs, ie, to get the price and delta of European options in a local volatility model. Compared to the algorithms based on the dynamic programming equation, ours provides regular solutions in time and space (which is coherent with the regularity of the option price and delta). To apply it to the pricing and hedging of American options, we use a technique introduced by El Karoui et al. (1997a), which consists in approximating a reflected BSDE by a sequence of standard BSDEs with penalisation.

The paper is organized as follows. In section 2, we briefly recall the link between BSDEs and PDEs which is the heart of our algorithm. In section 3, we describe the algorithm and in Section 4 we explain how the parallelization has been carried out. Section 5 describes how American options can be priced using BSDEs. Finally, in Section 6, we conclude the paper by some numerical tests of our parallel algorithm for pricing and hedging European and American basket options in dimension up to 10.

Definitions and Notations

• Let C k,l b be the set of continuously differentiable functions φ : (t, x) ∈ [0, T ] × R d with continuous and uniformly bounded derivatives w.r.t. t (resp. w.r.t. x) up to order k (resp. up to order l).

• C k p denotes the set of C k-1 functions with piecewise continuous k th order derivative.

• For α ∈]0, 1], C k+α is the set of C k functions whose k th order derivative is Hölder continuous with order α.

BSDEs 2.1 General results on standard BSDEs

Let (Ω, F, P) be a given probability space on which is defined a q-dimensional standard Brownian motion W , whose natural filtration, augmented with P-null sets, is denoted (F t ) 0≤t≤T (T is a fixed terminal time). We denote (Y, Z) the solution of the following backward stochastic differential equation (BSDE) with fixed terminal time T

-dY t = f (t, X t , Y t , Z t )dt -Z t dW t , Y T = Φ(X T ), (2.1)
where f : [0, T ] × R d × R × R q → R, Φ : R d → R and X is the R d -valued process solution of

X t = x + t 0 b(s, X s )ds + t 0 σ(s, X s )dW s , (2.2) with b : [0, T ] × R d → R d and σ : [0, T ] × R d → R d×q .
From now on, we assume the following Hypothesis, which ensures the existence and uniqueness of the solution to Equations (2.1)-(2.2).

Hypothesis 1

• The driver f is a bounded Lipschitz continuous function, ie, for all

(t 1 , x 1 , y 1 , z 1 ), (t 2 , x 2 , y 2 , z 2 ) ∈ [0, T ] × R d × R × R, ∃ L f > 0, |f (t 1 , x 1 , y 1 , z 1 ) -f (t 2 , x 2 , y 2 , z 2 )| ≤ L f (|t 1 -t 2 | + |x 1 -x 2 | + |y 1 -y 2 | + |z 1 -z 2 |).
• σ is uniformly elliptic on [0, T ] × R d , ie, there exist two positive constants σ 0 , σ 1 s.t. for any ξ ∈ R d and any

(t, x) ∈ [0, T ] × R d σ 0 |ξ| 2 ≤ d i,j=1 [σσ * ] i,j (t, x)ξ i ξ j ≤ σ 1 |ξ| 2 . • Φ is bounded in C 2+α , α ∈]0, 1]. • b and σ are in C 1,3 b and ∂ t σ is in C 0,1 b .

Link with semilinear PDEs

Let us also recall the link between BSDEs and semilinear PDEs. Although the relation is the keystone of our algorithm, as explained in Section 3, we do not develop it and refer to [START_REF] Pardoux | Backward Stochastic Differential Equations and Quasilinear Parabolic Partial Differential Equations[END_REF] or El [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF] for more details.

According to Pardoux and Peng (1992, Theorem 3.1), we can link the solution (Y, Z) of the BSDE (2.1) to the solution u of the following PDE:

∂ t u(t, x) + Lu(t, x) + f (t, x, u(t, x), (∂ x uσ)(t, x)) = 0, u(T, x) = Φ(x), (2.3)
where L is defined by

L (t,x) u(t, x) = 1 2 i,j [σσ * ] ij (t, x)∂ 2 x i x j u(t, x) + i b i (t, x)∂ x i u(t, x).
Theorem 1 [START_REF] Delarue | A forward-backward stochastic algorithm for quasi-linear PDEs[END_REF], Theorem 2.1). Under Hypothesis 1, the solution

u of PDE (2.3) belongs to C 1,2 b . Moreover, the solution (Y t , Z t ) 0≤t≤T of (2.1) satisfies ∀t ∈ [0, T ], (Y t , Z t ) = (u(t, X t ), ∂ x u(t, X t )σ(t, X t )).
(2.4)

3 Presentation of the Algorithm

Description

We present the algorithm introduced by [START_REF] Gobet | Solving BSDE with adaptive control variate[END_REF] to solve standard BSDEs.

It is based on Picard's iterations combined with an adaptive Monte-Carlo method. We recall that we aim at numerically solving BSDE (2.1), which is equivalent to solving the semilinear PDE (2.3). The current algorithm provides an approximation of the solution of this PDE. Let u k denote the approximation of the solution u of (2.3) at step k. If we are able to compute an explicit solution of (2.2), the approximation of (Y, Z) at step k follows from (2.4):

(Y k t , Z k t ) = (u k (t, X t ), ∂ x u k (t, X t )σ(t, X t ))
, for all t ∈ [0, T ]. Otherwise, we introduce X N the approximation of X obtained with a N -time step Euler scheme:

∀s ∈ [0, T ], dX N s = b(ϕ N (s), X N ϕ N (s) )ds + σ(ϕ N (s), X N ϕ N (s) )dW s , (3.1)
where ϕ N (s) = sup{t j : t j ≤ s} is the largest discretization time not greater than s and

{0 = t 0 < t 1 < • • • < t N = T } is a regular subdivision of the interval [0, T ]. Then, we write (Y k t , Z k t ) = (u k (t, X N t ), ∂ x u k (t, X N t )σ(t, X N t )), for all t ∈ [0, T ].
It remains to explain how to build the approximation (u k ) k of u. The basic idea is the following:

u k+1 = u k + Monte-Carlo evaluations of the error(u -u k ).
Combining Itô's formula applied to u(s, X s ) and u k (s, X N s ) between t and T and the semilinear PDE (2.3) satisfied by u, we get that the correction term c k is given by

c k (t, x) = (u -u k )(t, x) = E Ψ (t, x, f u , Φ, W ) -Ψ N t, x, -(∂ t + L N )u k , u k (T, .), W |G k where • L N u(s, X N s ) = 1 2 i,j [σσ * ] ij (ϕ(s), X N ϕ(s) )∂ 2 x i x j u(s, X N s )+ i b i (ϕ(s), X N ϕ(s) )∂ x i u(s, X N s ). • f v : [0, T ] × R d → R denotes the following function f v (t, x) = f (t, x, v(t, x), (∂ x vσ)(t, x)),
where f is the driver of BSDE (2.1), σ is the diffusion coefficient of the SDE satisfied by X and v :

[0, T ] × R d → R is C 1 w.r.t.
to its second argument.

• Ψ and Ψ N denote

Ψ(s, y, g 1 , g 2 , W ) = T s g 1 (r, X s,y r (W ))dr + g 2 (X s,y T (W )), Ψ N (s, y, g 1 , g 2 , W ) = T s g 1 (r, X N,s,y r (W ))dr + g 2 (X N,s,y T (W )),
where X s,y (resp. X N,s,y ) denotes the diffusion process solving (2.2) and starting from y at time s (resp. its approximation using an Euler scheme with N time steps), and W denotes the standard Brownian motion appearing in (2.2) and used to simulate X N , as given in (3.1).

• G k is the σ-algebra generated by the set of all random variables used to build u k . In the above formula of c k , we compute the expectation w.r.t. the law of X and X N and not w.r.t. the law of u k , which is G k measurable. (See Definition 2 for a rigorous definition of G k ).

Note that Ψ and Ψ N can actually be written as expectations by introducing a random variable U uniformly distributed on [0, 1].

Ψ(s, y, g 1 , g 2 , W ) = E U (T -s)g 1 (s + (T -s)U, X s,y s+(T -s)U (W )) + g 2 (X s,y T (W )) , Ψ N (s, y, g 1 , g 2 , W ) = E U (T -s)g 1 (s + (T -s)U, X N,s,y s+(T -s)U (W )) + g 2 (X N,s,y T (W )) .
In the following, let ψ N (s, y, g 1 , g 2 , W, U ) denote

ψ N (s, y, g 1 , g 2 , W, U ) = (T -s)g 1 (s + (T -s)U, X N,s,y s+(T -s)U (W )) + g 2 (X N,s,y T (W )) (3.2) such that Ψ N (s, y, g 1 , g 2 , W ) = E U [ψ N (s, y, g 1 , g 2 , W, U )].
From a practical point of view, the PDE (2.3) is solved on [0, T ] × D where D ⊂ R d such that sup 0≤t≤T |X t | ∈ D with a probability very close to 1.

Algorithm 1 We begin with u 0 ≡ 0. Assume that an approximated solution u k of class C 1,2 is built at step k -1. Here are the different steps to compute u k+1 .

• Pick at random n points (t k i , x k i ) 1≤i≤n uniformly distributed over [0, T ] × D. • Evaluate the Monte-Carlo correction c k at step k at the points (t k i , x k i ) 1≤i≤n using M independent simulations

c k (t k i , x k i ) = 1 M M m=1 ψ N t k i , x k i , f u k + (∂ t + L N )u k , Φ -u k , W m,k,i , U m,k,i . • Compute the vector (u k (t k i , x k i )) 1≤i≤n . Now, we know the vector (u k + c k )(t k i , x k i )) 1≤i≤n .
From these values, we extrapolate the function

u k+1 = u k + c k on [0, T ] × D. u k+1 (t, x) = P k (u k + c k )(t, x), for (t, x) ∈ [0, T ] × D, (3.3)
where P k is a deterministic operator, which only uses the values of the function at the points (t k i , x k i ) 1≤i≤n to approximate the function on the whole domain [0, T ] × D. The choice of the operator P k is discussed in Section 3.2. Since c k is computed using Monte-Carlo simulations instead of a true expectation, the values (c k (t k i , x k i )) 1≤i≤n are random variables. Therefore, u k+1 is a random function depending on the random variables needed to compute

u k and W m,k,i , U m,k,i , 1 ≤ m ≤ M, 1 ≤ i ≤ n.
In view of this comment, the σ-algebra G k has to be redefined to take into account this new source of randomness.

Definition 2 (Definition of the σ-algebra G k ). Let G k+1 define the σ-algebra generated by the set of all random variables used to build u k+1 . Using (3.3) yields

G k+1 = G k ∨ σ(A k , S k ),
where A k is the set of random points used at step k to build the estimator P k (see below),

S k = {W m,k,i , U m,k,i , 1 ≤ m ≤ M, 1 ≤ i ≤ n},
is the set of independent Brownian motions used to simulate the paths X m,k,N (x k i ), and G k is the σ-algebra generated by the set of all random variables used to build u k .

Choice of the operator

The most delicate part of the algorithm is how to extrapolate a function h and its derivatives when only knowing its values at n points (t i , x i ) i=1,...,n ∈ [0, T ] × D.

A kernel operator

In the first version of Algorithm 1 presented in [START_REF] Gobet | Solving BSDE with adaptive control variate[END_REF], a function h was extrapolated from the values computed on the grid by using a kernel operator of the form

h(t, x) = n i=1 u(t i , x i )K t (t -t i )K x (x -x i ),
where K t is a one dimensional kernel whereas K x is a product of d one dimensional kernels. Hence, evaluating the function h at a given point (t, x) requires O(n × d) computations.

The convergence result established by Gobet and Labart (2010, Theorem 5.1) is based on the properties of the operator presented in Gobet and Labart (2010, Section 4). Using the linearity and the boundedness of the operator, they managed to prove that the errors v -P k v and ∂ x v -∂ x (P k v) are bounded, which is a key step in proving the convergence of the algorithm. At the end of their paper, they present an operator based on kernel estimators satisfying the assumptions required to prove the convergence of the algorithm.

An extrapolating operator

The numerical properties of kernel operators are very sensitive to the choice of their window parameters which is quite hard to tune for each new problem. Hence, we have tried to use an other solution. Basically, we have borrowed the solution proposed by [START_REF] Longstaff | Valuing American options by simulation : A simple least-square approach[END_REF] which consists in extrapolating a function by solving a least square problem defined by the projection of the original function on a countable set of functions. Assume we know the values (y i ) i=1,...,n of a function h at the points (t i , x i ) i=1,...,n , the function h can be extrapolated by computing

α = arg min α∈R p n i=1 y i - p l=1 α l B l (t i , x i ) 2 , (3.4)
where (B l ) l=1,...,p are some real valued functions defined on [0, T ]×D. Once α is computed, we set ĥ(t, x) = p l=1 α l B l (t, x). For the implementation, we have chosen the (B l ) l=1,...,p as a free family of multivariate polynomials. For such a choice, ĥ is known to converge uniformly to h when p goes to infinity if D is a compact set and h is continuous on [0, T ] × D. Our algorithm also requires to compute the first and second derivatives of h which are approximated by the first and second derivatives of ĥ. Although the idea of approximating the derivatives of a function by the derivatives of its approximation is not theoretically well justified, it is proved to be very efficient in practice. We refer to [START_REF] Wang | Pricing and hedging american-style options: a simple simulationbased approach[END_REF] for an application of this principle to the computations of the Greeks for American options.

Practical determination of the vector α In this part, we use the notation d ′ = d + 1. It is quite easy to see from Equation (3.4) that α is the solution of a linear system. The value α is a critical point of the criteria to be minimized in Equation (3.4) and the vector α solves

p l=1 α l n i=1 B l (t i , x i )B j (t i , x i ) = n i=1 y i B j (t i , x i ) for j = 1, . . . , p Aα = n i=1 y i B(t i , x i ) (3.5)
where the p × p matrix A = ( n i=1 B l (t i , x i )B j (t i , x i )) l,j=1,...,p and the vector B = (B 1 , . . . , B p ) * .

The matrix A is symmetric and positive definite but often ill-conditioned, so we cannot rely on the Cholesky factorization to solve the linear system but instead we have to use some more elaborate techniques such as a QR factorization with pivoting or a singular value decomposition approach which can better handle an almost rank deficient matrix. In our implementation of Algorithm 1, we rely on the routine dgelsy from Lapack [START_REF] Anderson | LAPACK Users' Guide[END_REF], which solves a linear system in the least square sense by using some QR decomposition with pivoting combined with some orthogonalization techniques. Fortunately, the ill-conditioning of the matrix A is not fate; we can improve the situation by centering and normalizing the polynomials (B l ) l such that the domain [0, T ] × D is actually mapped to [-1, 1] d ′ . This reduction improves the numerical behaviour of the chaos decomposition by a great deal.

The construction of the matrix A has a complexity of O(p 2 nd ′ ). The computation of α (Equation 4.3) requires to solve a linear system of size p × p which requires O(p 3 ) operations. The overall complexity for computing α is then

O(p 3 + np 2 d ′ ).
Choice of the (B l ) l . The function u k we want to extrapolate at each step of the algorithm is proved to be quite regular (at least C 1,2 ), so using multivariate polynomials for the B l should provide a satisfactory approximation. Actually, we used polynomials with d ′ variates, which are built using tensor products of univariate polynomials and if one wants the vector space Vect{B l , l = 1, . . . , p} to be the space of d ′ -variate polynomials with global degree less or equal than η, then p has to be equal to the binomial coefficient d ′ +η η . For instance, for η = 3 and d ′ = 6 we find p = 84. This little example shows that p cannot be fixed by specifying the maximum global degree of the polynomials B l without leading to an explosion of the computational cost, we therefore had to find an other approach. To cope with the curse of dimensionality, we studied different strategies for truncating polynomial chaos expansions. We refer the reader to Chapter 2 of [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF] for a detailed review on the topic. From a computational point of view, we could not afford the use of adaptive sparse polynomial families because the construction of the family is inevitably sequential and it would have been detrimental for the speed-up of our parallel algorithm. Therefore, we decided to use sparse polynomial chaos approximation based on an hyperbolic set of indices as introduced by [START_REF] Blatman | Anisotropic parcimonious polynomial chaos expansions based on the sparsity-f-effects principle[END_REF].

A canonical polynomial with d ′ variates can be defined by a multi-index ν ∈ N d ′ν i being the degree of the polynomial with respect the variate i. Truncating a polynomial chaos expansion by keeping only the polynomials with total degree not greater than η corresponds to the set of multi-indices:

{ν ∈ N d ′ : d ′ i=1 ν i ≤ η}.
The idea of hyperbolic sets of indices is to consider the pseudo q-norm of the multi-index ν with q ≤ 1

   ν ∈ N d ′ : d ′ i=1 ν q i 1/q ≤ η    . (3.6)
Note that choosing q = 1 gives the full family of polynomials with total degree not greater than η. The effect of introducing this pseudo-norm is to favor low-order interactions.

Parallel approach

In this part, we present a parallel version of Algorithm 1, which is far from being embarrassingly parallel as a crude Monte-Carlo algorithm. We explain the difficulties encountered when parallelizing the algorithm and how we solved them.

Detailed presentation of the algorithm

Here are the notations we use in the algorithm.

• u k = (u k (t k i , x k i )) 1≤i≤n ∈ R n • c k = (c k (t k i , x k i )) 1≤i≤n ∈ R n
• n: number of points of the grid

• K it : number of iterations of the algorithm

• M : number of Monte-Carlo samples

• N : number of time steps used for the discretization of X

• p: number of functions B l used in the extrapolating operator. This is not a parameter of the algorithm on its own as it is determined by fixing η and q (the maximum total degree and the parameter of the hyperbolic multi-index set) but the parameter p is of great interest when studying the complexity of the algorithm.

• (B l ) 1≤l≤p is a family of multivariate polynomials used for extrapolating functions from a finite number of values.

• α k ∈ R p is the vector of the weights of the chaos decomposition of u k .

• d ′ = d + 1 is the number of variates of the polynomials B l .

Algorithm 1 Iterative algorithm

1: u 0 ≡ 0, α 0 ≡ 0. 2: for k = 0 : K it -1 do 3:
Pick at random n points (t k i , x k i ) 1≤i≤n .

4:

for i = 1 : n do 5:

for m = 1 : M do 6:

Let W be a Brownian motion with values in R d discretized on a time grid with N time steps.

7:

Let U ∼ U [0,1] . 8: Compute a i,k m = ψ N t k i , x k i , f u k + (∂ t + L N )u k , Φ -u k , W , U . / * We recall that u k (t, x) = p l=1 α k l B l (t, x) * / 9:
end for

c k i = 1 M M m=1 a i,k m (4.1)
u k i = p l=1 α k l B l (t k i , x k i ) (4.2) 10:
end for 11:

Compute

α k+1 = arg min α∈R p n i=1 (u k i + c k i ) - p l=1 α l B l (t k i , x k i ) 2 .
(4.3)

12: end for

Complexity of the algorithm

In this section, we study in details the different parts of Algorithm 1 to determine their complexities. Before diving into the algorithm, we would like to briefly look at the evaluations of the function u k and its derivatives. We recall that

u k (t, x) = p l=1 α k l B l (t, x)
where the B l (t, x) are of the form t β l,0 d i=1 x i β l,i and the β l,i are some integers. Then the computational time for the evaluation of u k (t, x) is proportional to p × d ′ . The first and second derivatives of u k write

∇ x u k (t, x) = p l=1 α k l ∇ x B l (t, x), ∇ 2 x u k (t, x) = p l=1 α k l ∇ 2 x B l (t, x),
and the evaluation of ∇ x B l (t, x) (resp. ∇ 2 x B l (t, x)) has a computational cost proportional to d 2 (resp. d 3 ).

• The computation (at line 6) of the discretization of the d-dimensional Brownian motion with N time steps requires O(N d) computations.

• The computation of each a k,i m (line 8) requires the evaluation of the function u k and its first and second derivatives which has a cost O(pd 3 ). Then, the computation of c k i for given i and k has a complexity of O(M pd 3 ).

• The computation of α (Equation 4.3) requires O(p 3 + np 2 d) operations as explained in Section 3.2.2.

The overall complexity of Algorithm 1 is

O(K it nM (pd 3 + dN ) + K it (p 2 nd + p 3 )).
To parallelize an algorithm, the first idea coming to mind is to find loops with independent iterations which could be spread out on different processors with a minimum of communications. In that respect, an embarrassingly parallel example is the well-known Monte-Carlo algorithm. Unfortunately, Algorithm 1 is far from being so simple. The iterations of the outer loop (line 2) are linked from one step to the following, consequently there is no hope parallelizing this loop. On the contrary, the iterations over i (loop line 4) are independent as are the ones over m (loop line 5), so we have at hand two candidates to implement parallelizing. We could even think of a 2 stage parallelism : first parallelizing the loop over i over a small set of processors and inside this first level parallelizing the loop over m. Actually, M is not large enough for the parallelization of the loop over m to be efficient (see Section 3.2). It turns out to be far more efficient to parallelize the loop over i as each iteration of the loop requires a significant amount of work.

Description of the parallel part

As we have just explained, we have decided to parallelize the loop over i (line 4 in Algorithm 1). We have used a Robbin Hood approach. In the following, we assume to have P + 1 processors at hand with n > P . We use the following master slave scheme:

1. Send to each of the P slave processors the solution α k computed at the previous step of the algorithm.

2. Spread the first P points (t k i , x k i ) 1≤i≤P to the P slave processors and assign each of them the computation of the corresponding c k i .

3. As soon as a processor finishes its computation, it sends its result back to the master which in turn sends it a new point (t k i , x k i ) at which evaluating ψ N to compute c k i and this process goes on until all the (c k i ) i=1,...,n have been computed.

At the end of this process, the master knows c k , which corresponds to the approximation of the correction at the random points (t k i , x k i ) 1≤i≤n . From these values and the vector u k , the master computes α k+1 and the algorithm can go through a new iteration over k.

What is the best way to send the data to each slave process?

• Before starting any computations, assign to each process a block of iterations over i and send the corresponding data all at once. This way, just one connection has to be initialised which is faster. But the time spent by the master to take care of its slave is longer which implies that at the beginning the slave process will remain longer unemployed. When applying such a strategy, we implicitly assume that all the iterations have the same computational cost.

• Send data iteration by iteration. The latency at the beginning is smaller than in the block strategy and performs better when all iterations do not have the same computational cost.

Considering the wide range of data to be sent and the intensive use of elaborate structures, the most natural way to pass these objects was to rely on the packing mechanism of MPI.

Moreover, the library we are using in the code (see Section 4.5) already has a MPI binding which makes the manipulation of the different objects all the more easy. The use of packing enabled to reduce the number of communications between the master process and a slave process to just one communication at the beginning of each iteration of the loop over i (line 4 of Algorithm 1).

Random numbers in a parallel environment

One of the basic problem when solving a probabilistic problem in parallel computing is the generation of random numbers. Random number generators are usually devised for sequential use only and special care should be taken in parallel environments to ensure that the sequences of random numbers generated on each processor are independent. We would like to have minimal communications between the different random number generators, ideally after the initialisation process, each generator should live independently of the others. Several strategies exist for that.

1. Newbies in parallel computing might be tempted to take any random number generator and fix a different seed on each processor. Although this naive strategy often leads to satisfactory results on toy examples, it can induce an important bias and lead to detrimental results.

2. The first reasonable approach is to split a sequence of random numbers across several processors. To efficiently implement this strategy, the generator must have some splitting facilities such that there is no need to draw all the samples prior to any computations. We refer to L' Ecuyer and Côté (1991); L 'Ecuyer et al. (2002) for a presentation of a generator with splitting facilities. To efficiently split the sequence, one should know in advance the number of samples needed by each processor or at least an upper bound of it. To encounter this problem, the splitting could be made in substreams by jumping ahead of P steps at each call to the random procedure if P is the number of processors involved. This way, each processor uses a sub-sequence of the initial random number sequence rather than a contiguous part of it. However, as noted by [START_REF] Entacher | Parallel random number generation: Long-range correlations among multiple processors[END_REF], long range correlations in the original sequence can become short range correlations between different processors when using substreams.

Actually, the best way to implement splitting is to use a generator with a huge period such as the Mersenne Twister (its period is 2 19937 -1) and divide the period by a million or so if we think we will not need more than a million independent substreams. Doing so, we come up with substreams which still have an impressive length, in the case of the Mersenne Twister each substream is still about 2 19917 long.

3. A totally different approach is to find generators which can be easily parametrised and to compute sets of parameters ensuring the statistical independence of the related generators. Several generators offer such a facility such as the ones included in the SPRNG package (see [START_REF] Mascagni | Some methods of parallel pseudorandom number generation[END_REF] for a detailed presentation of the generators implemented in this package) or the dynamically created Mersenne Twister (DCMT in short), see [START_REF] Matsumoto | chapter Dynamic Creation of Pseudorandom Number Generator[END_REF].

For our experiments, we have decided to use the DCMT. This generator has a sufficiently long period (2 521 for the version we used) and we can create at most 2 16 = 65536 independent generators with this period which is definitely enough for our needs. Moreover, the dynamic creation of the generators follows a deterministic process (if we use the same seeds) which makes it reproducible. The drawback of the DCMT is that its initialization process might be quite lengthy, actually the CPU time needed to create a new generator is not bounded. We give in Figure 1

The library used for the implementation

Our code has been implemented in C using the PNL library (see Lelong (2007Lelong ( -2011))). This is a scientific library available under the Lesser General Public Licence and it offers various facilities for implementing mathematics and more recently some MPI bindings have been added to easily manipulate the different objects available in PNL such as vectors and matrices.

In our problem, we needed to manipulate matrices and vectors and pass them from the master process to the slave processes and decided to use the packing facilities offered by PNL through its MPI binding. The technical part was not only message passing but also random number generation as we already mentioned above and PNL offers many functions to generate random vectors or matrices using several random number generators among which the DCMT. Besides message passing, the algorithm also requires many other facilities such as multivariate polynomial chaos decomposition which is part of the library. For the moment, three families of polynomials (Canonical, Hermite and Tchebichev polynomials) are implemented along with very efficient mechanism to compute their first and second derivatives. The implementation tries to make the most of code factorization to avoid recomputing common quantities several times. The polynomial chaos decomposition toolbox is quite flexible and offers a reduction facility such as described in Section 3.2.2 which is completely transparent from the user's side. To face the curse of dimensionality, we used sparse polynomial families based on an hyperbolic set of indices.

Pricing and Hedging American options in local volatility models

In this Section, we present a method to price and hedge American options by using Algorithm 1, which solves standard BSDEs.

Framework

Let (Ω, F, P) be a probability space and (W t ) t≥0 a standard Brownian motion with values in R d . We denote by (F t ) t≥0 the P-completion of the natural filtration of (W t ) t≥0 . Consider a financial market with d risky assets, with prices S 1 t , • • • , S d t at time t, and let X t be the d-dimensional vector of log-returns X i t = log S i t , for i = 1, • • • , d. We assume that (X t ) satisfies the following stochastic differential equation:

dX i t =   r(t) -δ i (t) - 1 2 d j=1 σ ij 2 (t, e Xt )   dt + d j=1 σ ij (t, e Xt )dW j t , i = 1, • • • , d (5.1)
on a finite interval [0, T ], where T is the maturity of the option. We denote by X t,x s a continuous version of the flow of the stochastic differential Equation (5.1). X t,x t = x almost surely.

In the following, we assume Hypothesis 2

1. r : [0, T ] -→ R is a C 1 b function. δ : [0, T ] -→ R d is a C 1 b function. 2. σ : [0, T ] × R d -→ R d×d is a C 1,2 b function.
3. σ satisfies the following coercivity property:

∃ε > 0 ∀(t, x) ∈ [0, T ] × R d , ∀ξ ∈ R d 1≤i,j≤d [σσ * ] i,j (t, x)ξ i ξ j ≥ ε d i=1 ξ 2 i
We are interested in computing the price of an American option with payoff Φ(X t ), where Φ :

R d -→ R + is a continuous function (in case of a put option, Φ(x) = (K -1 d (e x 1 + • • • + e x d ) + ).
From [START_REF] Jaillet | Variational inequalities and the pricing of American options[END_REF], we know that if Hypothesis 3 Φ is continuous and satisfies |Φ(x)| ≤ M e M |x| for some M > 0, the price at time t of the American option with payoff Φ is given by

V (t, X t ) = esssup τ ∈T t,T E e -τ t r(s)ds Φ(X τ )|F t ,
where T t,T is the set of all stopping times with values in [t, T ] and

V : [0, T ] × R d -→ R + defined by V (t, x) = sup τ ∈T t,T E e -τ t r(s)ds Φ(X t,x τ
) . We also introduce the amount ∆(t, X t ) involved in the asset at time t to hedge the American option. The function ∆ is given by ∆(t, x) = ∇ x V (t, x).

In order to link American option prices to BSDEs, we need three steps:

1. writing the price of an American option as a solution of a variational inequality (see Section 5.2)

2. linking the solution of a variational inequality to the solution of a reflected BSDE (see Section 5.3) 3. approximating the solution of a RBSDE by a sequence of standard BSDEs (see Section 5.4)

We refer to [START_REF] Jaillet | Variational inequalities and the pricing of American options[END_REF] for the first step and to El Karoui et al. (1997a) for the second and third steps.

American option and Variational Inequality

First, we recall the variational inequality satisfied by V . We refer to Jaillet et al. (1990, Theorem 3.1) for more details. Under Hypotheses 2 and 3, V solves the following parabolic PDE max(Φ(x)u(t, x), ∂ t u(t, x) + Au(t, x)r(t)u(t, x)) = 0, u(T, x) = Φ(x).

(5.2)

where

Au(t, x) = d i=1 (r(t)-δ i (t)-1 2 d j=1 σ 2 ij (t, e x ))∂ x i u(t, x)+ 1 2 1≤i,j≤d [σσ * ] ij (t, e x )∂ 2 x i x j u(t, x
) is the generator of X.

Variational Inequality and Reflected BSDEs

This section is based on El Karoui et al. (1997a, Theorem 8.5). We assume Hypothesis 4 Φ is continuous and has at most a polynomial growth (ie, ∃C, p > 0 s.t.

|Φ(x)| ≤ C(1 + |x| p )).
Let us consider the following reflected BSDE

     -dY t = -r(t)Y t dt -Z t dW t + dH t , Y T = Φ(X T ), Y t ≥ Φ(X t ), ∀t, T 0 (Y t -Φ(X t ))dH t = 0.
(5.3) Then, under Hypotheses 2 and 4, u(t, x) = Y t,x t is a viscosity solution of the obstacle problem (5.2), where Y t,x t is the value at time t of the solution of (5.3) on [t, T ] where the superscripts t, x mean that X starts from x at time t. We also have (∇ x u(t, x)) * σ(t, x) = Z t,x t ( * means transpose). Then, we get that the price V and the delta ∆ of the option are given by

V (t, X t ) = Y t , ∆(t, X t ) = (Z t σ(t, X t ) -1 ) * .

Approximation of a RBSDE by a sequence of standard BSDEs

We present a way of approximating a RBSDE by a sequence of standard BSDEs. The idea was introduced by El Karoui et al. (1997a, Section 6) for proving the existence of a solution to RBSDE by turning the constraint Y t ≥ Φ(X t ) into a penalisation. Let us consider the following sequence of BSDEs indexed by i (5.4) whose solutions are denoted (Y i , Z i ). We define

Y i t = Φ(X T ) - T t r(s)Y i s ds + i T t (Y i s -Φ(X s )) -ds - T t Z i s dW s ,
H i t = i T t (Y i s -Φ(X s )) -ds.
Under Hypotheses 2 and 4, the sequence (Y i , Z i , H i ) i converges to the solution (Y, Z, H) of the RBSDE (5.3), when i goes to infinity. Moreover, Y i converges increasingly to Y . The term H i is often called a penalisation. From a practical point, there is no use solving such a sequence of BSDEs, because we can directly apply our algorithm to solve Equation (5.4) for a given i. Therefore, we actually consider the following penalized BSDE

Y t = Φ(X T ) - T t r(s)Y s ds + ω T t (Y s -Φ(X s )) -ds - T t Z s dW s , (5.5)
where ω ≥ 0 is penalization weight. In practice, the magnitude of ω must remain reasonably small as it appears as a contraction constant when studying the speed of convergence of our algorithm. Hence, the larger ω is, the slower our algorithm converges. So, a trade-off has to be found between the convergence of our algorithm to the solution Y of (5.5) and the accuracy of the approximation of the American option price by Y .

European options

Let us consider a European option with payoff Φ(X T ), where X follows (5.1). We denote by V the option price and by ∆ the hedging strategy associated to the option. From El [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF], we know that the couple (V, ∆) satisfies

-dV t = -r(t)V t dt -∆ * t σ(t, X t )dW t , V T = Φ(X T ).
Then, (V, ∆) is solution of a standard BSDE. This corresponds to the particular case ω = 0 of (5.5), Y corresponding to V and Z to ∆ * t σ(t, X t ).

6 Numerical results and performance

The cluster

All our performance tests have been carried out on a 256-PC cluster from SUPELEC Metz. Each node is a dual core processor : INTEL Xeon-3075 2.66 GHz with a front side bus at 1333Mhz. The two cores of each node share 4GB of RAM and all the nodes are interconnected using a Gigabit Ethernet network. In none of the experiments, did we make the most of the dual core architecture since our code is one threaded. Hence, in our implementation a dual core processor is actually seen as two single core processors. The accuracy tests have been achieved using the facilities offered by the University of Savoie computing center MUST.

Black-Scholes' framework

We consider a d-dimensional Black-Scholes model in which the dynamics under the risk neutral measure of each asset S i is supposed to be given by

dS i t = S i t ((r -δ i )dt + σ i dW i t ) S 0 = (S 1 0 , . . . , S d 0 ) (6.1)
where W = (W 1 , . . . , W d ). Each component W i is a standard Brownian motion. For the numerical experiments, the covariance structure of W will be assumed to be given by W i , W j t = ρt1 {i =j} + t1 {i=j} . We suppose that ρ ∈ (-1 d-1 , 1), which ensures that the matrix C = (ρ1 {i =j} + 1 {i=j} ) 1≤i,j≤d is positive definite. Let L denote the lower triangular matrix involved in the Cholesky decomposition C = LL * . To simulate W on the time-grid 0 < t 1 < t 2 < . . . < t N , we need d × N independent standard normal variables and set

       W t 1 W t 2 . . . W t N-1 W t N        =         √ t 1 L 0 0 . . . 0 √ t 1 L √ t 2 -t 1 L 0 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . √ t N -1 -t N -2 L 0 √ t 1 L √ t 2 -t 1 L . . . √ t N -1 -t N -2 L √ t N -t N -1 L         G,
where G is a normal random vector in R d×N . The vector σ = (σ 1 , . . . , σ d ) is the vector of volatilities, δ = (δ 1 , . . . , δ d ) is the vector of instantaneous dividend rates and r > 0 is the instantaneous interest rate. We will denote the maturity time by T . Since we know how to simulate the law of (S t , S T ) exactly for t < T , there is no use to discretize equation (6.1) using the Euler scheme. In this Section N = 2.

European options

We want to study the numerical accuracy of our algorithm and to do that we first consider the case of European basket options for which we can compute benchmark price by using very efficient Monte-Carlo methods, see [START_REF] Jourdain | Robust adaptive importance sampling for normal random vectors[END_REF] for instance, while it is no more the case for American options.

In this paragraph, the parameter ω appearing in (5.5) is 0.

European put basket option. Consider the following put basket option with maturity T K -

1 d d i=1 S i T + (6.2)
Figure 2 presents the influence of the parameters M , n and q. The results obtained for n = 1000, M = 50, 000 and q = 1 (curve (+)) are very close to the true price and moreover we can see that the algorithm stabilizes after very few iterations (less than 10). • Influence of M : curves (+) and (×) show that taking M = 5000 is not enough to get the stabilization of the algorithm.

+ + + + + + + + + + + + + + + + + + + + × × × × × × × × × × × × × × × × × × × × * * * * * * * * * * * * * * * * * * * * ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ n=1000 M=50000 q=1 + n=1000 M=5000 q=1 × n=1000 M=50000
• Joined influence of n and M : curves (×) and ( ) show the importance of well balancing the number of discretization points n with the number of Monte-Carlo simulations M . The sharp behaviour of the curve ( ) may look surprising at first, since we are tempted to think that a larger numbers of points n will increase the accuracy. However, increasing the number of points but keeping the number of Monte-Carlo simulations constant creates an over fitting phenomenon because the Monte-Carlo errors arising at each point are too large and independent and it leads the approximation astray.

• Influence of q: we can see on curves (+) and ( * ) that decreasing the hyperbolic index q can lead a highly biased although smooth convergence. This highlights the impact of the choice of q on the solution computed by the algorithm.

To conclude, we notice that the larger the number of Monte-Carlo simulations is, the smoother the convergence is, but when the polynomial family considered is too sparse it can lead to a biased convergence.

European call basket option. Consider the following put basket option with maturity T Figure 3 illustrates the impact of the sparsity of the polynomial basis considered on the convergence of the algorithm. The smoothest convergence is achieved by the curve (+), ie when M = 30, 000, η = 5 and q = 0.6. The algorithm stabilizes very close to the true price and after very few iterations.

1 d d i=1 S i T -K + (6.3)
• Influence of η : for a fixed value of q, the sparsity increases when η decreases, so the basis with η = 4, q = 0.6 is more sparse than the one with η = 5, q = 0.6. We compare curves (+) (η = 5) and (×) (η = 4) for fixed values of q (= 0.6) and M (= 30, 000). We can see that for η = 4 (curve (×)) the algorithm stabilizes after 7 iterations, whereas for η = 6 (curve (+)) less iterations are needed to converge.

• Influence of M : for fixed values of η (= 5) and q (= 0.6), we compare curves (+) (M = 30000) and ( * ) (M = 5000). Using a large number of simulations is not enough to get a good convergence, as it is shown by curve ( * ).

Actually, when the polynomial basis becomes too sparse, the approximation of the solution computed at each step of the algorithm incorporates a significant amount a noise which has a similar effect to reducing the number of Monte-Carlo simulations. This is precisely what we observe on Figure 3: the curves (×) and ( * ) have a very similar behaviour although one of them uses a much larger number of simulations.

American options

In this paragraph, the penalisation parameter ω appearing in (5.5) is 1.

Pricing American put basket options. We have tested our algorithm on the pricing of a multidimensional American put option with payoff given by Equation (6.2) Figure 4 presents the influence of the parameters M and q.

1
+ + + + + + + + + × × × × × × × × × * * * * * * * * * ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ Δ Δ Δ Δ Δ Δ Δ Δ Δ M=5000
• Influence of M : when zooming on Figure 4, one can indeed see that the curves using 30000 Monte-Carlo simulations are a little smoother than the others but these extra simulations do not improve the convergence as much as in Figures 2 and3 (compare curves (+) and ( ), Figure 4). The main explanation of this fact is that put options have in general less variance than call options and in Figure 2 a maturity of T = 3 was used which leads to a larger variance than with T = 1.

• Influence of q : once again, we can observe that increasing the sparsity of the polynomial basis (ie, decreasing q) can lead to a biased convergence. When q = 0.6, we get a biased result (see curves ( * ) and (△)), even for M large (curve (△), M = 30000).

Then, it is advisable for American put options to use almost full polynomial basis with fewer Monte-Carlo simulations in order to master the computational cost rather than doing the contrary.

Hedging American put basket options. Now, let us present the convergence of the approximation of the delta at time 0. Table 1 presents the values of the delta of an American put basket option when the iterations increase. We see that the convergence is very fast (we only need 3 iterations to get a stabilized value). The parameters of the algorithm are the following ones: n = 1000, M = 5000, q = 1, η = 3 and ω = 1.

Iteration

∆ 1 ∆ 2 ∆ 3 ∆ 4 ∆ 5 1 -0.

Dupire's framework

We consider a d-dimensional local volatility model in which the dynamics under the riskneutral measure of each asset is supposed to be given by

dS i t = S i t ((r -δ i )dt + σ(t, S i t )dW i t ) S 0 = (S 1 0 , . . . , S d 0 )
where W = (W 1 , . . . , W d ) is defined and generated as in the Black-Scholes framework. The local volatility function σ we have chosen is of the form σ(t, x) = 0.6(1.2e -0.1t e -0.001(xe rt -s) 2 )e -0.05 √ t , (6.4) with s > 0. Since there exists a duality between the variables (t, x) and (T, K) in Dupire's framework, one should choose s equal to the spot price of the underlying asset. Then, the bottom of the smile is located at the forward money. The parameters of the algorithm in this paragraph are the following : n = 1000, M = 30000, N = 10, q = 1, η = 3.

Pricing and Hedging European put basket options. We consider the put basket option with payoff given by (6.2). The benchmark price and delta are computed using the algorithm proposed by [START_REF] Jourdain | Robust adaptive importance sampling for normal random vectors[END_REF], which is based on Monte-Carlo methods.

Concerning the delta, we get at the last iteration the following vector ∆ = (-0.062403 -0.061271 -0.062437 -0.069120 -0.064743). The benchmark delta is -0.0625.

Pricing and Hedging American put basket options. We still consider the put payoff given by (6.2). In the case of American options in local volatility models, there is no benchmark. However, Figure 6 shows that the algorithm converges after few iterations. We get a price around 6.30. At the last iteration, we get ∆ = (-0.102159 -0.102893 -0.103237 -0.110546 -0.106442).

Speed up.

Remark 3. Because in high dimensions, the sequential algorithm can run several hours before giving a price, we could not afford to run the sequential algorithm on the cluster to have a benchmark value for the reference CPU time used in the speed up measurements. Instead we have computed speed ups as the ratio speed up = CPU time for 8 processors / 8 CPU time for n processors × n (6.5)

This explains why we may get in the tables below some speed ups slightly greater than 1. Our goal in this paper was to design a scalable algorithm for high dimensional problems, so it is not surprising that the algorithm does not behave so well in relatively small dimension as highlighted in Table 2. In dimension 3, the speed ups are linear up to 28 processors but then they dramatically decrease toward zero: this can be explained by the small CPU load of the computation of the correction term at a given point (t i , x i ). The cost of each iteration of the loop line 4 of Algorithm 1 is proportional to M pd 3 and when d is small so is pthe number of polynomials of total degree less or equal than η. For instance, for d = 3 and η = 3, we have p = 20, which gives a small complexity for each iteration over i. Hence, when the number of processors used increases, the amount of work to be done by its processor between two synchronisation points decreases to such a point that most of the CPU time is used for transmitting data or waiting. This explains why the speed ups decrease so much. Actually, we were expecting such results as the parallel implementation of the algorithm has been designed for high dimensional problems in which the amount of work to be done by each processor cannot decrease so much unless several dozens of thousands of processors are used. This phenomena can be observed in Table 3 which shows very impressive speed ups: in dimension 6, even with 256 processors the speed ups are still linear which highlights the scalability of our implementation. Even though computational times may look a little higher than with other algorithms, one should keep in mind that our algorithm not only computes prices but also hedges, therefore the efficiency of the algorithm remains quite impressive. 2: Speed ups for the American put option with d = 3, r = 0.02, T = 1, σ = 0.2, ρ = 0, S 0 = 100, K = 95, M = 1000, N = 10, K = 10, n = 2000, r = 3, q = 1, ω = 1. See Equation (6.5) for the definition of the "Speed up" column.

Conclusion

In this work, we have presented a parallel algorithm for solving BSDE and applied it to the pricing and hedging of American option which remains a computationally demanding problem for which very few scalable implementations have been studied. Our parallel algorithm shows an impressive scalability in high dimensions. To improve the efficiency of the algorithm, we could try to refactor the extrapolation step to make it more accurate and less sensitive to the curse of dimensionality.

Nb proc. 3: Speed ups for the American put option with d = 6, r = 0.02, T = 1, σ = 0.2, ρ = 0, S 0 = 100, K = 95, M = 5000, N = 10, K = 10, n = 2000, r = 3, q = 1, ω = 1. See Equation (6.5) for the definition of the "Speed up" column.

Figure 1 :

 1 Figure 1: Distribution of the CPU time needed for the creation of one Mersenne Twister generator

Figure 2 :

 2 Figure 2: Convergence of the algorithm for a European put basket option with d = 5, ρ = 0.1, T = 3,S 0 = 100, σ = 0.2, δ = 0, r = 0.05 K = 100, η = 3, ω = 0.The benchmark price computed with a high precision Monte-Carlo method yields 2.0353 with a confidence interval of (2.0323, 2.0383).

Figure 3 :

 3 Figure 3: Convergence of the price of a European call basket option with d = 10, ρ = 0.2, T = 1,S 0 = 100, σ = 0.2, δ = 0, r = 0.05, K = 100, n = 1000 and ω = 0. The benchmark price computed with a high precision Monte-Carlo method yields 7.0207 with a confidence interval of (7.0125, 7.0288).

Figure 4 :

 4 Figure 4: Convergence of the price of an American put basket option with d = 5, ρ = 0.2, T = 1, S 0 = 100, σ = 0.25, δ = 0.1, K = 100, r = 0.05 and ω = 1.

  Convergence of the delta for an American put basket option with d = 5, ρ = 0.2, T = 1, S 0 = 100, σ = 0.25, δ = 0.1, K = 100, r = 0.05.

Figure 5 :

 5 Figure 5: Convergence of the algorithm for a European put basket option with d = 5, ρ = 0, T = 1,S 0 = 100, δ = 0, K = 100, r = 0.05, η = 3, ω = 0. The benchmark price computed with a high precision Monte-Carlo method yields 1.745899 with a confidence interval of (1.737899, 1.753899).

Figure 6 :

 6 Figure 6: Convergence of the algorithm for an American put basket option with d = 5, ρ = 0, T = 1, S 0 = 100, δ = 0.1, K = 100, r = 0.05 and ω = 1.

  the distribution.

	3.2
	2.8
	2.4
	2.0
	1.6
	1.2
	0.8
	0.4
	0
	0.4	0.8	1.2	1.6	2.0	2.4


Ines Ben Jemaa
email: ines.benjemaa@inria.fr

Manabu Tsukada
email: manabu.tsukada@inria.fr

Hamid Menouar
email: hamid.menouar@hitachi-eu.com

Thierry Ernst
email: thierry.ernst@inria.fr

Validation and evaluation of NEMO in VANET using geographic routing

The combination of geographic-based routing protocols (GeoNetworking) and IPv6 NEtwork MObility (NEMO) into a single communication architecture (IPv6 GeoNetworking) is key in Vehicular Ad-hoc Networks (VANET). While NEMO manages Internet access and session continuity between the vehicle and the Internet, geographically based data forwarding allows an efficient dissemination of the information between vehicles and the infrastructure. In this paper, we refer to the basic scenarios that led to the design of the IPv6 GeoNetworking architecture in the context of the GeoNet project. A prototype implementation of the modules that couple these two technologies is described, in particular the adaptation of IPv6 and C2CNet, a layer that ensures the geographic capabilities. Results of a light experimental performance evaluation are reported.

I. INTRODUCTION

Intelligent Transportation Systems (ITS) gained a lot of attention in the past few years. New access technologies (i.e. 802.11p) saw the light in order to meet Vehicular Adhoc NETworks (VANET) needs such as highly dynamic topologies. Moreover, promising road safety, traffic efficiency and infotainment services have to rely on a network layer that could cover most of the communication scenarios and types which is a challenging target in VANETs. Vehicles are expected to transmit information to the surrounding vehicles as well as to the infrastructure and peers reachable in the Internet. This requires, on one hand, specific routing mechanisms to disseminate efficiently critical information in their direct surrounding and on the other hand the capability to maintain permanent Internet connectivity. Regarding the first point, location capabilities such as the Global Position System (GPS) have to be deployed in order to distribute data based on geographic routing decisions. For continuous Internet reachability, which is highly required for value-added services and usual Internet applications, IPv6 mobility support functions such as NEMO Basic Support (RFC 3963) are required. In this scope, coupling VANETs geographicbased routing (GeoNetworking) and NEMO into a single communication architecture supporting both scenarios is key. A mobile router embedded in the vehicle could then be used to maintain Internet connectivity for all in-vehicle nodes (navigation system, PDAs, etc.) while a VANET geographic routing protocol would allow communication with neighbor vehicles.

A few recent studies have dealt with the combination of VANET routing protocol together with NEMO. A mobile gateway allows the in-vehicle nodes reachability using a permanent prefix assignment. [START_REF] Wakikawa | Design of vehicle network: mobile gateway for MANET and NEMO converged communication[END_REF] explains how a mobile gateway could maintain multiple paths to the Internet and to the destination using a mobile ad-hoc (MANET) routing protocol and lists the advantages of using MANET and NEMO converged communication that allows fault tolerance and scalability. [START_REF] Tsukada | Simultaneous usage of NEMO and MANET for vehicular communication[END_REF] presents experimental results of the cited concept by testing the simultaneous use of OLSR as a MANET routing protocol and NEMO in a vehicular communication through several access technologies and using specific routing policies. [START_REF] Baldessari | A MANET-centric solution for the application of NEMO in VANET using geographic routing[END_REF] analyses the requirements of an efficient communication in which NEMO complies with VANETs and points out the advantages of a MANET-centric approach that includes a reactive mechanism that manages the MANETs routes and can switch back to a NEMO route. In addition to the research studies, worldwide, many organizations and consortia work on the design of a communication architectures that consider NEMO as a protocol running on top of a VANET routing layer, for instance the ITS station architecture commonly specified by ETSI TC ITS [START_REF] Tc | Intelligent Transport Systems (ITS)[END_REF] and ISO TC204 [START_REF]Intelligent Transport Systems -Communication for Land Mobiles (CALM) -IPv6 Networking; ISO/FDIS 21210[END_REF] and partly implemented in the CVIS project [START_REF] Sorasen | CVIS: CALM Proof of Concept Preliminary Results[END_REF].

Among the above cited work and studies, few real experiments that couple NEMO and geographic-based routing decisions have been performed. This was actually achieved in the context of the GeoNet european project which specified and implemented a communication architecture combining IPv6 and GeoNetworking (IPv6 GeoNetworking) [START_REF][END_REF].

This paper thus presents work performed in the GeoNet project. Section 2 presents the reference architecture contributed by GeoNet. Section 3 explains the design and the implementation of the system. Section 4 reports the results of the experimental performance evaluation tests realized in our testbed. The last section summarizes and concludes this paper.

II. REFERENCE SYSTEM AND ARCHITECTURE

Safety and non-safety applications should rely on a single communication architecture that is expected to provide both Particularly, the GeoNet european project 1 focused on such communication scenarios including the infrastructure. It aimed at combining Internet Protocol Version 6 (IPv6) features with geographic routing capabilities for VANETs into a single communication stack, referred to as IPv6 GeoNetworking and detailed in [START_REF][END_REF]. [START_REF] Choi | IPv6 support for VANET with geographical routing[END_REF] provides further information on the motivation for combining IPv6 with GeoNetworking.

As a part of the IPv6 GeoNetworking architecture, the C2CNet layer capabilities, designed in the frame of the Carto-Car Communication Consortium 2 (C2C-CC) are combined with IPv6 in order to enable Vehicle-to-Vehicle, Vehicle-to-Infrastructure and Internet-based communication. Three types of nodes are considered in the architecture: GeoNet OBU embedded in the vehicle, GeoNet RSU (Access Routers) deployed on the roadside infrastructure and other IPv6 nodes running GeoAware applications. Only GeoNet OBU and GeoNet RSU comprise the C2CNet layer capabilities; they are forming a GeoNet domain, that may comprise nodes implementing the C2CNet layer capabilities but not IPv6 (C2CNet nodes). The 1 EU FP7 GeoNet european project: http://www.geonet-project.eu 2 C2C-CC: http://www.car-to-car.org GeoNet OBU (On Board Unit) is an IPv6 Mobile Router (MR) connecting the in-vehicle network to other vehicles, the roadside infrastructure or the Internet. The GeoNet RSU (Road Side Unit) is an IPv6 Access Router (AR) connected to the roadside infrastructure network and providing Internet access to the OBUs in its communication range.

In this paper, we consider two basic scenarios. The first one is when classic Internet services (e.g. infotainment, video-ondemand or weather status information) are used. In this case, the vehicle has to communicate with some stationary nodes in the Internet. As shown in Fig. 1, the packets are sent from the IPv6 node attached to the GeoNet OBU. The OBU implements NEMO Basic Support (NEMO BS) to maintain its reachability in the Internet when moving from one network to another. By means of the Neighbor Discovery Protocol [START_REF] Narten | Neighbor Discovery for IP version 6 (IPv6)[END_REF], the GeoNet OBU can select the access router to which it has to deliver each packet destined to the Internet. The packets are then forwarded on the C2CNet link. Intermediate C2CNet nodes relay the packet at the C2CNet layer until it reaches the GeoNet RSU (next IP hop from the GeoNet OBU) which forwards them to the Home Agent (HA). The HA routes them to the destination node. In the second scenario, a vehicle is expected to send alerting messages to other surrounding vehicles whenever it detects a road traffic hazard. One of the common cases to consider, is when the vehicle has no reachability with other vehicles. In this case the safety information has to be sent through the infrastructure to a control center or directly to the other vehicles belonging to the same geographic area (i.e. the same highway). The following sections give more details about the C2CNet layer functions, The NEMO BS and the Neighbor Discovery Protocol.

A. The C2CNet layer concept

C2CNet is a communication layer that enables geographical addressing and routing. C2CNet includes position-based routing mechanisms adapted to vehicular communications. The Greedy Perimeter Stateless Routing (GPSR) [START_REF] Karp | GPSR: greedy perimeter stateless routing for wireless networks[END_REF] algorithm has been adopted. It benefits from the reactive approaches of GeoRouting in wireless networks where the route determination is initiated on demand. GPSR is based on the greedy forwarding decisions using only information about immediate neighbors in the network topology. C2CNet defines a new network header (see Fig. 2) which carries the C2C identifier (C2CNet ID) of the source and the destination and their geographic locations. Each node in the vehicular ad-hoc domain is addressed by a unique C2CNet ID that it exchanged between one vehicle and its neighbors. Thus, the routing decision is based on geographic location of communication peers, source, destination and intermediary nodes.

C. The Neighbor Discovery Protocol in GeoNet

Each C2CNet egress interface of a GeoNet OBU (MR) or a GeoNet RSU (AR) must be configured with two different IPv6 addresses: a link-local address and a global address. All GeoNet nodes attached to the same IPv6 C2CNet link should be reachable using both addresses. The global address must be used when trying to reach other nodes not directly attached to the C2CNet link (i.e. nodes on the global Internet as well as nodes attached to the GeoNet OBUs and GeoNet RSUs). The mechanism used to configure the IPv6 C2CNet egress interfaces of GeoNet OBUs in an automatic way is based on the IPv6 Stateless Address Autoconfiguration protocol specified in RFC 4862 [START_REF] Narten | Neighbor Discovery for IP version 6 (IPv6)[END_REF]. This protocol basically enables a host to generate its own addresses using a combination of locally available information (interface identifier part of the address) and information advertised by routers (prefixes that identify the subnets associated with a link). The concept used in GeoNet is the same: GeoNet RSUs advertise the prefix information by sending Router Advertisements, and the GeoNet OBUs use that prefix information and their C2CNet IDs, to generate a valid global IPv6 address assigned to the C2CNet egress interfaces. A link-local address is also generated on the C2CNet egress interface, using the same C2CNet interface identifier and the link-local IPv6 prefix (FE80://64). A prefix length of 64 bits is used and thus the length of the C2CNet ID is 64 bits. IPv6 multicast packets produced by Neighbor Discovery, i.e. the all-nodes multicast address (FF02::1) and the all-routers multicast address (FF02::2) are mapped to a geographic area of delivery (GeoDestination) when transmitted to the C2CNet layer. Since Router Advertisements are typical IPv6 multicast packets, they are encapsulated at the C2CNet layer and forwarded as GeoBroadcast within a well delimited geographical area. The GeoDestination of the Router Advertisement is set up at the C2CNet level.

III. DESIGN AND IMPLEMENTATION

To test our contribution, we have implemented and integrated our work to Hitachi's C2CNet layer which was implemented within the GeoNet project. This C2CNet layer implementation is also used in other research projects, such as PRE-DRIVE C2X 3 . The specifications that have been used as basis in ETSI TC ITS standards are also inline with the V2X standard drafts. In GeoNet, modules and Service Access Points (SAPs) between layers have been defined. In this paper, we focus mainly on the C2C-IPv6 SAP between the IPv6 and C2CNet layers. The whole system is implemented in linux based on the 2.6.29 kernel version. The main modules involved are:

• The C2CNet module which implements the geographical routing functionalities as explained in the above section; • The PositionSensor module which provides geographical coordinates to the On-Board-Unit through GPS;

• The Lower Layer module whose role is to adapt the C2CNet to the MAC layer;

• The Service Access Point between the IPv6 layer and the C2CNet layer that performs the transmission and the encapsulation of the IPv6 packets through the C2CNet link.

As C2CNet is implemented in the user space, a virtual interface is used to deliver packets from IP layer implemented in kernel space to the C2CNet module. We used the NEPL [START_REF]Nautilus6[END_REF] implementation of NEMO which is installed in OBUs as well as in the Home Agent. The radvd4 software is installed in both the RSU to enable Router Advertisement in the GeoNet domain and in the OBU to advertise the in-vehicle prefix in the in-vehicle network. We used the madwifi driver for our Atheros cards in the 802.11g standard. Applications providing HMI to the users in the in-vehicle network are running in the nodes attached to the OBU. Once the virtual interface is up, it gets a link local address built from the FE80::/64 prefix and the C2CNet ID of the OBU. A patch is developed for NEMO to allow OBUs to configure their CoA using the C2CNet ID. This implementation is necessary to enable IPv6 address autoconfiguration using the C2CNet ID instead of the MAC address. The Binding Update refresh time is set up to 15 seconds. When RSU sends a Router Advertisement, it uses the multicast address ff02::1 (the all-node link local address). The RA packets are then GeoBroadcasted in an radius area of 500m from the RSU. OBUs receive the RA and update their routing tables setting the link local address of the RSU as a default gateway address. They also autoconfigure their CoA using the prefix advertised in the RA packet. The Correspondent Node traffic is encapsulated first into the NEMO tunnel and then into the C2CNet tunnel. Fig. 3 illustrates the internal functioning of the GeoNet implementation. When an Application Unit (AU) sends a packet to the OBU, the NEMO daemon intercepts the packet, encapsulates it into the NEMO tunnel and sends it to the C2CNet virtual interface. The C2CNet layer, listening on the tun0 interface, gets the packet, executes the routing algorithm based on the previous knowledges of its neighbors, decides where to send the packet and then forwards it locally to the lower layer. The lower layer sends finally the packet on the wireless link. A new ethernet type for the C2CNet header is defined to allow its encapsulation into the MAC header.

IV. TESTS AND EVALUATION

A. Evaluation environment

IPv6 GeoNetworking is evaluated in the indoor testbed. The indoor test environment is designed to evaluate the pure performance of IPv6 GeoNetworking avoiding interferences due to unexpected radio perturbations and difficulties to trace the movements of the GeoNet OBUs. The scenario is shown in Fig. 4. The GPS data is not obtained from actual GPS device but is statically recorded in a configuration file. The advantage of this method is that the same test scenario can be repeated several times with various parameters. The traffic is generated by the iperf tool. The two communication end-points are one AU attached to the GeoNet OBU that sends the traffic to be evaluated and a Correspondent Node, considered to be located in the Internet. In these tests, packets transit via the Home Agent. We evaluate two types of traffic: i) UDP traffic which is a unidirectional transmission flow from the source to the destination end-nodes where the considered metrics are the packet loss rate and the throughput, and ii) ICMPv6 traffic which is a bi-directional communication flow between the two end-nodes where the considered metrics are RTT and the packet loss rate.

B. Latency evaluation

To evaluate the latency, we measured the Round Trip Time between the two end-points. The AU sends ICPMv6 Request every 0.1 second. The ICMPv6 packet is increased by 20 bytes. The packet size is varying from 20 bytes to 1500 bytes. From the obtained results, we extract the maximum, the minimum and the average RTT as well as the packet loss for each packet size. A previous evaluation of IPv6 over C2CNet is already presented in [START_REF] Tsukada | Experimental Evaluation for IPv6 over VANET Geographic routing[END_REF] which doesn't include the use of NEMO. As depicted in Fig. 5, we evaluate the average RTT between the Correspondent node and the GeoNet Mobile Router and the packet loss. The maximum RTT is around 110 ms which corresponds to the maximum packet loss (45 percent) for 420 bytes of packet size. As we can see in Fig. 5, packets with size exceeding 1300 bytes cannot be delivered by C2CNet due to the MTU of the packet. At the time of writing this paper, the packet fragmentation operation was not yet implemented in the C2CNet layer.

C. Packet delivery ratio and bandwidth evaluation

In this tests, we evaluate the packet loss ratio in a UDP communication. The packet delivery ratio is the percentage of packets arriving at the receiver divided by the packets sent by the sender. The UDP packets are generated in the AU attached to the OBU, sent through the C2CNet link to the HA and finally to the Correspondent Node. The sender sends UDP packets to the receiver with fixed rate. The UDP client geographical routing. and server save the log file traces. After the tests, the log files of both the client and the server are parsed through pointers (the port number) and the packet loss results are plotted. In these tests, the bandwidth is varying from 1 to 6 Mb/s. For each bandwidth value, the read-write buffer is increased from 20 bytes to 1900 bytes. The throughput is shown on the receiver side. As illustrated in Fig. 6, when the packet has a small size, the packet delivery ratio is weak. The best values are obtained when the packet size is between 800 and 1300 bytes for the lowest sending rate; when the bandwidth is 1 M and the packet size is 1300 bytes, the packet delivery ratio is almost 100 percent. The maximum throughput is around 2500 Kbits/second. It reaches its maximum when the packet size is 1300 bytes packet. It corresponds to a 5M sending rate.

D. Result interpretation

Sub-optimal routing is caused by the packets being forced to pass via the HA. This leads to performance degradation due to increased delay and is undesirable for some applications. Packet Encapsulation of additional 40 bytes header increases packets overhead and may result into packet fragmentation. This turns the results into an increased processing delay for every packets being encapsulated and decapsulated in both the GeoNet OBU and the HA. Bottlenecks in the HA are a severe issue because significant traffic to and from MNNs is aggregated in the HA when it supports several GeoNet OBUs acting as gateways for several MNNs. This may cause congestion at the HA that would lead to additional packet delays, or even packet losses. This issue is subject to bringing further enhancements to this specification (Route optimization solutions) although it is not peculiar to IPv6 GeoNetworking.

V. CONCLUSION AND FUTURE WORK This paper shows how network mobility support capabilities using NEMO Basic Support can be combined with the GeoNetworking capabilities offered by the C2CNET sub-layer. We performed validation tests on our in-door testbed. The results show that a route optimization solution is highly required to resolve the performance issues and bottleneck in the Home Agent. As a future work we plan to conduct field tests with real vehicles and realistic mobility scenarios. Regarding the C2CNet layer implementation, further enhancements are required to avoid the packet processing performance issues.

Fig. 1 .

 1 Fig. 1. Reference scenario and architecture

Fig. 2 .

 2 Fig. 2. GeoNet packet encapsulation

Fig. 3 .

 3 Fig. 3. Implementation of IPv6 over C2CNet

Fig. 4 .

 4 Fig. 4. Network configuration of the indoor test

Fig. 5 .

 5 Fig. 5. RTT between AU and CN

Fig. 7

 7 Fig.7presents the measured throughput in the network.

Fig. 6 .

 6 Fig. 6. Packet delivery ratio between AU and CN

EU FP7 Project PRE-DRIVE C2X: http://www.pre-drive-c2x.eu

Linux Router Advertisement Daemon: http://www.litech.org/radvd/

ACKNOWLEDGMENT

The work reported in this paper has been financed by the European Commission under the GeoNet project (grant agreement 216269). The authors would also like to thank all the members of the GeoNet project and more specifically NEC Europe and IMDEA who provide us support for developing the code.