Antoine Beugnard
email: antoine.beugnard@telecom-bretagne.eu

Jean-Marc Jézéquel
email: jezequel@irisa.frplouzeau@irisa.fr

Noël Plouzeau

Contract Aware Components, 10 years after

The notion of contract aware components has been published roughly ten years ago and is now becoming mainstream in several fields where the usage of software components is seen as critical. The goal of this paper is to survey domains such as Embedded Systems or Service Oriented Architecture where the notion of contract aware components has been influential. For each of these domains we briefly describe what has been done with this idea and we discuss the remaining challenges.

Introduction

In a special issue of IEEE Computer published in July 1999, we published an article entitled Making Components Contract Aware. The goal of this article was to discuss the extension of Meyer's Design-By-Contract idea to the world of software components. Our finding were that component contracts would have to deal with specific concerns that we classified into four categories, from syntactic to semantic to synchronization, up to Quality of Service.

There is empirical evidence that this article has got some lasting impact on the field of Component Based Software Engineering. For instance, on June 2010 it has been cited 545 times by other papers (according to Google Scholar). Several of these papers themselves got quite a sizable number of citations. Impacted domains include1 :

Components: Software reuse strategies and component markets [START_REF] Ravichandran | Software reuse strategies and component markets[END_REF] (81 citations), Specification, implementation, and deployment of components [START_REF] Crnkovic | Specification, implementation, and deployment of components[END_REF] (80 citations), Interface Compatibility Checking for Software Modules [START_REF] Chakrabarti | Computer Aided Verification[END_REF] (78 citations)

Adaptation: Composing Adaptive Software [START_REF] Mckinley | Composing Adaptive Software[END_REF] (281 citations), A Taxonomy of Compositional Adaptation [START_REF] Philip | A Taxonomy of Compositional Adaptation[END_REF] (53 citations)

Internet (beyond components): Information agent technology for the Internet: A survey [START_REF] Klusch | Information agent technology for the Internet: A survey[END_REF] (199 citations), WSOL -Web Service Offerings Language [START_REF] Tosic | Web Services, E-Business, and the Semantic Web[END_REF] (106 citations)

Real-time: Monitoring, Testing and Debugging of Distributed Real-Time Systems [START_REF] Thane | Monitoring, Testing and Debugging of Distributed Real-Time Systems[END_REF] (90 citations)

High Performance Computing: A Component Architecture for High-Performance Scientific Computing [START_REF] Bernholdt | A Component Architecture for High-Performance Scientific Computing[END_REF] (86 citations)

This has prompted us to try to reflect upon the evolution of the contract-aware components domain roughly ten years after the publication of our article in IEEE Computer (we will refer to this article as the MCCA paper in the following sections). The goal of the present paper is thus to highlight domains where the notion of contract aware components has been influential, to briefly describe what has been done with this idea and to discuss the remaining challenges. The rest of this paper is organized as follows.

In Section 2, we briefly recall the main ideas discussed in our original paper. In Section 3 we give an overview of the state of the art with respect to contract-aware components in domains such as Embedded Systems or Service Oriented Architecture. We conclude the paper with a discussion on the possible evolutions of the contract concept.

Original paper summary

Ten years ago, we proposed in [START_REF] Beugnard | Making components contract aware[END_REF] to apply B. Meyer's Design by Contract [START_REF] Meyer | Applying "Design by Contract[END_REF] principles to components with an attempt at generalization. We introduced a classification of contracts in four levels, and beyond a simple verification use of contract, the possibility to manage and negotiate contracts at runtime.

The classification, like all classifications, expresses a point of view that we discuss a little further (section 4). It can be summarized as:

Syntactic (or basic) The goal is to make the system work. It is generally specified with Interface Definition Languages (IDLs), as well as typed object-based or object-oriented languages. It ensures the components can be assembled.

Behavioral

The goal is to specify each operation. It is generally specified with a couple of assertions: a precondition and a postcondition. It ensures the operations offered and required are not only syntactically compatible but also semantically.

Synchronization

The goal is to specify the coordination of operations. It can be specified with an automaton labelled with operations. It ensures the operations are used in the proper order.

Quality of Service

The goal is to quantify a few features associated to operations. Performance, availability and quality of result can be specified and negotiated at that level.

This classification clearly helps to structure the specification and to understand the coverage of requirements a component has.

The structural decomposition that we proposed also introduces a temporal decomposition with a contract management life-cycle. The main idea is to reify contracts (defining contracts as objects). The life-cycle can be summarized as follows:

Define Describe the component features with all required contracts. A component is considered in isolation.

Subscribe Select from all contracts those that are useful in the context of component use, and configure them. Components are assembled (used) with an intent.

Check Evaluate contracts and react accordingly. The moment the checking is carried out depends on the level of contract: usually levels 1 and 2 can be statically checked while levels 3 and 4 require some runtime monitoring. In some cases, level 3 can be statically checked. In the case of a contract violation we proposed 4 kinds of reactions: ignore, reject, wait or negotiate.

Terminate Decide when to stop evaluating contracts.

The 4-level structural classification is oriented towards the nature of the constraints that have to be specified. As we will see later (section 4) space and time information are other classification dimensions that can be considered.

In the next section, we focus on some application domains in order to analyze how contracts were used.

Related work

Contracts have been applied in many fields. They help to design and structure both the product and the production process. In this section we focus on examples which range from the less flexible architectures (closed embedded systems) to the most flexible ones (service oriented architectures) through an intermediate flexible architecture (component based systems). This shows that "contract awareness" is not limited to components, but can efficiently be extended to other paradigms such as services.

Embedded Systems

In recent years, the notion of contract aware components has clearly percolated in the domain of Real-Time Embedded Systems. For instance, the Artist network of excellence supported a research cluster on Component-Based Software Development where the notion of contract was central. However, a key characteristic of component-based embedded systems is the wide heterogeneity of component models. This heterogeneity can be found in different execution models (synchronous, asynchronous, vs. timed), different communication models (synchronous vs. asynchronous), as well as different scheduling paradigms. Designing heterogeneous embedded systems from diverse types of components, and allowing the prediction and optimization of functional and non-functional properties of the designed systems is still an open challenge [START_REF] Basu | Modeling Heterogeneous Real-time Components in BIP[END_REF]. Thus there is a need to develop an innovative theory for building complex heterogeneous systems which better addresses composability and compositionality issues [START_REF] Sifakis | Component-Based Construction of Real-Time Systems in BIP[END_REF]. Such a comprehensive theory is still missing today, thereby making it difficult to understand how to build systems that combine, e.g., synchronously and asynchronously executing components and reason about non-functional properties, but this is the subject of very active research (e.g. as published in the Emsoft conferences).

It is also now acknowledged that a key issue in component-based embedded software development is about handling non-functional properties (including real-time and QoS properties). The notion of a rich component model [START_REF] Damm | Controlling Speculative Design Processes Using Rich Component Models[END_REF][START_REF] Benveniste | Multiple Viewpoint Contract-Based Specification and Design[END_REF] is thus gaining momentum to help engineers to model, specify, and predict timing, QoS, and resources properties of components and of systems composed from components. Still, typical support for handling QoS and resource usage is rather limited.

Specific non-functional (level 4) contracts are dealt with in ad hoc manners. For timing properties, different variants of timed automata have been used, as in, e.g., the Omega component model [START_REF] Damm | A discrete-time UML semantics for concurrency and communication in safety-critical applications[END_REF]. For properties relating to queuing and performance, models based on queueing networks, Markov chains, etc. have been used. These approaches offer a precise mechanism for specifying and analyzing QoS properties, but they suffer from scalability problems.

In Real-Time Embedded Component-Based Systems, level 3 contracts are also of prime importance to manage both expectations of the component about its environment and guarantees offered in return by the component to its environment. An example of level 3 contracts is the notion of interface automata [START_REF] Luca De Alfaro | Timed Interfaces[END_REF], viewed as enriched type systems (the so-called behavioral type systems), which capture the ordering aspects of software component interactions. Interface automata provide a semantically well-founded, built-in notion of refinement: a component refines another one if it imposes weaker constraints about the environment and offers stronger guarantees in return [START_REF] Chakrabarti | Computer Aided Verification[END_REF].

UML emerged in recent years as a modeling standard for software, including software for embedded systems for which specific UML profiles have been developed. A first attempt was made with the UML Profile for Schedulability, Performance, and Time (SPT) to model real-time concerns. However, the SPT profile suffered from several shortcomings, and was quickly superseded by the MARTE profile (Modeling and Analysis of Real-Time and Embedded systems), which better addresses issues such as compliance with the UML2.0 specification of not only real-time constraints but also other embedded In this domain, SPEEDS (SPEculative and Exploratory Design in Systems Engineering) is another example of a recently completed European project that developed contracts and multi-viewpoint as core concepts (see Figure 1).

SPEEDS has delivered a solution based on enriching functional decomposition with contracts, i.e. textual and formal descriptions of functional and non-functional aspects. The SPEEDS version of contracts give abstract descriptions of the exposed properties of a component, classified as assumptions (requirements of the component with respect to its environment) and promises (guaranteed delivered properties of the component provided its assumptions have been fulfilled). According to the documentation produced by the SPEEDS project2 , typical contracts specifying the behaviour of a component (used as promises) could be for instance:

• The output out is the sum of the inputs in1 and in2.

• Every request will be served.

• A request will be served within 10 ms.

Conversely, examples of behavioural contracts used as assumptions could be:

• Input data in1 will never show a negative value.

• The data rate of in1 will not exceed one message in 10 ms. SPEEDS has provided a complete framework for modeling, combining, analyzing and managing such multi-viewpoint, enriched models.

Component architectures

Contracts have gained first class status in recent years. The set of concepts used to describe and manage them is now very rich. More than ten years ago, the first approaches relied on a static model of contract. For instance, an operation definition may provide a non negotiable precondition (the caller cannot negotiate the precondition) and also a non negotiable postcondition. While this may be sufficient for level 2 contracts (especially if one considers level 2 contracts as a type mechanism) this scheme is too static to maximise reuse and adaptation to evolving needs. Hence a negotiation model is required. In turn this implies that contracts have complex features and that they are much more than a pair of precondition/postcondition predicates.

In this section we use the Confract approach [START_REF] Collet | Composite contract enforcement in hierarchical component systems[END_REF] as an example architecture to illustrate elaborate means of contract management (especially negotiation and monitoring). Confract is a contract management framework for the Fractal [START_REF] Bruneton | The fractal component model and its support in java[END_REF] component model. Fractal has a rich component life cycle, which in turn allows for rich contract lifecycle with high possibilities regarding dynamic contracts. As in most modern contract based systems, a Confract contract is a negotiated agreement between several parties whose responsibilities are clearly stated for each contract feature. Confract allows for several kind of contracts, depending on the location that receives them:

• interface contracts describe agreement between a service provider and a service client;

• external composition contracts express compound usage rules for a component, defining how a component can be used globally (such contracts encompasses several interfaces);

• internal composition contracts define how internal components are connected to external interfaces.

The various contract locations mentioned above show that the cooperation properties are more central than in the case of initial design by contract approach.

The set of specifications is at the heart of Confract: a contract is a "collaboration hub". Through a set of operations entities (component instances), the partners select the specification items that they are interested in (either as service provider of service consumer) with adequate values for the specification parameters. Once partners have reached an agreement on a specification term, this term is "closed". A closed term includes a contract manager (which is in charge of managing the violation of this term), contributors (which strive at ensuring the contract term by providing services) and beneficiaries (which are customers of the contract bound services). When all terms are closed, the contract can be instantiated and used.

In a typical architecture, an instantiated contract is a reference to monitor the activities of the partner and detect contract violations. This contract instance is at the center of the collaboration management. In service based systems and also in evolving component based systems configurations can be computed on the fly and therefore contract models must also take these evolutions into account. Techniques such as model driven engineering help to master reconfiguration and they play a key role in contract management. More precisely, the "model at runtime" technique [START_REF] Morin | Models at Runtime to Support Dynamic Adaptation[END_REF] allows for a precise reconfiguration control of a running system. This includes the analysis of contract validity on an evolving architecture.

Service Oriented Architecture

Contract management is central to Service Oriented Architecture design. With respect to component based architectures, service oriented ones are characterized by the highly dynamic configuration and reconfiguration capabilities that they require. As in component based architecture, SOA contracts are the key concept to setup cooperation between entities. But SOA contracts have a richer lifecycle: they are not set in stone at design time but are much more malleable than design time component contracts. Conditions that were favorable to the definition of a cooperation at a given time and place may change at any moment and this change lead to a deal break. In turn this implies a reaction from the contract partners, from a basic renegotiation of bounds to a contract cancellation. In the spirit of service oriented architectures, such contract volatility is not a bug but an essential feature.

This high level of dynamicity has a profound impact on the way to create, configure and supervise contract. This leads to a large set of concepts (see [START_REF] Keller | The WSLA framework: Specifying and monitoring service level agreements for web services[END_REF] for an example of such conceptual framework). Hereafter we will mention some points on configuration and monitoring of contracts in a SOA context.

The terminology of the SOA domain includes the notion of service level agreement (SLA). There is a connection between the contract notion (in the sense described in the previous section on components) and the SLA concept. Both concepts share the idea of mutual agreement on provided and consumed services. However, while component contracts are based on precise notions of quality, referring to low level items of execution (e.g. maximum time between service invocation and service termination events), SLA properties refer to higher level properties, abstracting many execution details and often using stochastic descriptions. SLA deals with the global performance of a system as a service provider, and SLA descriptions include various metrics concerns and dimensions, such as business metrics: financial properties, for instance a billing computation formula based on stochastic behavior of a compound metrics (e.g. mean execution time of a service during the last hour, etc).

Contract creation and configuration

Strictly speaking, contracts are created from offers, once offers have been read, chosen and configured. In the component design world, components are interconnected through ports and in many component metamodels, the relationship between a provided port and a required port is bidirectional and represented by a connector (explicitly or implicitly). A contract is built from a required properties part bound to the requiring/client port, matched with provided properties bound to the provider port.

Regarding negotiation aspects, SOA offers a wide range of possibilities for every contract level, in the sense of initial MCCA article:

1. level 1 negotiation is realized by the service discovery mechanism inherent to service based systems; the notions of type compatibility are managed there;

2. level 2 negotiation is also connected to service discovery: pre and postconditions are a form of type definition that extends data type definitions;

3. level 3 is realized by the definition of an orchestration; usually negotiation is limited because orchestrations are more static than the services they rely upon;

4. level 4 is realized by service level agreement mechanisms.

Monitoring services This dynamic approach requires that SOA contracts have to be supervised at all times, contrary to a more static view of contract composition analysis from the embedded or component systems design. In practice, the monitoring system is generated and configured by the contract management infrastructure. Monitoring can be performed by external parties in order to obtain a neutral point of view of the compliance of the interested parties to the contract. With respect to the classification of the MCCA paper, monitors can manage different contract levels: monitoring systems such as the one in [START_REF] Baresi | Applying 'Design by Contract[END_REF] rely on an infrastructure that can accept plugins for various monitoring tactics, with a rule-based engine to trigger reaction strategies. Reaction strategies may adapt monitors and monitoring rules, as well as strategies that retry, rebind or reorganize.

Regarding monitoring, service based architectures address the four MCCA levels with very different techniques:

1. level 1 contracts can be managed by WSDL alike service descriptions, with a type system that is loose enough to deal with heterogeneity and fast paced evolution; contract compliance is computed by loose type computation and checking;

2. level 2 contracts can be described by pre and postcondition and can be checked by weaving monitors, inserting adapters, etc [START_REF] Baresi | Applying 'Design by Contract[END_REF];

3. level 3 contracts can be addressed by orchestrations, which somewhat can ease contract monitoring because the control is centralized;

4. level 4 contracts can be addressed by SLA validation and monitoring Monitoring techniques are different for embedded architectures and service based systems:

• service based designs have limited possibilities of static analysis tools because of the highly dynamic nature of these systems. Checking can be done at run time only: at contract negotiation time, during the execution of service requests, at system reconfiguration time;

• the consequences of contract violation differ between traditional component based systems and service based ones. The goals of monitoring are therefore different: for embedded architectures, monitoring is a complement to more static analysis and validation techniques; for service based systems monitoring is a first class concept that is an integral part of the contract life cycle.

• monitors for embedded systems can be generated statically and injected in the main system code infrastructure to allow for efficient violation detections and fast error management [START_REF] Saudrais | From formal specifications to QoS monitors[END_REF]; on the contrary dynamic systems such as service based architectures are more difficult to generate in advance and violation management policies can also be defined dynamically.

Other work

The initial application field of the MCCA paper was components. A wide interpretation of the word component is allowed and leads to several applications of contracts. Moreover, as discussed in section 4, a contract can be attached to many kinds of component: from operation (tiny "component") to system (huge "component"). In this section we describe briefly other work using contracts. A recent interesting work has been realized by [START_REF] Waignier | Canevas de développement agile pour l'évolution fiable de systèmes logiciels à composants et orientés services[END_REF]. The author proposed an agile framework for the evolution of systems based on components and services. Figure 3 shows the whole development process. Four types of contracts (gray boxes) are used at different stages of the process with different tools. Some contracts are used with static checkers in order to prove architectural static properties while other contracts are transformed and integrated into the code to enable observations and runtime monitoring. This work shows that contracts can be a natural interface in order to mix and assembly components and services as in the Service Component Architecture initiative.

In 2000, the OMG introduces the MDA (Model Driven Architecture) initiative [START_REF] Soley | Model Driven Architecture[END_REF]. Beyond the UML or MOF dependency of this approach, this initiative promotes a systematic use of models (whatever the modeling language) and of model transformations (whatever the transformation language) that describes design decisions. In [START_REF] Kaboré | Contribution à l'automatisation d'un processus de construction d'abstractions de communication par transformations successives de modèles[END_REF], E. Kaboré used contracts to specify all transformations; constraints on input and output models are considered as pre and postconditions on transformations.

More recently, A. Koudri proposed Modal [START_REF] Koudri | MODAL: A SPEM Extension to Improve Co-design Process Models[END_REF] a process model to interpret parts of a software development process as "process components" which, hence, allows to attach contracts to them. Figure 2: A full process with contracts from [START_REF] Waignier | Canevas de développement agile pour l'évolution fiable de systèmes logiciels à composants et orientés services[END_REF] 4 Discussion on the classification As for any classification, our four levels contract classification expresses a point of view and can be discussed. The point of view of the MCCA paper is to consider the component as a whole and not to consider it as a set of parts such as interfaces, ports or other "components". A contract is then attached to a component, but this position is easily reconsidered with contracts attached to components parts. Basic and semantic contracts are naturally linked to operations, themselves grouped into interfaces or ports whatever the boundary of the component is composed of. Synchronization (or coordination) and quality of service contracts describe more global features and are easily attached to the component itself. In his thesis, G. Waignier [START_REF] Waignier | Canevas de développement agile pour l'évolution fiable de systèmes logiciels à composants et orientés services[END_REF] proposes a systematic approach using contracts at many structural levels: operation, port, whole component, assembly of components (this is made possible by the Confract component model [START_REF] Collet | Composite contract enforcement in hierarchical component systems[END_REF])

The levels that are the most subject to interpretation are the behavior and synchronization ones. This is probably because synchronization is a part of the component behavior and that languages used to specify operation semantics in preconditions and postconditions are rich enough to describe control states. In a strict interpretation of the behavioral contract, the operation specification is stateless and makes no assumption on the state of the component. For instance, a deposit in a BankAccount is balance += amount; whatever the BankAccount status (non-existing, open, closed, deleted, etc.) is. Unfortunately, languages such as OCL allows references to variables that can encode a state, introducing a synchronization-level issue in the behavioral level. In order to remove this moving borderline, the Accord project [START_REF] Beugnard | Assemblage de composants par contrats. Le modèle de composants ACCORD[END_REF] proposes a classification in only three levels: syntactic, semantic and pragmatic. But, this usual classification in language theory introduces a new kind of ambiguity; pragmatic aspects usually denote the way things are used. Therefore, is the synchronization contract that describes the way operations are used a pragmatic contract (from the operation point of view) or a semantic contract (from the component -as a unit -point of view)?

Finally, a completely different and complementary classification is based on the moment a contract is designed within the design process. All contracts are defined and refined during the whole development process. A more precise classification can define requirement contracts, specification contracts, design contracts, implementation contracts and runtime contracts. UML Components [START_REF] Cheesman | UML Components[END_REF] and G. Waignier [START_REF] Waignier | Canevas de développement agile pour l'évolution fiable de systèmes logiciels à composants et orientés services[END_REF] rely on such an approach.

As a conclusion for this discussion on contract classifications we can distinguish three dimensions :

1. the nature of the contract, as in the MCCA initial proposal 2. the location of the contract;

3. the process moment (time position) of the contract.

The nature of a contract depends on its attachment and its position in the development cycle of a component. All these contracts have compatibility constraints that need to be analyzed along the process. The classification of the nature of the contract is finally a convenient tool to decompose constraints, requirements or rules and is probably business dependent. In fact, the classification is just a tool that helps improve the trust we can have in the coverage (not to say the completeness) of a contract-based specification.

Conclusion

Beyond the primary idea of Design by Contracts [START_REF] Meyer | Applying "Design by Contract[END_REF] which was a pragmatic and systematic way of using pre/post conditions in a functional (or object-oriented) world, the idea to generalize the use of contracts to components leads to many interesting works.

Contracts are a natural way to decompose and structure specifications. It appears that a contract can be attached to almost any software entity. One of the main challenges is probably to deal with the many kinds of contracts and of contracts checkers. But even if contacts are widely used, they are not yet a first-class entity in UML or even standard profiles.

A fundamental question that arises is the role of contracts in the specification of software entities. Can we consider the set of contracts as "the" specification of a service or a component? Another important question is the life cycle of contracts. Some contracts can be directly used and statically checked while others need to be refined and give birth to monitors, checkers and negotiation artifacts. The refinement process of contracts depends on the level of contract and is tightly coupled to the software entity development process itself with the underlying question; what has to be exposed in the contract, what has to be hidden or, in other words, what is the level of grayness of the specification. 3Beyond the specification of a single entity, the contract concept is now central to many software architectures. Yet many complex questions need to be answered: which contract metamodels are relevant for a given system, how can designers conciliate different metamodels stemming from various parts that are combined to build a system, how can one integrate dynamic contract negotiation, etc.

So, ten years later, the concept of contract still appears attractive but has probably not reached a level of maturity that would allow a large use in industry. Rendezvous in ten more years . . .

Figure 1 :

 1 Figure 1: The SPEEDS vision: Design enriched by contracts

Number in parenthesis is the number of citations of this paper according to Google Scholar.

http://www.speeds.eu.com/

On this topic, the WCSI paper[START_REF] Ouederni | Tau Be or not Tau Be? A Perspective on Service Compatibility and Substituability[END_REF] provides an interesting view.


Marie Bernhart
email: marie-externe.bernhart@edf.fr

Peter Tankov
email: peter.tankov@polytechnique.org

Xavier Warin
email: xavier.warin@edf.fr

91g20

A finite dimensional approximation for pricing moving average options

Keywords: American options, indexed swing options, moving average, finite-dimensional approximation, Laguerre polynomial, least squares Monte Carlo

We propose a method for pricing American options whose pay-off depends on the moving average of the underlying asset price. The method uses a finite dimensional approximation of the infinite-dimensional dynamics of the moving average process based on a truncated Laguerre series expansion. The resulting problem is a finite-dimensional optimal stopping problem, which we propose to solve with a least squares Monte Carlo approach. We analyze the theoretical convergence rate of our method and present numerical results in the Black-Scholes framework.

Introduction

We introduce a new method to value American options whose payoff at exercise depends on the moving average of the underlying asset price. The simplest example (sometimes known as surge option) is a variable strike call or put, whose strike is adjusted daily to the moving average of the underlying asset over a certain fixed-length period preceding the current date. American-style options on moving average are widely used in energy markets. In gas markets, for example, these options are known as indexed Swing options and allow the holder to purchase an amount of gas at a strike price, which is indexed on moving averages of various oil-prices: typically gas oil and fuel oil prices are averaged over the last 6 months and delayed in time with a 1 month lag.

We shall denote by X the moving average of the underlying S over a time window with fixed length δ > 0:

X t = 1 δ t t-δ S u du, ∀t ≥ δ. (1)
moving average and additional partial averages of the price over the rolling period. Their number is computed heuristically and as it tends to the number of time steps within the rolling period, the computed price tends to the true price of the moving average option. The same kind of approach is used by Grau [START_REF] Grau | Applications of Least-Squares Regressions to Pricing and Hedging of Financial Derivatives[END_REF], but the author improves its numerical efficiency by a different choice of basis functions in the regressions used for the conditional expectations estimation. Kao and Lyuu [START_REF] Kao | Pricing of moving-average-type Options with Applications[END_REF] introduce a tree method based on the CRR model to price moving average lookback and reset options. Their method can handle only short averaging windows: the numerical results that are shown deal at most with 5 discrete observations in the averaging period. Indeed, this tree-based approach leads to an algorithm complexity (number of tree nodes) which exponentially increases with the number of time steps in the averaging period. Finally, Dai et al. [START_REF] Dai | A lattice algorithm for pricing moving average barrier options[END_REF] introduce a lattice algorithm for pricing Bermudan moving average barrier options. The authors propose a finite dimensional PDE model for such options and solve it using a grid method.

The pricing of moving average options is closely related to high-dimensional optimal stopping problems. It is well-known that deterministic techniques such as finite differences or approximating trees are made inefficient by the so-called curse of dimensionality. Only Monte Carlo type techniques can handle American options in high dimensions. Bouchard and Warin [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF] and references therein shall give to the interested reader a recent review of this research field.

More generaly, a related problem is that of optimal stopping of stochastic differential equations with delay. With the exception of a few cases where explicit dimension reduction is possible [START_REF] Federico | Optimal Stopping of Stochastic Differential Equations with Delay Driven by Lévy Noise[END_REF][START_REF] Gapeev | Optimal stopping problem in a diffusion-type model with delay[END_REF], there is no numerical method for solving such problems, and the Laguerre approximation approach of the present paper is a promising direction for further research.

The rest of the paper is structured as follows. In Section 2, we introduce the mathematical context and provide a general result which links the strong error of approximating one moving average process with another to a certain distance between their weighting measures. We then introduce an approximation of the weighting measure as a series of Laguerre functions truncated at n terms, which leads to (n + 1)-dimensional Markovian approximation to the initial infinite dimensional problem. The properties of Laguerre functions combined with our strong approximation result then enable us to establish a bound on the pricing error introduced by our approach as n goes to infinity. In Section 3, our numerical method, based on least squares Monte Carlo algorithm, is presented. The final section of this paper reports the results of numerical experiments in the Black-Scholes framework which include pricing moving average options with time delay.

Throughout the paper we assume that the price of the underlying asset S = (S t ) t≥0 is a nonnegative continuous Markov process defined on the probability space (Ω, F, F, P), where P is a martingale probability for the financial market and F = (F t ) t≤T is the natural filtration of S.

For the sake of simplicity, we present our results in the framework of a 1-dimensional price model but they are directly generalizable to a multi-asset model or to a model with unobservable risk factors such as stochastic volatility.

As usual, we denote by L 2 ≡ L 2 ([0, +∞)) the Lebesgue space of real-valued square-integrable functions f on [0, +∞) endowed with its norm:

f 2 := ∞ 0 |f (x)| 2 dx 1 2
.

A finite dimensional approximation of moving average options price

Strong approximations of moving average processes Consider a general moving average process of the form1

M t = ∞ 0 S t-u µ(du)
where µ is a finite possibly signed measure on [0, ∞). Throughout the paper, we shall adopt the following convention for the values of S on the negative time-axis:

S t = S 0 , ∀t ≤ 0. ( 2
)
We shall use an integrability assumption on the modulus of continuity of the price process: there exists a constant C < ∞ such that

E sup t,s∈[0,T ]:|t-s|≤h |S t -S s | ≤ Cε(h), ε(h) := h ln 2T h . (3)
Fischer and Nappo [START_REF] Fischer | On the moments of the modulus of continuity of Itô processes[END_REF] show that this holds in particular when S is a continuous Itô process of the form

S t = S 0 + t 0 b s ds + t 0 σ s dW s with E sup 0≤t≤T |b s | < ∞ and E sup 0≤t≤T |σ s | 1+γ < ∞
for some γ > 0.

The following lemma provides a tool for comparing moving averages with different weighting measures.

Lemma 2.1. Let Assumption (3) be satisfied, and let µ and ν be finite signed measures on [0, ∞) with Jordan decompositions µ = µ + -µ -and ν = ν + -ν -, such that µ + (R + ) > 0. Define

M t = ∞ 0 S t-u µ(du), N t = ∞ 0 S t-u ν(du). Then E sup 0≤t≤T |M t -N t | ≤ C|µ(R + ) -ν(R + )| + C µ + ([0, T ]) + ν -([0, T ]) + |µ([0, T ]) -ν([0, T ])| ε 1 µ + ([0, T ]) T 0 |F µ (t) -F ν (t)|dt (4)
for some constant C < ∞ which does not depend on µ and ν, where

F ν (t) := ν([0, t]) and F µ (t) := µ([0, t]).
Proof.

Step 1. We first assume that µ and ν are probability measures. Let F -1 µ and F -1 ν be generalized inverses of µ and ν respectively. Then,

E sup 0≤t≤T |M t -N t | = E sup 0≤t≤T 1 0 |S t-F -1 µ (u) -S t-F -1 ν (u) |du ≤ 1 0 E sup 0≤t≤T |S t-F -1 µ (u) -S t-F -1 ν (u) | du ≤ C 1 0 ε |F -1 µ (u) ∧ T -F -1 ν (u) ∧ T | du ≤ Cε 1 0 |F -1 µ (u) ∧ T -F -1 ν (u) ∧ T |du ,
where the last inequality follows from the concavity of ε(h). The expression inside the brackets is the Wasserstein distance between the measures µ and ν truncated at T . Therefore, from the Kantorovich-Rubinstein theorem we deduce

E sup 0≤t≤T |M t -N t | ≤ Cε T 0 |F µ (t) -F ν (t)|dt . Step 2. Introduce μ = µ1 [0,T ] and ν = ν1 [0,T ] + (µ([0, T ]) -ν([0, T ]))δ 2T
, where δ 2T is the point mass at the point 2T . Then,

E sup 0≤t≤T |M t -N t | ≤ C|µ(R + ) -ν(R + )| + E sup 0≤t≤T ∞ 0 S t-u μ(du) - ∞ 0 S t-u ν(du) ≤ C|µ(R + ) -ν(R + )| + (μ + (R + ) + ν-(R + ))E sup 0≤t≤T ∞ 0 S t-u μ+ (du) + ν-(du) μ+ (R + ) + ν-(R + ) - ∞ 0 S t-u μ-(du) + ν+ (du) μ+ (R + ) + ν-(R + ) .
Since μ+ (R + ) + ν-(R + ) = μ-(R + ) + ν+ (R + ), both measures under the integral sign are probability measures, and we can apply Step 1, which gives Introducing Laguerre approximation The aim of this paragraph is to provide heuristic arguments which lead to Laguerre approximation of the moving average. We would like to find a finitedimensional approximation to M , that is, find n processes Y 1 , . . . , Y n such that (S, Y 1 , . . . , Y n ) are jointly Markov, and M t is approximated in some sense to be made precise later by M n t which depends deterministically on S t , Y 1 t , . . . , Y n t . Since M is linear in S, it is natural to require that the approximation also be linear. Therefore, we assume that Y = (Y 1 , . . . , Y n ) satisfies the linear SDE

E sup 0≤t≤T |M t -N t | ≤ C|µ(R + ) -ν(R + )| + C(μ + (R + ) + ν-(R + ))ε 1 μ+ (R + ) + ν-(R + ) T 0 |F μ+ +ν -(t) -F μ-+ν + (t)|dt = C|µ(R + ) -ν(R + )| + C(μ + (R + ) + ν-(R + ))ε 1 μ+ (R + ) + ν-(R + ) T 0 |F µ (t) -F ν (t)|dt ,
dY t = -AY dt + 1(αS t dt + βdS t ), ( 5
)
where A is an n × n matrix, 1 is a n-dimensional vector with all components equal to 1 and α and β are constants. Similarly, the approximation is given by a linear combination of the components of Y : M n = B ⊥ Y , where B is a vector of size n and ⊥ denotes the matrix transposition. The solution to (5) can be written as

Y t = e -At Y 0 + t 0 e -A(t-s) 1(αS s ds + βdS s )
or, assuming stationarity, as

Y t = t -∞
e -A(t-s) 1(αS s ds + βdS s ) and

M n t = t -∞
B ⊥ e -A(t-s) 1(αS s ds + βdS s ).

Integration by parts then yields:

M n t = βB ⊥ 1S t + t -∞ B ⊥ (α -Aβ)e -A(t-s) 1S s ds := K n S t + t -∞ h n (t -u)S u du,
Recalling the structure of the matrix exponential, it follows that the function h n is of the form

h n (t) = K k=1 e -p k t n k i=0 c k i t i , (6)
where n 1 + . . . + n K + K = n (K is the number of Jordan blocks of A). Therefore, the problem of finding a finite-dimensional approximation for M boils down to finding an approximation of the form K n δ 0 (dt) + h n (t)dt for the measure ν. This problem is well known in signal processing, where the density h of ν is called impulse response function of a system, and h n is called Hankel approximation of h. For arbitrary µ and n, Hankel approximations may be very hard to find, and in this paper we shall focus on a subclass for which K = 1, that is, the function h n is of the form

h n (t) = e -pt n-1 i=0 c i t i . (7)
This is known as Laguerre approximation, because for a fixed p, the first n scaled Laguerre functions (defined below) form an orthonormal basis of the space of all functions of the form (7) endowed with the scalar product of L 2 ([0, ∞)) which will be denoted by •, • . See [START_REF] Wahlberg | On approximation of stable linear dynamical systems using Laguerre and Kautz functions[END_REF] for a discussion of optimality of Laguerre approximations among all approximations of type (6).

Definition 2.1. Fix a scale parameter p > 0. The scaled Laguerre functions (L p k ) k≥0 are defined on [0, +∞) by

L p k (t) = 2p P k (2pt)e -pt , ∀k ≥ 0 (8)
in which (P k ) k≥0 is the family of Laguerre polynomials explicitly defined on [0, +∞) by

P k (t) = k i=0 k k -i (-t) i i! , ∀k ≥ 0 (9)
or recursively by

     P 0 (t) = 1 P 1 (t) = 1 -t P k+1 (t) = 1 k+1 ((2k + 1 -t)P k (t) -kP k-1 (t)) , ∀k ≥ 1. ( 10
)
The scaled Laguerre functions (L p k ) k≥0 form an orthonormal basis of the Hilbert space L 2 ([0, ∞)), i.e.

∀(j, k), L p j , L p k = δ j,k . Fix now an order n ≥ 1 of truncation of the series. In view of Lemma 2.1, we propose the following Laguerre approximation of the moving average process M :

• Let H(x) = µ([x, +∞)).
• Compute the Laguerre coefficients of the function H:

A p k = H, L p k . Set H p n (t) = n-1 k=0 A p k L p k (t) and h p n (t) = -d dt H p n (t).
In view of Lemma A.2, the function h p n can be written as

h p n = n-1 k=0 a p k L p k , a p k = pA p k + 2p n-1 i=k+1 A p i . (11)
• Approximate the moving average M with

M n,p t = (H(0) -H p n (0))S t + +∞ 0 h p n (u)S t-u du, ∀t ≥ 0. ( 12
)
Remark 2.1. The approximation proposed in [START_REF] Kao | Pricing of moving-average-type Options with Applications[END_REF] (and in particular the correction coefficient in front of S t ) is chosen so that the total mass of the weighting measure of the approximate moving average M n,p t is equal to the total mass of the weighting measure µ of the exact moving average. In particular, such an approximation becomes exact for a constant asset price S.

From definitions [START_REF] Grau | Applications of Least-Squares Regressions to Pricing and Hedging of Financial Derivatives[END_REF] and [START_REF] Kao | Pricing of moving-average-type Options with Applications[END_REF], it seems natural to introduce n random processes X p,0 , . . . , X p,n-1 defined by

X p,k t = +∞ 0 L p k (v)S t-v dv, ∀t ≥ 0, ∀k = 0, . . . , n -1. ( 13
)
They will be called Laguerre processes associated to the process S throughout this paper and are related to the moving average approximation by

M n,p t = (H(0) -H p n (0))S t + n-1 k=0 a p k X p,k t , ∀t ≥ 0. ( 14
)
Proposition 2.1. Let n ≥ 1 and p > 0. The (n+1)-dimensional process (S, X p,0 , X p,1 , . . . , X p,n-1 ) is Markovian. The n Laguerre processes follow the dynamics:

               dX p,0 t = √ 2pS t -pX p,0 t dt dX p,1 t = √ 2pS t -2pX p,0 t -pX p,1 t dt . . . dX p,n-1 t = √ 2pS t -2p n-2 k=0 X p,k t -pX p,n-1 t dt with initial values X p,k 0 = S 0 (-1) k √ 2p p , ∀k ≥ 0. ( 15
)
Proof. Immediate, from Equations ( 13) and ( 8) and properties A.1-(i) and A.1-(ii).

Convergence of the Laguerre approximation Proposition 2.2. Let Assumption (3) be satisfied, and suppose that the moving average process M is of the form

M t = K 0 S t + ∞ 0 S t-u h(u)du (16)
where K 0 is a constant and the function h has compact support, finite variation on R, is constant in the neighborhood of zero and is not a.e. negative on [0, T ]. Then the error of approximation (14) admits the bound 4 ).

E sup 0≤t≤T |M t -M n,p t | ≤ Cε(n - 3
Proof. We shall use Lemma 2.1. The measures µ and ν are defined by µ(dx) = K 0 δ 0 (dx) + h(x)dx and ν(dx) = (H(0) -H p n (0))δ 0 (dx) + h p n (x)dx. Therefore, these measures have the same mass, and the first term in estimate (4) disapears. In addition,

|µ([0, T ]) -ν([0, T ])| = |H(T ) -H p n (T )|,
which remains bounded by Lemma A.4. Let us show that ν -([0, T ]) is bounded as well. For this it is enough to prove that (H p n ) 2 is bounded on n. A straightforward computation using Lemma A.2 shows that

c p k = 2pH(0) -2p k-1 i=0 A p i -pA p k , (17)
where c p k := h, L p k are Laguerre coefficients of h. By definition of H p n and a p k in [START_REF] Grau | Applications of Least-Squares Regressions to Pricing and Hedging of Financial Derivatives[END_REF], this leads to

a p k = c p k - √ 2p [H(0) -H p n (0)] .
We have thus:

(H p n ) 2 2 = k≤n-1 |a p k | 2 ≤ 2 k≤n-1 |c p k | 2 + 2 k≤n-1 2p [H(0) -H p n (0)] 2 = O(n -1 2 )
by Lemma A.4 and using √ 2p [H(0) -H p n (0)] = c p n + pA p n issued from [START_REF] Mäkilä | Laguerre series approximation of infinite dimensional systems[END_REF]. Therefore, there exists a constant C < ∞, which does not depend on n, such that

E[ sup 0≤t≤T |M t -M n,p t |] ≤ Cε 1 T 0 h + (t)dt T 0 |H(t) -H p n (t)|dt .
By Cauchy-Schwartz inequality and Lemma A.4,

T 0 |H(t) -H p n (t)|dt ≤ √ T H -H p n 2 = √ T   k≥n |A p k | 2   1 2 = O(n -3 4 ),
from which the result follows using the properties of ε and the fact that h is not a.e. negative on [0, T ] (which means that

T 0 h + (t)dt > 0).

Approximating option prices

The price of the American option whose pay-off depends on the moving average M and the price of the underlying is given by

sup τ ∈T E [φ (S τ , M τ )]
where T is the set of F-stopping times and φ is the payoff function. It can then be approximated by the solution to

sup τ ∈T E [φ (S τ , M n,p τ )] . (18)
Proposition 2.3. Let Assumption (3) be satisfied, and suppose that the payoff function φ is Lipschitz in the second variable and that the moving average process M satisfies the assumptions of Proposition 2.2. Then the pricing error admits the bound

E pricing (n, p) := sup τ ∈T E [φ (S τ , M τ )] -sup τ ∈T E [φ (S τ , M n,p τ )] ≤ Cε(n -3 4 ).
where C > 0 is a constant independent of n.

Proof. We have first:

∀τ, E [φ (S τ , M τ )] = E [φ (S τ , M n,p τ )] + E [φ (S τ , M τ ) -φ (S τ , M n,p τ )] =⇒ sup τ E [φ (S τ , M τ )] = sup τ E [φ (S τ , M n,p τ )] + E [φ (S τ , M τ ) -φ (S τ , M n,p τ )] ≤ sup τ E [φ (S τ , M n,p τ )] + sup τ E [φ (S τ , M τ ) -φ (S τ , M n,p τ )]
In consequence,

sup τ E [φ (S τ , M τ )] -sup τ E [φ (S τ , M n,p τ )] ≤ sup τ E |φ (S τ , M τ ) -φ (S τ , M n,p τ )| .
By symmetry and (A2), we get

sup τ E [φ (S τ , M τ )] -sup τ E [φ (S τ , M n,p τ )] ≤ sup τ E |M τ -M n,p τ | ≤ E sup 0≤t≤T |M t -M n,p t |
and the result follows from Proposition 2.2.

Uniformly-weighted moving average The uniform weighting measure

µ(dx) = h(x)dx = 1 δ 1 [0,δ] dx (19)
satisfies the assumptions of Proposition 2.2. In particular,

H(x) = 1 δ (δ -x) + . From Lemma A.2, the Laguerre coefficients A δ,p k = H, L p k are related to the Laguerre coefficients of h, c δ,p k = h, L p k 2 , via A δ,p k = (-1) k √ 2p p - 1 p c δ,p k - 2 p k-1 i=0 (-1) k-i c δ,p i , (20)
and the coefficients of h can be computed from the values of Laguerre polynomials:

c δ,p n = √ 2p δp (1 -e -pδ P n (2pδ)) + 2 n k=1 (-1) k (1 -e -pδ P n-k (2pδ)) . (21)
Given the length of the averaging window δ > 0 and an order n ≥ 1 of approximation (number of Laguerre functions), we determine the optimal scale parameter p opt (δ, n) as

p opt (δ, n) = arg min p>0 H -H p n 2 = arg min p>0 δ 3 - n-1 k=0 A δ,p k 2 . ( 22
)
Finding an explicit formula to p opt (δ, n) does not seem to be possible, but finding a numerical solution is easy using the explicit expressions ( 20), ( 21) and [START_REF] Shreve | Stochastic Calculus for Finance: continuous time models[END_REF]. In addition, once p opt is computed for a couple (1, n), the scaling property of Laguerre functions [START_REF] Federico | Optimal Stopping of Stochastic Differential Equations with Delay Driven by Lévy Noise[END_REF] gives the value of p opt (δ, n), for any δ > 0:

p opt (δ, n) = p opt (1, n) δ .
Table 1 gives the values for p opt (1, n) for the first 10 values of n computed with an accuracy of 10 -3 . Figure 1 Proposition 2.2 also applies when the moving average is delayed by a fixed time lag l ≥ 0:

X t = 1 δ t-l t-l-δ S u du, ∀t ≥ δ + l. (23)
In this case, the weighting measure is µ(dx) = 1 δ 1 [l,l+δ] dx and H(x) = 1 δ {(δ + l -x) + -(l -x) + }. Figure 2 shows the approximation of H by H popt(n) n

for n = 1, 3, 5, 7 Laguerre basis functions with δ = 1 and l = 0.5 as well as the L 2 -error made as a function of n (in the same way as above, we numerically compute p opt (n) for minimizing the L 2 -error made by Laguerre approximation). In comparison to the previous case, it appears that the number of Laguerre basis functions necessary to approximate H is greater for an equivalent accuracy of the approximation: the error is less than 5% from n = 5. This is due to the fact that the density of the weighting measure has two points of discontinuity.

A Monte Carlo-based numerical method

In this section we present a numerical method for computing the solution to the approximate problem [START_REF] Wahlberg | On approximation of stable linear dynamical systems using Laguerre and Kautz functions[END_REF], which is a (n + 1)-dimensional optimal stopping time problem. For the sake of simplicity, we restrict ourselves to uniformly weighted moving averages. Since the dimension of the problem may be high, we use a Monte Carlo technique. Our numerical approach corresponds to the one from Longstaff and Schwartz [START_REF] Longstaff | Valuing American options by simulation: A simple Least Squares approach[END_REF] and the computation of conditional expectations is done with a regression based approach. In particular, we shall use the technique of adaptative local basis proposed by Bouchard and Warin [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF].

Forward simulation in discrete time

We compute the price of the discrete time version of the American option [START_REF] Wahlberg | On approximation of stable linear dynamical systems using Laguerre and Kautz functions[END_REF] in which the moving average X has been replaced by its approximation M n,popt defined in [START_REF] Lebedev | Special functions and their applications[END_REF] and the exercise is possible on an equidistant discrete time grid with N ≥ 1 time steps ∆t = T N : π = {t 0 = 0, t 1 , . . . , t N = T }. We assume that there are exactly N δ ≥ 1 time steps within the averaging window of length δ: N δ = δ ∆t = δ T N , and that the spot price S can be simulated on π either exactly or using the Euler scheme. We denote this simulated discrete time price by S π t0 = S 0 , S π t1 , . . . , S π t N , extend this definition to [0, T ] by

S π t = S π ti+1 , ∀t ∈ (t i , t i+1 ] (24)
and shall also apply convention (2) to S π . We define the discrete time version X π of the moving average process X by

X π ti = 1 δ ti ti-δ S π t dt = 1 N δ i j=i-N δ +1 S π tj , ∀t i ∈ π. (25)
Similarly, the discrete time versions of the Laguerre processes are defined by

X p,k,π ti = ti -∞ L p k (t i -v)S π t dt = i j=1 S π tj -S π tj-1 (i -j + 1)∆t c (i-j+1)∆t,p k + S 0 (-1) k √ 2p p , ∀t i ∈ π.
Backward resolution of the optimal stopping time problem The resolution is based on the well-known backward American dynamic programming principle. We adopt a Longstaff and Schwartz-style approach which consists in estimating the optimal exercise time (or equivalently the optimal cashflows generated by the optimal exercise rule) instead of focusing on the computation of the option value processes (as for example in Tsitsiklis and Van Roy [START_REF] Tsitsiklis | Regression Methods for Pricing Complex American-Style Options[END_REF]). Throughout the paper, the approach presented below will be called (Lag-LS). The optimal payoffs are evaluated using the approximate value of moving average X π derived from ( 14):

M n,popt,π ti = (H(0) -H popt n (0))S π ti + n-1 k=0 a popt k X popt,k,π ti , ∀t i ∈ π. ( 26
)
Denote by (τ π i ) i=N δ ,...,N the sequence of discretized optimal exercise times: τ π i is the optimal exercise time after t i ∈ π. The backward algorithm works as follows:

1. Initialization: τ π N = T 2. Backward induction for i = N -1, . . . , N δ :

τ π i = t i 1 Ai + τ π i+1 1 Ai with A i = φ S π ti , M n,popt,π ti ≥ E ti φ S π τ π i+1 , M n,popt,π τ π i+1
3. Estimation of the option price at time 0:

V π 0 = E φ S π τ π N δ , M n,popt,π τ π N δ
in which:

E ti [•] = E •| S π ti , X popt,0,π ti , . . . , X popt,n-1,π ti .
Estimators of the conditional expectations are constructed with a Monte-Carlo based technique. It consists in using M ≥ 1 simulated paths on π of the (n + 1)-dimensional state process:

S π,(m) , X popt,0,π,(m) , . . . , X popt,n-1,π,(m) , m = 1, . . . , M.

The corresponding paths of the approximate moving average are denoted by: M n,popt,π,(m) , m = 1, . . . , M.

Conditional expectations estimators E M ti are then computed by regression on local basis functions (see the precise description of the procedure in Bouchard and Warin [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF]). We shall denote by b S , b X 0 , . . . , b X n-1 the numbers of basis functions used in each direction of the state variable: b S for S π , b X 0 for X popt,0,π , b X 1 for X popt,1,π , etc. The Monte-Carlo based backward procedure becomes thus:

1. Initialization: τ π,(m) N = T , m = 1, . . . , M 2. Backward induction for i = N -1, . . . , N δ , m = 1, . . . , M :    τ π,(m) i = t i 1 A (m) i + τ π,(m) i+1 1 A (m) i A (m) i = φ S π,(m) ti , M n,popt,π,(m) ti ≥ E M ti φ S π τ π i+1 , M n,popt,π τ π i+1
3. Estimation of the option price at time 0:

V π 0 = 1 M M m=1 φ S π,(m) τ π,(m) N δ , M n,popt,π,(m) τ π,(m) N δ
Remark 3.1. We will use a numerical improvement to this standard backward induction algorithm, which might seem rather natural for practitioners. It consists in evaluating the optimal payoffs using the exact value (25) of the moving average. In particular, the optimal stopping frontier becomes:

A * i = φ S π ti , X π ti ≥ E ti φ S π τ π i+1 , X π τ π i+1 .
This improved method will be called (Lag-LS*) and unlike (Lag-LS) will exhibit a monotone convergence as n goes to infinity.

"Non Markovian" approximation for moving average options Motivated by a reduction of dimensionality, the numerical method that is most often used in practice to value moving average options consists in computing the conditional expectations in the Longstaff-Schwartz algorithm using only the explanatory variables (S, X): namely, the price and the moving average appearing in the option payoff. The resulting exercise time is thus suboptimal, but the approximate option price is often close to the true price. To assess the improvement offered by our method, we systematically compare our approximation to this suboptimal approximate price, also computed using a Longstaff and Schwartz approach and referred to as (NM-LS). Let (θ π i ) i=N δ ,...,N denote the discrete time sequence of the estimated optimal exercise times (θ π i being the optimal exercise time after t i ∈ π). (NM-LS) algorithm works as follows:

1. Initialization: θ π N = T 2. Backward induction for i = N -1, . . . , N δ :

θ π i = t i 1 Ai + θ π i+1 1 Ai with A i = φ S π ti , X π ti ≥ E φ S π θ π i+1 , X π θ π i+1 | S π ti , X π ti 3.
Estimation of the option price at time 0:

U π 0 = E φ S π θ π N δ , X π θ π N δ
Similarly to other methods, the conditional expectations are computed with the adaptative local basis regression-based technique from [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF]. The numbers of basis functions used in each direction will be denoted by b S for S π and b X for X π .

Numerical examples

For our examples, we use the single-asset Black and Scholes framework. We study standard moving average options for different values of the averaging window δ as well as moving average options with delay [START_REF] Tanguy | Improved Method for Optimum Choice of Free Parameter in Orthogonal Approximations[END_REF].

With the same notations as in Section 3, recall that the dimension of the discrete time version of moving average option pricing problem is equal to N δ with a Markovian state:

S π ti , S π ti-1 , . . . , S π t i-N δ+1 , ∀t i ∈ π, t i ≥ t N δ .
We use the standard Longstaff and Schwartz algorithm for such a Bermudan option in dimension N δ as the benchmark method. This method will be called (M-LS) and our Monte Carlo regression based approach (see more details in [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF]) allows to deal with cases up to dimension 8. For applications in which N δ is larger, this method becomes computationally unfeasible. We provide at the end of this section a numerical comparison between the convergence rate of our Laguerre-based approximation and (M-LS) with respect to the state dimension.

Moving average options: benchmark prices Consider a standard moving average American option with value at time 0:

sup τ ∈T [δ,T ] E e -rτ φ (S τ , X τ ) , X τ = 1 δ τ τ -δ S u du
where the asset price S is assumed to follow the risk-neutral Black and Scholes dynamics:

dS t = S t (rdt + σdW t ) , S 0 = s
and W is a standard Brownian motion. We shall consider call options with pay-off φ(s, x) = (s-x) + . Unless specified otherwise, the following parameters are used below: Maturity T = 0.2 Risk free interest r = 5% Volatility σ = 30% Initial spot value s = 100 and we consider a Bermudan option with exercise possible every day (when T = 0.2, the time interval [0, T ] is divided into N = 50 time steps).

Table 2 shows the prices of moving average call options computed by (NM-LS) and (M-LS) for various averaging periods δ, with M = 10 million of Monte Carlo paths and b S = b X = 2. The prices are averages over 5 valuations and the relative standard deviation is given in brackets. For reasonable volatility coefficients of the underlying price process and relatively small averaging window δ, (NM-LS) seems to provide a very good approximation (from below) to moving average options prices. This justifies the approximation made by practitioners and (among others) by Broadie and Cao [START_REF] Broadie | Improved Lower and Upper Bound Algorithms for Pricing American Options by Simulation[END_REF]. Moving average American options: Laguerre approximation Figure 3 shows a simulated trajectory of the underlying price S π , its moving average X π with δ = 0.04 and the corresponding Laguerre-based moving average approximation M n,popt,π with n = 1, 3 and 7 Laguerre basis functions. Already for n ≥ 3 Laguerre basis functions, M n,popt,π accurately mimics the exact moving average dynamics of X π and this approximation seems to be almost exact when n = 7. Table 3 reports the prices of moving average call options computed using the Laguerre-based method presented in Section 3 (Lag-LS) and its improved version (Lag-LS*) (with the same parameters as above, in particular δ = 0.04). The price values are means over 5 valuations, the relative standard deviation is given in brackets and we used M = 5 million Monte Carlo paths for n = 1, . . . Remark 4.1. When the averaging window is large, the variance of the Laguerre states (X p,k k ) k≥0 is small, and at least much smaller than the variance of the price S. In consequence, increasing the numbers b X k of basis functions in the directions of these states does not have a strong impact on the conditional expectation estimators and the resulting option price. On the contrary, the number b S of basis functions in the direction of the spot price S should be sufficiently large.

N δ (NM-LS) (M-
Whereas (Lag-LS) oscillates as n increases (this is due to the non monotone approximation of the moving average X by M n,popt ), (Lag-LS*) shows a monotone convergence when increasing n, as shown in Figure 4. The limiting value (almost 4.277) is around 0.2% above the value computed by the practitioner's approximation (NM-LS). Figure 5 presents the prices of moving average call options computed by (Lag-LS*) and (NM-LS) when varying δ from 0 to T with the same parameters as above. 7 Laguerre basis functions were used with method (Lag-LS) as soon as N δ ≥ 8. For smaller N δ , we take n = N δ -1: n must satisfy the condition n ≤ N δ -1 because otherwise the estimation of the conditional expectation at time t N δ leads to a degenerate linear system. In the limit case when δ = T , we retrieve the price of the Asian option with payoff S T -1 T T 0 S t dt + . For large averaging periods, the price that we obtain with 7 Laguerre functions is about 0.30% above the benchmark value given by the (suboptimal) non Markovian approximation. Moving average American options with time delay Consider now a moving average American option with time delay l ≥ 0 whose value at time 0 is:

sup τ ∈T [δ+l,T ] E [φ (S τ , X τ )] , X τ = 1 δ τ -l τ -l-δ S u du.
With the same option characteristics and parameters as above and an averaging period equal to δ = 0.02 (number of time steps N δ = 5), Figure 6 presents the prices of delayed moving average call options computed by (Lag-LS*) and (NM-LS) when varying l from 0 to T -δ. In the limit case when l = T -δ we retrieve the price of the Asian option with payoff

S T -1 T -l T -l 0 S t dt + .
The relative difference between the option values given by (Lag-LS*) and (NM-LS) is significant (bigger than 5%) for time lags such that l ∈ [0.04, 0.152] (corresponding to 10 ≤ N l ≤ 38). For example, when l = 0.1 (corresponding to N l = 25 time steps), the relative difference is around 11%. As expected, for pricing delayed moving average options, the non Markovian approximate method (NM-LS) seems to be a worse approximation than in the case without time lag: the error increases with the time lag but decreases with the length of the the averaging period. For large time lags and relatively small averaging period, our improved Laguerre-based method (Lag-LS*) gives option prices up to 11% above the benchmark value given by the suboptimal approximation (NM-LS) (cf. N l = 25 on Figure 6). For a good accuracy of (Lag-LS*), the required number of Laguerre functions is however bigger than in the case without time lag as explained at the end of Section 2.

Pricing of moving average Bermudean options: a convergence rate improvement To compute the exact price of a moving average Bermudan option, one can either use the classical method (M-LS) taking a sufficient number of steps within the averaging window or use the Laguerrebased method with sufficient number of state processes. In this last example we illustrate the fact that our method (Lag-LS*) converges much faster than classical method (M-LS) with respect to the state dimension for pricing the same Bermudan option.

Let us consider a moving average call option with maturity T = 0.5 and moving window δ = 0.1. Figure 8 provides a comparison between price values given by (Lag-LS*) for a time step ∆t = 1 80 when varying the number of Laguerre functions from 1 to 7 and by (M-LS) when varying the number of time steps within the averaging period from 2 to 8, that is ∆t = 1 20 , 1 30 , . . . , 1 80 (the state dimension varies in both cases from 2 to 8). M = 20 million of Monte Carlo paths were used in both cases and b S = 2, b X = 1.

x 0

I n (x )dx = 4 ∞ k=0 A k x ν (k+2)/2 J k+2 ( √ νx), (31)
where A 0 = 1, A 1 = 0, A 2 = 1 2 and other A i -s satisfy the equation (m + 2)A m+2 = (m + 1)A mν 2 A m-1 . The series converge uniformly in x on any compact interval.

Proof. The first formula is from [START_REF] Erdélyi | Higher transcendental functions[END_REF]. The other two follow readily using the integration formula for Bessel functions: In addition, for x > 0 fixed, e -x/2 P n (x) = O(n -1/4 ).

Proof. This result follows from Lemma A.3, using the asymptotic expansion for Bessel functions

J n (x) = ( 1 2 πx) -1/2 cos x - π 2 n - π 4 + O(x -3/2 ),
which holds uniformly [START_REF] Erdélyi | Higher transcendental functions[END_REF] on bounded domains outside a neighborhood of zero.

  because μ coincides with µ and ν coincides with ν on [0, T ]. Using the properties of the function ε and the definition of μ and ν, we then get (4) with a different constant C.

  (left graph) illustrates the approximation of H by the truncated Laguerre expansion H popt(n) n for n = 1, 3, 7 Laguerre basis functions (with δ = 1). The corresponding error ||H -H popt(n) n || 2 as a function of n is shown in the right graph. The error is less than 5% already with n = 3. A simple least squares estimation by a power function gives a behavior in O(n -1.06

Figure 1 :

 1 Figure 1: Left: Laguerre approximation of the function H(x) = 1 δ (δ -x) + . Right: L 2 error of the approximation.

Figure 2 :

 2 Figure 2: Left: Laguerre approximation of the function H(x) = 1 δ {(δ + l -x) + -(l -x) + }. Right: L 2 error of the approximation.

7 Figure 3 :

 73 Figure 3: Simulated trajectory of the moving average process and its Laguerre approximations.

, 3

 3 Laguerre functions and M = 10 million Monte Carlo paths for n = 4, . . . , 7 Laguerre functions, with b S = 4 and b X k = 1, ∀k ≥ 0. With M = 10 million Monte Carlo paths, b S = 4 and b X = 1, (NM-LS) gives an option value equal to 4.268.

Figure 4 :

 4 Figure 4: Convergence of the improved Laguerre-based approximation.

Figure 5 :

 5 Figure 5: Moving average option price as function of the averaging window δ.

Figure 6 :

 6 Figure 6: Price of the moving average option with time delay as function of the lag l.

  LS*) with n = 7 (Lag-LS*) with n = 5 (Lag-LS*) with n = 3 (Lag-LS*) with n = 2

Figure 7 :

 7 Figure 7: Price of the moving average option with time delay as function of the averaging window δ.

Figure 8 :

 8 Figure 8: Convergence of the improved Laguerre-based approximation and the benchmark method for pricing a Bermudan option.

1 0x

 1 ν+1 J ν (ax)dx = a -1 J ν+1 (a). Lemma A.4. Let µ be a finite signed measure on [0, ∞), with bounded support which does not contain zero. Let {c n } denote the Laguerre coefficients of the function h(x) := µ([x, ∞)) and {A n } denote the Laguerre coefficients of the function H(x) := ∞ x h(t)dt. Then c n = O(n -3/4 ) and A n = O(n -5/4 ).

Table 1 :

 1 ). opt (1, n) 2.149 4.072 6.002 4.234 5.828 7.473 9.155 10.866 9.153 10.726 Optimal scaling parameters for approximating H(x) = 1

	n	1	2	3	4	5	6	7	8	9	10
	p

δ (δ -x) + .

Table 2 :

 2 Moving average options pricing with (NM-LS) and (M-LS).

	LS)

Table 3 :

 3 Moving average options pricing with (Lag-LS) and (Lag-LS*).

	n	(Lag-LS*)	(Lag-LS)
	1 4.266 (0.020 %) 4.092 (0.017 %)
	2 4.273 (0.022 %) 4.302 (0.019 %)
	3 4.276 (0.023 %) 4.182 (0.018 %)
	4 4.276 (0.022 %) 4.227 (0.020 %)
	5 4.277 (0.023 %) 4.275 (0.020 %)
	6 4.277 (0.024 %) 4.287 (0.022 %)
	7 4.277 (0.024 %) 4.258 (0.022 %)

In the literature (see[START_REF] Basse | Lévy driven moving averages and semimartingales[END_REF] and references therein), moving averages are usually defined via the stochastic integral of S. Our definition as an ordinary integral with respect to a weighting measure is closer to the financial specifications.

A Appendix

Lemma A.1. The Laguerre polynomials (P k ) k≥0 belong to C ∞ ([0, +∞)) and:

Proof. (i) can be found for example in Szegö [START_REF] Szegö | Orthogonal polynomials[END_REF] and (ii) is a consequence of (10).

Lemma A.2. The definite integrals and derivatives of Laguerre functions can be computed using the following formulas:

∞ t e -s/2 P n (s)ds = 2e -t/2 P n (t) + 4e -t/2 n k=1 (-1) k P n-k (t). ( 27)

Proof. This follows, after some computations, from the contour integral representation of Laguerre polynomials:

e -ts 1-s (1 -s)s n+1 ds.