Cerea
email: korsakissok@cerea.enpc.fr

This paper details a number of existing formulations used in Gaussian models in a clear and usable way, and provides a comparison within a single framework-the Gaussian plume and puff models of the air quality modeling system Polyphemus. The emphasis is made on the comparison between (1) the parameterizations to compute the standard deviations and (2) the plume rise schemes. The Gaussian formulae are first described and theoretically compared. Their evaluation is then ensured by comparison to the observations as well as to several well-known Gaussian and CFD models performance. The model results compare well with the other Gaussian models for two of the three parameterizations for standard deviations, Briggs's and similarity-theory, while Doury's shows a tendency to underestimate the concentrations because of a large horizontal spread. The results with Kincaid experiment point out the sensitivity to the plume rise scheme and the importance of an accurate modeling of the plume interactions with the inversion layer. Using three parameterizations for the standard deviations and the same number of plume rise schemes, we were able to highlight a large variability in the model outputs.

Introduction

Gaussian models are widely used to model dispersion at local scale, in spite of their wellknown limits, since they are based on simple analytical formulae, and computationally cheap to use. However, the use of empirical schemes to compute the Gaussian standard deviations involves large uncertainties-hence many different formulations exist. When plume rise occurs, additional parameterizations are introduced, which results in still more variability. In 1 this paper, the Gaussian models developed on the Polyphemus platform [START_REF] Mallet | Technical Note: The air quality modeling system Polyphemus[END_REF] are used to evaluate the results sensitivity to the physical parameterizations for standard deviations and plume rise. This evaluation is ensured by means of theoretical study, comparison to experimental data, and comparison to other well-known Gaussian models performance, like US EPA models ISCST3 and AERMOD, and UK's widely used Gaussian model ADMS.

Since Polyphemus is not a single air quality model but a platform aimed at handling different models, it provides an appropriate framework to benchmark parameterizations and features of Gaussian models, thus evaluating the spread in model outputs. In addition, one of the aims of this paper is also to put together a number of different parameterizations that were previously described in many different contexts, in a clear and usable way, which also makes comparisons easier.

We first present the Polyphemus Gaussian models in Section 2, with an emphasis on the parameterizations for standard deviations and plume rise. A first comparison of the parameterizations is given in Section 3. Then, the model evaluation is carried out with Prairie Grass (Section 4) and Kincaid (Section 5) field-experiments datasets.

Description of the Polyphemus Gaussian models

a. The Gaussian plume and puff models Polyphemus hosts both a Gaussian plume model and a Gaussian puff model. [START_REF] Arya | Air pollution meteorology and dispersion[END_REF] reviews the main assumptions for the relevance of the two models. A Gaussian plume model assumes a Gaussian distribution of mean concentration in the horizontal (crosswind) and vertical directions, in steady-state and homogeneous meteorological conditions. The dispersion in the downwind direction is supposed to be negligible compared to the transport by the mean wind in the plume model. A Gaussian puff model discretizes the plume into a series of puffs, each of them having a Gaussian shape in all three directions. The puff model can also take into account different meteorological data for each puff, therefore allowing a more accurate simulation of non-stationary and non-homogeneous conditions. In Polyphemus, the plume and puff models share the same parameterizations. Thus, in stationary and homogeneous conditions-as assumed in the Prairie Grass and Kincaid experiments-the Gaussian puff model gives the same solution as the Gaussian plume model, if averaged over a long enough period. Henceforth, only the Gaussian plume model is therefore used (it was verified that the Gaussian puff model gave the same results).

For a stationary plume, the concentration C is given by the Gaussian plume formula:

C(x, y, z) = Q 2πσ y σ z ū exp - (y -y s ) 2 2σ 2 y × exp(- (z -z p ) 2 2σ 2 z ) + exp(- (z + z p ) 2 2σ 2 z ) . (1)
Q is the source emission rate, given in mass per second, ū is the mean wind velocity, and σ y and σ z are the Gaussian plume standard deviations in the horizontal (crosswind) and vertical directions. The coordinate y refers to the crosswind horizontal direction, and y s is the source coordinate in that direction. The coordinate z refers to the vertical coordinate and z p is the plume height above ground-sum of the source height and the plume rise (see Section c). The second term of the last factor represents the ground reflection; similarly, additional terms can take into account the reflections on the elevated inversion layer. Also, when the plume is mixed enough in the boundary layer (σ z > 1.5 z i ), the formula is modified to consider that the plume is vertically homogeneous. In addition, other factors can be applied to this formula in order to take loss processes into account, such as dry deposition, scavenging, radioactive or biological decay. In Polyphemus, full gaseous chemistry may also be handled by the Gaussian puff model.

b. Dispersion parameterization schemes

The expression of the total plume standard deviation is the sum of the spread due to turbulence (σ y turb , σ z turb ), the additional spread due to plume rise (σ ypr , σ zpr ), and the initial spread due to the diameter d s of the source:

σ 2 y = σ 2 y turb + σ 2 ypr + d 2 s /4, σ 2 z = σ 2 z turb + σ 2 zpr (2)
For the sake of clarity, σ y turb and σ z turb are named in this section σ y and σ z . These standard deviations are estimated through empirical schemes, as functions of the downwind distance to the source and the stability, based on a few dispersion field experiments. Here, three parameterizations are described. Two of them are based on a discrete description of the atmospheric boundary layer: the Briggs formulae and the Doury formulae. The third one is based on similarity theory, and involves parameters such as the wind velocity fluctuations, the Monin-Obukhov length, the mixing height and the friction velocity.

1) Briggs's formulae

The Briggs formulae are based on the Pasquill-Turner stability classes [START_REF] Turner | Workbook of atmospheric diffusion estimates[END_REF]) and fitted on the Prairie Grass experiment. This parameterization is born from an attempt to synthesize several widely used schemes by interpolating them for open country and for urban areas, for which they are particularly recommended. The full formulae can be found in [START_REF] Arya | Air pollution meteorology and dispersion[END_REF] for instance. The general form is given by:

σ y = α x √ 1 + β x , σ z = α x (1 + β x) γ , (3)
with x the downwind distance from source, and α, β and γ coefficients depending on Pasquill stability class. There are six stability classes corresponding to different states of the atmosphere [START_REF] Pasquill | The estimation of the dispersion of windborne material[END_REF], from A (extremely unstable) to F (stable).

2) Doury's formulae

An alternative parameterization is described in [START_REF] Doury | Une méthode de calcul pratique et générale pour la prévision numérique des pollutions véhiculées par l'atmosphère[END_REF]. This parameterization has been developed for the specific application of radionuclide dispersion, and fitted on a wider experimental field than the Prairie Grass field. The formulae use only two stability situations, corresponding to "normal" and "low" dispersion, determined by the vertical temperature gradient. By default, it is assumed that low dispersion occurs during nighttime with low wind speed (u ≤ 3 m s -1 ). The standard deviations are given in both cases in the general form:

σ y = (A h t) K h , σ z = (A z t) Kz , ( 4
)
where t is the transfer time since release time. In the case of a steady-state plume, t = x/u, where x is the distance from the source and u is the wind speed. The formulae for the normal situation can be found in [START_REF] Demaël | Comparative evaluation of an eulerian CFD and gaussian plume models based on prairie grass dispersion experiment[END_REF].

3) Similarity theory

If enough meteorological data is available, σ y and σ z can be estimated using the standard deviations of wind velocity fluctuations in crosswind horizontal direction σ v and in vertical direction σ w . Following [START_REF] Irwin | Scheme for estimating dispersion parameters as a function of release height[END_REF], dispersion coefficients are investigated in the form:

σ y = σ v t F y , σ z = σ w t F z , ( 5
)
where t is the time in seconds, and F y and F z are functions of a set of parameters that specify the characteristics of the atmospheric boundary layer. These functions are determined from experimental data. Various expressions of F y and F z have been proposed (e.g., [START_REF] Irwin | Scheme for estimating dispersion parameters as a function of release height[END_REF], [START_REF] Weil | Dispersion in the convective boundary layer[END_REF]). Appendix 7 details the formulations used in Polyphemus for F y , F z , σ v and σ w . A modified formulation from [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF] is also available for the case of elevated sources (above 100 m height). Two different formulations are also given for the standard deviations above the boundary layer.

c. Plume rise schemes

When emissions are hotter than the ambient air, or have a significant ejection speed, the plume generally reaches a certain height before behaving like the surrounding air: this difference in heights can be significant (up to a few hundreds of meters) and is called "plume rise" (see Figure 1). Estimates of the plume rise can be provided, based on the source initial height, diameter, ejection temperature and ejection velocity. In all formulae described below, it is assumed that the plume instantaneously reaches its final height at the emission time, then travels downwind. More accurate formulations allow to determine the plume rise as a function of the downwind distance, but are not investigated here.

If the plume rise is ∆h, the plume final height is given by z p = z s + ∆h, with z s the source height. The initial spread due to plume rise is added to the total standard deviations at the emission time according to Equation 2. It is estimated following [START_REF] Irwin | Scheme for estimating dispersion parameters as a function of release height[END_REF] for the crosswind spread and [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF] for the vertical spread:

σ ypr = ∆h/3.5, σ zpr = ∆h/2. ( 6
)
Three different formulae to compute plume rise are described in this section. While the Briggs-HPDM formulae depend on the stability of the atmosphere, the two other formulae depend only on the source dynamics and buoyancy and are based on the source heat rate.

1) Briggs-HPDM plume rise

The plume rise computation is based on the Briggs formulae-detailed in Seinfeld and Pandis (1998)-assuming a buoyancy-driven plume rise (Equations 8 to 12) . In the formulae, ∆h is the plume rise, u is the wind velocity, w * is the convective velocity and u * is the friction velocity. In addition, s p is the Briggs static stability parameter and F b is the initial buoyancy flux parameter. They are defined as:

s p = gdθ/dz T , F b = g v s d 2 s T s -T T s (7)
with g the gravity acceleration, dθ/dz the vertical gradient of potential temperature, T the ambient temperature, v s the vertical velocity of the emission, d s the source diameter and T s the source temperature.

(i) Stable cases

The formulae for stable cases come from [START_REF] Seinfeld | Atmospheric chemistry and physics: from air pollution to climate change[END_REF] and are also used in the HPDM model [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF] . The final plume rise is ∆h = min(∆h 1 , ∆h 2 ).

∆h 1 = 2.6

F b u sp 1/3 , ∆h 2 = 4 F 1/4 b s -3/8 p . (8)
(ii) Unstable and neutral cases

The following formulae come from [START_REF] Seinfeld | Atmospheric chemistry and physics: from air pollution to climate change[END_REF] and are originally from Briggs. ∆h = 21.4

F 3/4 b u if F b < 55 (9) ∆h = 38.71 F 3/5 b u if F b ≥ 55 (10)
In the above formulae, ambient turbulence is supposed to be negligible, compared to the plume internal turbulence. However, in many unstable and neutral cases, the plume eventually reaches a point where ambient turbulence is large enough to stop the plume vertical progress. In [START_REF] Briggs | Some recent analyses of plume rise observations[END_REF], this so-called "breakup" height is determined as the point where both the plume and the ambient eddy dissipation rates are equal. It is reported in [START_REF] Hanna | Applications in air pollution modeling[END_REF] and used here for neutral cases, whereas the unstable breakup formula comes from [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF]. When there is the choice between a breakup formula and a Briggs formula, the one giving the minimal plume rise is chosen. ∆h = 4.3

F b u w 2 * 3/5 z 2/5 i unstable (11) ∆h = 1.54 F b u u 2 * 2/3 z 1/3 s neutral (12)
2) Holland formula

This plume rise formula is the same for all stabilities. It comes originally from [START_REF] Holland | A meteorological survey of the oak ridge area[END_REF] and was modified by [START_REF] Stümke | Vorschlag einer empirischen Formel für die Schornsteinüberhöhung[END_REF].

∆h = 1.5 d s v s u + 65 d 1.5 s u T s -T T s 1/4 (13)
3) Concawe formula

In [START_REF] Brummage | The calculation of atmospheric dispersion from a stack[END_REF], a comparison study was made between several plume rise formulae, including the Holland-Stümke formula which was among the best. An attempt was made to synthesize the multiple formulations into a relationship in the form ∆h = K

Q α h u β , where Q h
is the source heat rate and K, α and β are empirical constants. It resulted in the Concawe formulae:

∆h = 0.071 Q 0.55 h u 0.67 (14)
Here, the ejection heat rate is computed using the source diameter, its ejection velocity and the temperature difference between the plume and ambient air:

Q h = 228.19 v s d 2 s (T s -T ) (15)
This is an approximation, based on the assumption that the calorific capacity of the ejected gas is not too different from that of air.

4) Partial penetration in the inversion layer

In the case of an elevated source with a large plume rise, the plume can partially or totally penetrate the inversion layer (Figure 2). Thus, the part of the plume which is above the inversion height is supposed to be trapped there and cannot re-enter the boundary layer.

Overlooking this phenomenon can lead to false estimations of the ground concentration, either by overestimating concentrations if the plume is always trapped below the inversion layer, or by underestimation if the plume is trapped above the inversion layer. A partial penetration of the plume is therefore assumed when the plume centerline (after plume rise) is close enough to the inversion height. The formulae come from [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF]. If the plume rise ∆h exceeds (z i -z s )/1.5, the top of the plume impacts the capping inversion, and a partial or total penetration of the inversion layer occurs. The fraction of the plume that has penetrated the inversion layer, called the penetration factor, is equal to P = 1.5-∆z/∆h if 0.5 < ∆z/∆h < 1.5 (partial penetration), and P = 1 if ∆z/∆h ≤ 0.5 (total penetration),

where ∆z = z i -z s . In the case of a partial penetration, the fraction of the plume that remains within the boundary layer has a new emission rate Q (1 -P ) and a new value of plume rise ∆h ′ = (0.62 + 0.38P )(z i -z s ) (shown in Figure 2). In the case of a total penetration, the remaining emission rate is equal to zero and there is no plume left in the boundary layer.

Comparison of the dispersion formulae a. Comparison set-up

The evolution of the Gaussian standard deviations (without a possible additional spread due to the diameter and/or plume rise) can be plotted for the parameterizations described in the previous section, in order to evaluate the spread. It must be noted that the Briggs These preliminary results highlight a large variability in the standard deviations, depending on the stability and the parameterization. Thus, a substantial variability can also be expected in the output results. Besides, because of the differences in estimations of σ z , the plume might penetrate the inversion layer with some parameterizations and not with others, which would increase the variability. It should also be pointed out that the stabil-ity diagnosis depends on the standard deviation formula used. For instance, in the Prairie Grass experiments, the Pasquill-Turner diagnosis used with the Briggs formulae considers 30% of the experiments to be unstable, 20% stable and 50% neutral, whereas the diagnosis based on Monin-Obukhov length (Table 8)-used for similarity theory-gives about 48% of unstable cases and 52% of stable cases. This might compensate some of the differences in the standard-deviation estimations: for instance, similarity theory in unstable cases and Briggs's formulae in neutral cases give similar estimations.

Evaluation with Prairie Grass experiment a. Model evaluation criteria

The statistical evaluation method is used for all models in this comparison, and described in [START_REF] Chang | Air quality model performance evaluation[END_REF]. It is also the method used in the Model Validation Kit [START_REF] Olesen | Consolidating tools for model evaluation[END_REF], which has been developed and tested during a series of workshops and conferences on Harmonization within Atmospheric Dispersion Modeling for Regulatory purpose.

In many field experiments, the sensors are placed in arcs at various distances downwind of the source, and the maximum predicted and observed concentrations for each arc ("maximum arcwise concentrations") are compared. Therefore, all other concentration values on the arc are not used in the comparison. The maximum observed value is determined on timeaveraged concentration measurements (the averaging time for Prairie Grass experiments is 10 minutes). The comparison study is carried out using maximum arcwise concentrations normalized by the emission rate Q-so-called "centerline values"-in µg m -3 /(g s -1 ).

The results are discussed using scatter diagrams as well as statistical performance measures such as the fractional bias (FB), the normalized mean square error (NMSE), the fraction of predictions within a factor two of observations (FAC2) and the correlation coefficient were also evaluated using Prairie Grass data.

Table 1 shows A more detailed analysis of the results with Polyphemus Gaussian models is provided below.

d. Statistics for each arc

The model performance is now analyzed for the five monitoring arcs, so as to compare the behavior of the three standard-deviation parameterizations. Figure 5 shows the scatter plots for the three parameterizations. While the behavior of the model with the Briggs parameterization is quite well fitted with the observations, the overprediction tendency at the two further arcs is clearly shown for similarity theory.

The scatter plot for the Doury parameterization confirms its poor performance especially in term of FAC2: the results for one arc are almost flat, which means that the modeled concentrations at one arc do not vary much with the different experiments. This could come from the representation of the atmosphere stability in only two classes (normal and low diffusion) which may be insufficient in some meteorological situations. Section 3 also showed that the Doury parameterization can give a very high value of σ y compared to the other two parameterizations at a given distance, especially at low wind speeds.

Evaluation with Kincaid experiment a. Experiment and modeling set-up

The Kincaid experiment took place at the Kincaid power plant in Illinois (USA) in 1980 and 1981. The plant stack is 187 m high with a diameter of 9 meters and is surrounded by flat farmland and some lakes. Meteorological data came from an instrumented tower, and the conditions during the experiments ranged from neutral to convective. The tracer material was a buoyant plume of SF 6 gas, released from the stack and measured by monitors deployed on arcs between 0.5 and 50 km from the source at ground level. There were 200 monitors that were shifted day-to-day to correspond to the wind direction, and the arcs distances also varied with the meteorological situation. The recorded measurements are the arc maximum concentrations averaged over 1 hour. The database contains 171 hours of measurements. This dataset is often used to validate plume rise parameterizations because of the highly buoyant release that led to substantial plume rise. All data used for the following simulations is given with the Model Validation Kit (http://www.harmo.org/kit/ Download.asp). In this dataset, a quality indicator was assigned to each measurement value, denoting how reliable the maximum concentration should be considered [START_REF] Olesen | The model validation exercise at mol: overview of results[END_REF]. Only the most reliable data, of quality 2 and 3, were used here. The statistics are computed on 165 experiments.

It has to be pointed out that in the following study, the simulations with Polyphemus Gaussian model use the observed inversion height when available and, for some cases, the predicted height. On the contrary, other models presented here like ADMS, AERMOD, and most likely OML and HPDM, recomputed the inversion height with their own preprocessor.

reaches very high values for low wind speeds. The Holland and Concawe plume rises are similar, which is not surprising considering that the former was used in the design of the latter. Hence, we can suppose that the model behavior will be similar for those two schemes, whereas using Briggs-HPDM could lead to a model underestimation, if the plume is predicted to be too high. In addition, Table 2 shows the proportion of cases in Kincaid experiments with total or partial penetration in the inversion layer: as expected, it is higher for Briggs-HPDM plume rise, with almost 60% of the cases. The Holland and Concawe plume rises give about the same number of cases with total penetration. However, there are less cases with partial penetration given by Holland than by Concawe formulae.

Hereafter, evaluation results are analyzed for each of the three plume-rise parameterizations (statistics on all arcs and experiments). For each plume rise parameterization, the scheme for standard deviations which performs best is in bold in the tables. The similaritytheory parameterization is labeled "Sim.th.".

2) Results with Briggs-HPDM plume rise

Table 3 shows the model results with Briggs-HPDM plume rise, for the three standarddeviation schemes. As expected, the tendency is to underpredict the concentrations, by a factor of two for the Briggs formulae and three for the Doury formulae. On the contrary, the similarity-theory parameterization performs quite well for all indicators (since there is nearly no bias), except for the correlation, which is quite low for all standard deviations. This is confirmed by Figure 7 where the scatter plot for the Briggs formulae clearly indicates the underestimation tendency (Figure 7(b)), whereas the one with similarity theory (Figure 7(a)) is much more centered. However, in both cases, one can observe "vanishing concentrations"-concentrations that are observed but are not reproduced by the model.

This can occur when the plume rise is too high, if the plume totally penetrates the inversion layer, or if it touches the ground later than the observed plume. In the first case, simulated concentrations are equal to zero at all arcs for the experiment, since there is no re-entering of the plume below the inversion layer. In the second case, only the concentrations at the nearest arcs to the source are missing, that is, the highest concentrations. This phenomenon could be responsible for the low correlations that are globally observed.

3) Results with Holland plume rise

Results with Holland plume rise are given in Table 4. The modeled concentrations are globally higher than with the HPDM plume rise, since the estimated plume rise is lower.

However, the Doury parameterization still underestimates considerably the concentrations, as with the Prairie Grass experiment. The bias and NMSE are considerably improved for the Briggs parameterization whereas the concentrations overestimation with similarity theory has increased, leading to a higher bias and NMSE than with HPDM plume rise.

The "vanishing concentrations" can be seen for the scatter plots of all plume rise parameterizations (not shown here), not only for Briggs-HPDM. Hence, one cannot point out that a particular scheme overestimates the plume rise, and this may also come from an inaccurate estimation of the observed boundary layer height.

4) Results with Concawe plume rise

Here, the best results are given by the Briggs parameterization, as in the case of the Holland formulae. For the Briggs-HPDM plume rise, the best results were with similarity theory. As expected, the results with Concawe plume rise are not very different from the Holland results, especially when looking at the mean concentrations and the biases. However, the correlation and the proportion of values within a factor two of the observation have increased. This slightly better performance is not surprising, considering that this formula was built from a set of different plume rise formulae, including Holland, and adjusted on a set of cases, especially neutral and convective cases with buoyant plumes, as explained in [START_REF] Brummage | The calculation of atmospheric dispersion from a stack[END_REF].

c. Comparison with other Gaussian models

For this comparison, we use the results obtained with the Concawe parameterization which gives the best results for all standard-deviation formulae. In addition, the results with similarity theory and HPDM plume rise are also provided since they are good despite the lower correlation. Results for the other models come from CERC (2007), as in the case of Prairie Grass, for quality 2 and 3 (Table 6). As noted in Section a, we used 165 experiments instead of 171, hence the slight differences in the observed mean and standard deviations.

Results with Polyphemus models are between those of AERMOD and ADMS4, except for the correlations that are lower as already observed, and comparable to these of ISCST3.

Using data of quality 3 only, we can also compare with HPDM (USA, [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF]), OML (Denmark), SAFE AIR (Italy, [START_REF] Canepa | About the plume rise description in the dispersion code SAFE AIR[END_REF]), as well as the Lagrangian model NAME (UK, [START_REF] Webster | Validation of a lagrangian model plume rise scheme using the kincaid data set[END_REF]). The statistics for HPDM and OML come originally from [START_REF] Olesen | The model validation exercise at mol: overview of results[END_REF]. Table 7 shows the results for all these models. ISCST3, AER-MOD and Polyphemus with the Doury formulae show a significant underestimation of the mean values, hence the high bias and NMSE. All other models slightly underestimate the mean values except SAFE AIR and Polyphemus with similarity theory and Concawe plume rise which slightly overpredicts the observations. The correlation results are surprisingly against this trend: for instance, AERMOD has a correlation of 40% which is comparable to the best results. On the contrary, models that have an otherwise good performance such as SAFE AIR, OML, and Polyphemus models with Concawe formulae have relatively low correlations. The correlation for Polyphemus with similarity theory and HPDM plume rise is very low (5%) but this configuration is among the best models for all other indicators. This correlation is not surprising since the HPDM model evaluation in [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF] also showed low correlations with this experiment. Further improvements were made in HPDM, such as a correction to take into account the plume lofting at the inversion layer during convective conditions, and this seems to have greatly improved the results. Table 7 also shows that the spread in the model outputs is very large: the predicted mean concentrations range from about 20 ng m -3 /g s -1 to 70 ng m -3 /g s -1 , which means more than a factor of 3 between some models. It is interesting to note that, with the different parameterizations available in Polyphemus for standard deviations and plume rise, we can reproduce the same spread, the model configuration being otherwise unchanged. Only the correlations are still low, no matter the formula used, which tends to prove that it comes from the shared inversion-height modeling (plume penetration and lofting).

Conclusion

Polyphemus Gaussian plume and puff models were described, with an emphasis on the parameterizations used to estimate the standard deviations and the plume rise. The Gaussian plume model was then used to compare these various parameterizations within a single framework, all other model features being equal. A preliminary study on the standard deviation values (Section 3) highlighted a large variability: the estimated values could differ from up to a factor 30. Although the variability of the estimated plume rise is not so high, except for low wind speeds (Section 5), it produces an additional variability in the output results, especially since the plume might penetrate the inversion layer or not, depending on the parameterization.

The model evaluation against Prairie Grass and Kincaid experiments gave very satisfactory overall results, compared to other well-known Gaussian models. Similarity theory is the best standard-deviation scheme, which is not surprising since it relies on a more detailed description of the atmosphere. However, the Briggs formulae also compare well with other model results, not only with Prairie Grass experiment-on which they were adjusted-but also with Kincaid. Only the Doury parameterization appears to be unsatisfactory for the two experiments used here, since it tends to overestimate the plume spread, leading to underestimation of the concentrations. There is no striking difference in performance-for these cases-between "old generation" schemes, based on a discrete representation of the stability (Briggs, Pasquill-Guifford used in ISCST3) and the "new generation" estimations based on similarity theory, although the latter performs slightly better.

The results with Kincaid experiment point out the additional variability due to plume rise estimation. In particular, the "best" standard-deviation parameterization depends on the plume rise scheme used (similarity theory with Briggs-HPDM plume rise, and Briggs's standard deviations with Concawe and Holland formulae). All combinations perform well, except with the Doury formulae. The results also highlight the importance of some features such as a non-Gaussian shape for convective conditions, as well as plume lofting at the inversion layer. These missing features are most likely responsible for the relatively low correlations. A meteorological preprocessor that would recompute the inversion height could also help improving the results for Kincaid experiment.

With three parameterizations to compute the standard deviations, and the same number of plume rise schemes, we were able to show a large spread in the model outputs, representative of the variability of the Gaussian models. The main perspective is now to extend that work in order to include more parameterizations and to compare with more datasets, providing a full comparison between parameterizations in the same conditions. This study also suggests that ensemble modeling with Gaussian models could be interesting in the future. For that part we use formulae given in [START_REF] Hanna | Handbook on atmospheric diffusion[END_REF].

σ v = u * 12 -0.5 h L 1/3 unstable σ v = 1.3 u * exp -2 f z u * neutral σ v = max(1.3 u * 1 -z h , 0.2) stable
where h is the boundary layer height, L is the Monin-Obukhov length, f is the Coriolis parameter and z is the height where concentrations are computed.

The function F y is given by the general form

F y = 1 + 0.5 t τ L -1/2 (A1)
where τ L is the Lagrangian time scale. To compute this parameter, we use formulae given in [START_REF] Hanna | Applications in air pollution modeling[END_REF].

τ L = 0.15 h σv unstable τ L = 0.5 z σv / 1 + 15 f z u * neutral τ L = 0.07 h σv z h stable (A2)
where h is the boundary layer height, L is the Monin-Obukhov length, u * is the friction velocity, f is the Coriolis parameter and z is the height where concentrations are computed.

2) Vertical spread σ z

The vertical wind standard deviation is computed using formulae given in [START_REF] Venkatram | A semiempirical model to estimate dispersion of elevated releases in the stable boundary layer[END_REF] for stable and neutral cases and in [START_REF] Weil | Dispersion in the convective boundary layer[END_REF] for convective cases.

σ w = 0.6 w * unstable σ w = 1.3 u * 1 -z h 3/4 stable/neutral
where w * is the convective velocity, u * is the friction velocity, h the boundary layer height and L the Monin-Obukhov length.

The function F z is given by [START_REF] Irwin | Scheme for estimating dispersion parameters as a function of release height[END_REF] for stable and neutral case and by [START_REF] Weil | Dispersion in the convective boundary layer[END_REF] for convective cases.

• Unstable cases:

F z = 1 + 0.5 t τ Lz -1/2 (A3)
• Stable/Neutral cases:

F z =        1 + 0.9 ( t 50 ) -1 if z < 50 1 + 0.945 (0.1 t) 0.806 -1 if z ≥ 50
where h is the boundary layer height, L is the Monin-Obukhov length, z is the height where concentrations are computed. The following expression is advocated by [START_REF] Hanna | Handbook on atmospheric diffusion[END_REF] In [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF], an alternative set of formulae is proposed for sources of height greater than 100 m, since the other formulations based on similarity theory were calibrated on experiments closer to the ground and did not give satisfactory results in the case of the Kincaid experiment, for instance. These formulae are used for stable or neutral cases, where dispersion was found to be greater than given by the classical formulation.

• The wind vertical standard deviation is computed as:

σ w = 0.5 1.2 u 2 * + 0.35 w 2 * L ≥ 100 σ w = 1.3 u * 0 ≤ L < 100
• The wind horizontal standard deviation is computed as:

σ v = 0.7 3.6 u 2 * + 0.35 w 2 * L ≥ 100 σ v = max 1.5 3.6 u 2 * + 0.35 w 2 * , 0.5 0 ≤ L < 100
For the vertical standard deviation, the function F z does not change, but the vertical time scale is computed with a different set of formulae for stable and slightly unstable cases. In the following formulae, s p stands for the Briggs static stability parameter defined in Section 2.

• In the case of stable/neutral situations (L > 0.):

τ Lz = z/σ w if z ≤ L τ Lz = 0.27 √ s p if L ≤ 10 τ Lz = (z/σ w ) (L -10)/(z -10) + 0.27 √ s p (z -L)/(z -10) if 10 < L < z

  Research Laboratory, ENPC / EDF R&D, Université Paris Est Vivien Mallet INRIA

  Figure 3(a) shows the evolution of the horizontal standard deviations up to about three

(

  Corr). The formulae to compute these statistics are given in Appendix 8. It should be noted that these formulae give negative values of FB for model overprediction and positive values if the model tends to underpredict concentrations. A perfect model would have Corr and FAC2 equal to 1.0 and FB and NMSE equal to 0.0. As there is no such thing as a perfect model, it is useful to determine typical performance measures for acceptable models. Those criteria have been summarized in[START_REF] Chang | Air quality model performance evaluation[END_REF] and in[START_REF] Hanna | FLACS CFD air quality model performance evaluation with Kit Fox, MUST, Prairie Grass, and EMU observations[END_REF] based on the evaluation of many models with many field datasets. A model is deemed acceptable if the mean bias is within ±30% of the mean (-0.3 ≤ FB ≤ 0.3), if the random scatter is about a factor of two to three of the mean (NMSE ≤ 4), and if the fraction of predictions within a factor two of the observations is about 50 % (FAC2 ≥ 0.5).b. Experiment and modeling set-upThe Prairie Grass experiment has become a standard database on which parameterizations have been fitted, and has been used for many model-evaluation studies. The experiment took place in O'Neil, Nebraska, during summer 1956. The site was a flat terrain of short cut grass. A continuous plume of SO 2 was released, without plume rise, near the ground (at 0.46 m). Measurements were taken on five arcs at 50, 100, 200, 400 and 800 m from the source. A set of 43 of the 68 trials is usually used for model validations, including a wide range of stability conditions during day and night. In our simulations, the stability classes are taken from the database used in[START_REF] Hanna | FLACS CFD air quality model performance evaluation with Kit Fox, MUST, Prairie Grass, and EMU observations[END_REF] for FLACS evaluation, and the other experimental data are taken from the US EPA (Environmental Protection Agency) database (http://www.epa.gov/scram001/dispersion_prefrec.htm).c. Comparison with other Gaussian modelsMany Gaussian models and a few CFD models were evaluated with Prairie Grass experiment. Three of the most used Gaussian models were compared in CERC (2007): ISCST3 and AERMOD are recommended by the US EPA, and ADMS4 is developed in the UK and widely used in Europe. ISCST uses discrete categories to describe the atmospheric stability (Pasquill-Guifford parameterization for rural areas and Briggs's formulae for urban area), whereas AERMOD and ADMS use boundary layer depth and Monin-Obukhov length to describe the state of the atmospheric boundary layer. The CFD models FLACS (FLame ACceleration Simulator)[START_REF] Hanna | FLACS CFD air quality model performance evaluation with Kit Fox, MUST, Prairie Grass, and EMU observations[END_REF], and Mercure[START_REF] Demaël | Comparative evaluation of an eulerian CFD and gaussian plume models based on prairie grass dispersion experiment[END_REF]

  the Prairie Grass statistics for Polyphemus plume model with several parameterizations, along with the aforementioned Gaussian models. Results for ADMS4, AERMOD and ISCST3 come from(CERC 2007). The statistics for these three models are computed with a subset from the usual set of 43 experiments (used for Polyphemus evaluation), hence the differences in the mean observed centerline values. The statistics presented here are computed for the centerlines of the five arcs. Most values in Table1are within the acceptable range defined in Section a. AERMOD, ISCST3 and Polyphemus with the Briggs formulae show good results, as expected since the parameterizations they use for standard deviations are fitted on Prairie Grass experiment. On the contrary, the results for Polyphemus with Doury's formula and ADMS4 are hardly within the range of acceptable model performance for all indicators: ADMS4 underpredicts the mean concentrations by about 30%, and only 29% of the values predicted by the Doury model are within a factor two of the observations. Polyphemus with similarity theory shows very good results-there is almost no bias, and the correlation and NMSE are the best of all models presented here.

  Figure 4 shows the statistical indicators, described in Section a, computed on the centerline values for each arc separately, for the 43 experiments. The "acceptable range" as defined in Section a is delimited by the dashed lines. Most indicators show a decreasing performance at farther distances from the source, except the correlation. Whereas Briggs's and Doury's parameterizations underpredict the concentrations-this tendency increases with the distance (FB grows higher)-similarity theory overpredicts the concentrations, especially at the farthest arcs (400 m and 800 m).This can be related to the small standard deviation values given by this parameterization, pointed out in Section 3. At the first two arcs (50 m and 100 m), similarity theory gives the best results for all indicators. Finally, the Doury parameterization shows a poor performance in terms of FAC2, as well as for NMSE and FB, for all distances. In comparison to these results, the FLACS results reported in[START_REF] Hanna | FLACS CFD air quality model performance evaluation with Kit Fox, MUST, Prairie Grass, and EMU observations[END_REF] underestimated the concentra-tions at the shortest distances (up to 200 m), and slightly overpredicted the results at the farthest arcs.[START_REF] Demaël | Comparative evaluation of an eulerian CFD and gaussian plume models based on prairie grass dispersion experiment[END_REF] also found the underestimation tendency for Mercure CFD model at short distances. They attributed it to the Eulerian dispersion which corresponds to a plume behavior in σ ∝ √ t (with t the travel time) whereas the plume dispersion at short range should be in σ ∝ t. For C, D and E stability classes, the overall performance of the Mercure CFD model was shown to be between that of the Briggs and Doury parameterizations.

APPENDIX 7 .

 7 Formulae for standard deviations with similarity theory In all the formulae presented here, the stability is determined from the Monin-Obukhov length, as given in Table 8, which also shows the corresponding Pasquill class. a. Within the boundary layer 1) Transverse spread σ y First, we have to determine the horizontal standard deviation of wind velocity fluctuations σ v .

67Fig. 1 .Fig. 2 .Fig. 4 .

 124 Fig. 1. Representation of an elevated source: z s is the source height and z p is the effective source height, with plume rise. The modeled point of emission is at z p but the downwind travel is neglected during plume rise. The initial plume spread due to plume rise is σ zpr .

Fig. 5 .

 5 Fig. 5. Scatter plot of maximum arc concentrations for observations and simulations with various sigma parameterizations. Prairie Grass experiment, 43 trials. Unit is µg m -3 /(g s -1 ). The results are shown at five monitoring arcs (50 m,100 m, 200 m, 400 m and 800 m).

Table 1 .

 1 Statistics for several Gaussian models: comparison of maximum arc concentration for simulation and observation. Polyphemus results are presented for three formulae for standard deviations and 43 trials in Prairie Grass Experiment.

	(A4)

• In the case of slightly unstable situations (|L| ≤ 100):

τ Lz = 0.27 (z/σ w ) (0.55 -0.38 z/|L|) if z ≤ |L|

Table 2 .

 2 Proportion of cases with total or partial penetration in the inversion layer for all three plume rise parameterizations. Proportion with respect to 165 experiments.

	Model Mean	std	FB NMSE Corr FAC2
	Obs	39.89 39.22 0.00	0.00	1.00	1.00
	Briggs 33.01 28.26 0.19	1.46	0.21	0.42
	Doury 15.87 13.27 0.86	3.55	0.04	0.28
	Sim.th. 61.97 47.52 -0.43	1.38	0.24	0.38

Table 4 .

 4 Statistics for Polyphemus Gaussian models: comparison of centerline values C/Q for simulations and observations-Kincaid data of quality 2 and 3. There are 165 experiments. Unit is ng m -3 /(g s -1 ). Plume rise is computed with Holland formulae.

	Model Mean	std	FB NMSE Corr FAC2
	Obs	39.89 39.22 0.00	0.00	1.00	1.00
	Briggs 32.88 32.75 0.19	1.49	0.27	0.48
	Doury 16.87 17.74 0.81	3.29	0.12	0.31
	Sim.th. 61.82 53.79 -0.43	1.57	0.24	0.43

Table 5 .

 5 Statistics for Polyphemus Gaussian models: comparison of centerline values C/Q for simulations and observations-Kincaid data of quality 2 and 3. There are 165 experiments. Unit is ng m -3 /(g s -1 ). Plume rise is computed with Concawe formulae.

	Model	Mean	std	FB NMSE Corr FAC2
	Obs	41.0	39.3	0.00	0.00	1.00	1.00
	ADMS4	40.4	31.1 -0.05	0.80	0.50	0.58
	AERMOD	20.3	24.1	0.68	2.30	0.35	0.33
	ISCST3	23.1	53.3	0.56	3.83	0.26	0.26
	Obs	39.89 39.22 0.00	0.00	1.00	1.00
	Briggs/Concawe 32.83 32.80 0.19	1.50	0.27	0.47
	Doury/Concawe 16.75 17.81 0.82	3.32	0.12	0.30
	Sim.th./Concawe 61.70 53.91 -0.43	1.58	0.25	0.43
	Sim.th/HPDM	41.99 34.76 -0.05	1.33	0.19	0.44

Table 6 .

 6 Statistics for several Gaussian models: comparison of centerline values C/Q for simulations and observations-Kincaid data of quality 2 and 3. Unit is ng m -3 /(g s -1 ).

	Model	Mean	std	FB NMSE Corr FAC2
	Obs	54.34 40.25 0.00	0.00	1.00	1.00
	ADMS4	48.5	31.5	0.11	0.60	0.45	0.68
	AERMOD	21.8	21.8	0.86	2.07	0.40	0.29
	ISCST3	30.0	60.0	0.57	2.80	0.26	0.28
	HPDM	44.84 38.55 0.19	0.75	0.44	0.56
	SAFE AIR	59.73 52.54 -0.09	1.10	0.20	0.50
	OML	47.45 45.48 0.13	1.24	0.15	0.55
	NAME	38.7	47.2	0.33	1.45	0.27	0.56
	Obs	53.69 40.72 0.00	0.00	1.00	1.00
	Briggs/Concawe 37.95 33.06 0.34	1.20	0.21	0.58
	Doury/Concawe 17.36 14.21 1.02	3.19	0.18	0.30
	Sim.th./Concawe 68.91 51.25 -0.25	1.03	0.17	0.58
	Sim.th/HPDM	47.92 33.93 0.11	1.05	0.05	0.59

Table 7 .

 7 Statistics for several Gaussian models: comparison of centerline values C/Q for simulation and observation-Kincaid data of quality 3. Unit is ng m -3 /(g s -1 ).

	L: Monin-Obukhov length (m) Stability Pasquill class
	-100 ≤ L < 0 -10 5 ≤ L < -100 |L| > 10 5 10 ≤ L < 10 5 0 ≤ L < 10	unstable unstable neutral stable stable	B C D E F

Table 8 .

 8 Matching between Monin-Obukhov length and stability, from[START_REF] Seinfeld | Atmospheric chemistry and physics: from air pollution to climate change[END_REF].

This could lead to significant differences, since this height is of great importance to know whether the plume touches the ground or is trapped above the inversion layer, as explained in Section 2. Simulations were carried out with three parameterizations for standard deviations (similarity-theory, Briggs's and Doury's) and three parameterizations for plume rise described in Section 2. The chimney diameter as well as the initial plume spread due to plume rise were taken into account in the total standard deviations (Equation 2). When similarity theory is used, the alternative HPDM formulae for elevated release detailed in Appendix 7 are used. Compared to the usual similarity-theory parameterization, it slightly increases the overestimation but it also improves the correlation coefficient, especially with Briggs-HPDM plume rise. The model evaluation consists, first, in comparing the results with the different parameterizations, with an emphasis on the plume rise description. In a second part, we compare the results with those of other models.

b. Model evaluation

1) Plume rise

A first analysis of the plume rise values is done in this section, so as to compare the plume rise formulations and to give an insight on the additional variability due to plume rise. Figure 6 shows the plume rise values computed with the three parameterizations and 165 experiments, against the wind speed. As expected, the plume rise decreases with higher wind speeds, smoothing discrepancies between the different schemes. It was pointed out, for example in [START_REF] Canepa | About the plume rise description in the dispersion code SAFE AIR[END_REF], that the Briggs scheme has a tendency to overestimate plume rise at low wind speeds. Indeed, the plume rise given by the Briggs-HPDM scheme

c. Above the boundary layer 1) HPDM formulae for wind standard deviations All the formulae given by similarity theory are valid only within the boundary layer. There are not many studies about the dispersion above the boundary layer, so the formulae used for that case are quite simple. They come from [START_REF] Hanna | Hybrid plume dispersion model (HPDM) development and evaluation[END_REF].

They are applied for all stabilities. They ensure that the dispersion parameters at least are always strictly positive values (no division by zero or square root of a negative value...). However, since in most case studies the observed concentrations are inside the boundary layer, and there is no mixing through the inversion layer, these are not the most used formulae in practice.

2) Gillani formula for vertical standard deviation

An alternative way to compute σ z above the boundary layer is given by [START_REF] Gillani | Science algorithms of the EPA models-3 Community Multi-scale Air Quality (CMAQ) modeling system[END_REF]:

assuming that t is the time travel of the plume, as usual, and the initial time is equal to 0. One can note that this formula can be applied above the boundary layer for all parameterizations.

The initial plume spread σ zpr due to plume rise is computed differently from the formula given in Section 2 since the initial spread is often overestimated that way. An alternative way proposed by [START_REF] Gillani | Science algorithms of the EPA models-3 Community Multi-scale Air Quality (CMAQ) modeling system[END_REF] is:

where dT /dz is the vertical ambient temperature gradient.

Statistical indicators

where C p and C 0 are the predicted and measured maximum arcwise concentrations respectively, α is the average of value α over the dataset (that is, for all experiments and arcs), and σ α is the standard deviation over the dataset. ). Plume rise is computed with Briggs-HPDM formulae.

List of Tables

List of Figures


Luciano Campi
email: campi@ceremade.dauphine.fr

Sara Biagini

José M Corcuera

A note on market completeness with American put options

Keywords: market completeness, American put, Tanaka's formula, European call, marginals MS Classification (2000): 91B28, 60G40, 60G44, 60G48

We consider a non necessarily complete financial market with one bond and one risky asset, whose price process is modelled by a suitably integrable, strictly positive, càdlàg process S over [0, T ]. Every option price is defined as the conditional expectation under a given equivalent (true) martingale measure P, the same for all options. We show that every positive contingent claim on S can be approximately replicated (in L 2 -sense) by investing dynamically in the underlying and statically in all American put options (of every strike price k and with the same maturity T ). We also provide a counter-example to static hedging with European call options of all strike prices and all maturities t ≤ T .

Introduction

The aim of this paper is to investigate the issue of hedging in a market where agents are allowed to invest continuously on time in a risky asset and statically in American put options of all strike prices and with a fixed maturity T . Such additional investment opportunities are of particular interest in incomplete markets, in which case the payoffs cannot, in general, be replicated by a trading strategy in the underlying.

Since in real financial markets many types of options are becoming more and more liquid, it is very natural to reformulate hedging and optimal investment problems incorporating those larger trading opportunities. Indeed, in the recent years many papers have treated problems like absence of arbitrage, hedging, optimal portfolio choice in a financial market where investors are allowed to trade in the underlying assets as well as to assume static positions in some class of derivatives. Here, we recall only few of them: Campi [START_REF] Campi | Arbitrage and completeness in financial markets with given N -dimensional distributions[END_REF] for no-arbitrage and completeness issues, the papers by Ilhan et al. [START_REF] Ilhan | Optimal static-dynamic hedges for barrier options[END_REF][START_REF] Ilhan | Optimal investment with derivative securities[END_REF] and Carr et al. [START_REF] Carr | Optimal investment in derivative securities[END_REF] for optimal investment problems, and the more recent papers by Schweizer and Wissel [START_REF] Schweizer | Term structures of implied volatilities: Absence of arbitrage and existence results[END_REF][START_REF] Schweizer | Arbitrage-free market models for option prices: The multi-strike case[END_REF] and by Jacod and Protter [START_REF] Jacod | Risk Neutral Compatibility with Option Prices[END_REF] where an HJM approach for European call options is developed.

Our paper is much closer in spirit to that part of literature initiated by Ross [START_REF] Ross | Options and efficiency[END_REF], where it has been shown that in a single period model with a finitely many states and n stocks, simple options (i.e. written on only one stock) can be replicated by trading in European call (or put) options on that stock. Ross [START_REF] Ross | Options and efficiency[END_REF] established also the existence of a portfolio in all n stocks such that call (or put) options written on that portfolio span the set of all options, simple as well as complex, i.e. written on all stocks at the time. This result has been generalized in an infinite state space setting by Green and Jarrow [START_REF] Green | Spanning and completeness in markets with contingent claims[END_REF] (see also Nachman [16]). For a similar result in a dynamic discrete-time setting with a general state space, we quote only a paper by Rogge [START_REF] Rogge | Call completeness implies completeness in the n-period model of a financial market[END_REF], where the market completeness is characterized in terms of a so-called "call-completeness": there exists a strike price k > 0 such that if a claim is attainable then also a call written on it with strike k is.

A more explicit result has been established in Carr and Madan [START_REF] Carr | Optimal positioning in derivative securities[END_REF], where it is shown that a European option with pay-off f (S T ) (f sufficiently regular) can be hedged by a unique initial position of f (S 0 ) -f ′ (S 0 )S 0 unit discount bonds, f ′ (S 0 ) stock shares, and f ′′ (k)dk out-of-the-money call and put options of all strikes. Even more importantly, this result is model-free, so that it holds also in a continuous-time framework. For static hedging of exotic options in a diffusion setting, we refer to Carr et al. [START_REF] Carr | Static Hedging of Exotic Options[END_REF]. Moreover, in a discontinuous setting, Corcuera et al. [START_REF] Corcuera | Completion of a Lévy market by power-jump assets[END_REF] proved that a Lévy market can be completed by trading statically in power-jump assets, which are somewhat related to contracts on realized variance. In [START_REF] Corcuera | Dynamic complex hedging in additive markets[END_REF], a link between power-jump assets and European call options is established, allowing Corcuera and Guerra to prove that a Lévy market can be completed by trading in a continuum of European call options along strikes. Finally, in the paper [START_REF] Davis | Market completion using options[END_REF], Davis and Obloj provide necessary and sufficient conditions for a market model to be complete by trading in given finite set of European-type derivatives driven by finitely many factors modelled as diffusion processes.

To sum up, all these papers deal with the problem of finding a good class of derivatives capable to span all contingent claims written on the same underlying. By "good class of derivatives" we mean a class of options sufficiently liquid in real financial markets in order to reproduce with a certain precision the hedging strategies suggested in these papers.

The present paper identifies American put options of all strike prices and with the same maturity as a good class of derivatives allowing an investor to replicate approximately (in L 2 -sense) every positive contingent claim written on the same underlying in a very general frictionless financial market with one risk-free asset and one stock. Our main result, stated in Theorem 6, can be viewed as a generalization to a continuous-time setting of previously quoted Ross's result in [START_REF] Ross | Options and efficiency[END_REF] -notice that in Ross's single period setting there is no need to distinguish between European and path-dependent options like American puts. Moreover, our result is model-free. Indeed, we will make only the very mild assumption that the price process is a square-integrable and possibly discontinuous martingale with no jump at maturity. The price we have to pay is that the agent has to trade in infinitely many securities, in contrast with, e.g., Davis and Obloj [START_REF] Davis | Market completion using options[END_REF] where, since agents can trade only in finitely many assets, the authors need to assume more regularity as, e.g., analiticity of the coefficients together with a non-singularity condition, to get the market completeness.

The present paper is structured as follows: In Section 2, we give a short description of the model, Section 3 contains the main result on hedging with American put options. Finally, in Section 4, we exhibit a two-periods market model where static investments in all European call (equivalently, put options) of all strike prices and all maturities are not enough for replicating all contingent claims, so justifying the use of their American counterparts.

The model

We consider a financial market composed by a riskless asset and a stock. More precisely, let (Ω, F, P) be a probability space equipped with a filtration (F t ) t∈[0,T ] , where T > 0 is a finite horizon, F = F T and F 0 is trivial. (F t ) t∈[0,T ] is the natural filtration -right-continuous and P-saturated -generated by a strictly positive, càdlàg process S = (S t ) t∈[0,T ] modelling the price of the risky asset. We assume without loss of generality that the price at time t of the riskless asset satisfies S 0 t = e rt where r > 0 is the spot interest rate (see Remark 7 for a straightforwad generalization). In the sequel E[•] will denote expectation with respect to P and, for any process X, we set X := X/S 0 and ∆X t = X t -X t-for any t ∈]0, T ]. Finally, H 2 (P) will denote the space of all martingales bounded in L 2 (P) := L 2 (P, F T ).

Our first assumption concerns the behaviour of the risky asset process S.

Assumption 1 S = S/S 0 belongs to H 2 (P) and it does not jump at T , i.e. ∆S T = 0 a.s..

Since by assumption (F t ) t∈[0,T ] is the natural filtration of S, the equality ∆S T = 0 a.s. is equivalent to F T = F T -, i.e. the underlying filtration is left-continuous at T . We notice that this property is satisfied, e.g., by all diffusion models and also by all exponential Lévy models. Indeed, it is well-known that the natural filtration of any Lévy process is quasi-left continuous, i.e. it does not jump at predictable stopping times (see, e.g., Protter [START_REF] Protter | Stochastic integration and differential equations[END_REF], p. 150, Exercises 8 and 9, and p. 191 for details).

In this setting, a contingent claim pay-off on S with maturity T is naturally modelled by a random variable f ∈ L 2 (P). This model is not necessarily complete, so there may exist infinitely many equivalent (local) martingale measures different from P and so, equivalently, not every contingent claim can be hedged by trading only in the underlying, so that allowing an investor to trade also in some class of options can enlarge considerably his hedging opportunities.

The following assumption can be viewed as a kind of no-arbitrage consistency among prices of different options written on the same underlying S.

Assumption 2 In this financial market every option price comes from the same equivalent martingale measure P.

Observe that Assumption 2 has a very natural financial motivation. Indeed, let us forget about integrability issues for a while and consider two contingent claims with pay-offs f and f ′ and corresponding risk-neutral prices f t and f ′ t at time t. Assume that they come from two equivalent martingale measures Q and Q ′ respectively, i.e.

f t = E Q [f |F t ] and f ′ t = E Q ′ [f ′ |F t ].
If we want the enlarged market (S, f, f ′ ) to be arbitrage-free, there must be an equivalent martingale measure Q ′′ for (S, f, f ′ ), so that one has

f t = E Q [f |F t ] = E Q ′′ [f |F t ]
and the same for f ′ . This means that one can assume without loss of generality that the risk-neutral prices of f and f ′ come from the same measure Q ′′ .

Hedging with American put options

In this section, we will prove that additional static investments in American put options of every strike price and with the same maturity T allow an investor to hedge any contingent claim written on the underlying S. We make the following model assumption:

Assumption 3 American put options with any maturity T and any strike price k > 0 are available for trading and all are issued at t = 0.

Let us denote T the set of all [0, T ]-valued stopping times and τ k the smallest optimal exercise time corresponding to a given strike price k. Such stopping times do exist (see Remark 4 for details) so that we have

P 0 (k, τ k ) := sup τ ∈T E[e -rτ (k -S τ ) + ] = E[e -rb τ k (k -S b τ k ) + ], (1)
which is the price that an agent must pay for buying at time 0 an American put option with strike k and maturity T . More generally, we set

Pt (k, τ ) := E[e -rτ (k -S τ ) + | F t ],
for any stopping time τ ∈ T . Notice that under our assumptions each price process P (k, τ ) is a martingale in H 2 (P).

Remark 4 Observe that τ k exists for all strike prices k > 0. Indeed, the underlying discounted price process S is a strictly positive martingale in H 2 (P) over a finite time interval [0, T ], so that it is of class (D) and, obviously, constant in expectations. Thus, Théorème 2.43 in El Karoui St. Flour's Lecture Notes [START_REF] El Karoui | Les aspects probabilistes du contrôle stochastique[END_REF] applies, so giving existence of τ k for all k > 0.

Now, we consider the set R a of all discounted contingent claims which can be approximately replicated by investing dynamically in the underlying S and statically in finitely many American put options as follows: the American puts that the agent buys at time t = 0 can be exercised at any stopping time τ ∈ T , while the American puts that he sells will be exercised at their corresponding optimal times by their buyers.

Mathematically speaking, R a is the set of all F T -measurable random variables of the form

x + T 0 θ t d St + n i=1 α i ( PT (k i , τ i ) -P 0 (k i , τ ki )) - m j=1 β n+j ( PT (k n+j , τ kn+j ) -P 0 (k n+j , τ kn+j )) (2)
where

• x ∈ R is the initial endowment of the agent,
• θ is a real-valued S-integrable predictable process such that θd S is a martingale in H 2 (P) modelling the dynamic investment strategy in the stock S,

• n ≥ 0 is the number of American puts that the agent buys at time 0, while m ≥ 0 is the number of American puts sold at time 0,

• each weight α i ≥ 0 is a nonnegative real number representing the number of American puts with strike k i (i = 1, . . . , n) bought by the agent at time 0 paying the price P 0 (k i , τ ki ), while β n+j ≥ 0 (j = 1, . . . , m) represents the number of American puts with strike k n+j sold by the agent at time 0 receiving the price P 0 (k n+j , τ kn+j ).

We adopt the convention that any summation over an empty set of indexes is equal to zero. Notice that R a is not a linear subspace, it is nonetheless a convex cone. We will denote R * a its (positive) dual, i.e. R * a := {f ∈ L 2 (P) : E[f g] ≥ 0, ∀g ∈ R a }, and by R * * a its bidual, i.e. R * a := (R * a ) * . We recall that the bidual C * * of any convex cone C in a vector topological space coincides with the closure of C, and that for any pair of convex cones C and C ′ such that C ⊂ C ′ one has C * ⊃ (C ′ ) * (see, e.g, [START_REF] Zȃlinescu | Convex analysis in general vector spaces[END_REF], Chapter 1).

Remark 5 Notice that in the first summation appearing in (2), denoting the final gain coming from a long position taken at time t = 0 in American puts with strikes k i , the puts are not necessarily exercised at their optimal exercise times τ ki . The agent, willing to hedge against the risk of a given final pay-off f , can in principle exercise his puts at any stopping time between today and the maturity T . The question if hedging purposes may lead an agent to exercise such options at sub-obtimal times remains open.

Our main result is that R a ∩ L2 + (P) is dense in L 2 + (P) which denotes the set of all positive random variables in L 2 (P). In financial terms, it means that every positive contingent claim can be (approximately) replicated by a mixed investment as in (2): dynamic in the underlying and static in American put options.

Theorem 6 Under our assumptions, the closure of R a contains L 2 + (P).

Proof. We recall once more that the closure of a convex cone C in a topological vector space equals its bidual C * * , so that all we need to prove is that the bidual R * * a contains L 2 + (P). In order to do that, it suffices to show that L 2 + (P) contains the dual cone R * a , i.e. that any random variable

N T ∈ L 2 (P) which is positive over R a , i.e. E[N T f ] ≥ 0 for all f ∈ R a , is a.s. positive itself. Denote N the martingale in H 2 (P) associated to N T , i.e. N is the càdlàg version of the martingale E[N T |F t ], for t ∈ [0, T ].
First observe that, since R a contains R, we have E[N T ] = N 0 = 0. From the fact that E[N T f ] ≥ 0 for any f equal to the static hedging part appearing in (2), we can deduce in particular that

α 1 E[N T PT (k 1 , τ 1 )] -α 2 E[N T PT (k 2 , τ k2 )] ≥ 0
for all nonnegative real numbers α 1 , α 2 , all strikes k 1 , k 2 , and all stopping times τ 1 ∈ T . Taking α 2 = 0, we get that

0 ≤ E[N T PT (k 1 , τ 1 )] = E[N T Pτ1 (k 1 , τ 1 )] = E[N τ1 e -rτ1 (k 1 -S τ1 ) + ] = E[ Ñτ1 (k 1 -S τ1 ) + ],
for every strike price k 1 and every stopping time τ 1 ∈ T . Doob's optional sampling theorem implies that the process Ñ (k -S) + is a P-submartingale for every k > 0. The integration by parts formula gives

Ñt (k -S t ) + = Ñt (S t -k) -= t 0 (S u--k) -d Ñu + t 0 Ñu-d(S -k) - u + [ Ñ , (S -k) -] t
for all t ∈ [0, T ]. By Tanaka's formula for discontinuous semimartingales (e.g. Protter [START_REF] Protter | Stochastic integration and differential equations[END_REF], Theorem 68, p. 216) and since dS u = rS u du + e ru d Su , we have that

t 0 Ñu-d(S -k) - u = - t 0 Ñu-1 {Su-≤k} (rS u du + e ru d Su ) (3)
+ 0<u≤t Ñu-1 {Su->k} (S u -k) -+ 1 {Su-≤k} (S u -k) + + 1 2 t 0 Ñu-dL k u (S), (4)
for every instant t ∈ [0, T ] and strike price k > 0.

The fact that N T belongs to the dual of R a implies also that N T is weakly orthogonal to S and so strongly orthogonal to S as well,1 i.e. N S = Ñ S belongs to H 1 (P), the space of all martingales bounded in L 1 (P) (use, e.g., Lemma 2, Section IV, in [START_REF] Protter | Stochastic integration and differential equations[END_REF]). Moreover, being

d[ Ñ , (S -k) -] t = e -rt d[N, (S -k) -] t and L k (S) a continuous increasing process, one has [ Ñ , (S -k) -] t = 0<u≤t ∆ Ñu ∆(S -k) - u ,
which is a pure jump process. Notice that Ñ (k -S) + is a P-submartingale with Doob-Meyer decomposition Ñ (k -S) + = M + B, where M is a local martingale and B is a predictable increasing process. The local martingale part M include certainly the following term

t 0 (S u--k) -d Ñu - t 0
Ñu-1 {Su-≤k} e ru d Su and, since any càdlàg local martingale with finite variation must be purely discontinuous (see, e.g., [15, Lemma 4.14 b)]), M cannot contain the finite variation terms with continuous paths appearing in ( 3) and ( 4). Thus, those terms have to belong to the increasing part B, which could in principle contains some additional terms of pure jump type. An important consequence of it is that the process 1 2

t 0 Ñu-dL k u (S) - t 0 rk Ñu 1 {Su≤k} du, t ∈ [0, T ],
must be increasing, i.e.

t s rk Ñu 1 {Su≤k} du ≤ 1 2 t s Ñu-dL k u (S), (5)
for all s, t ∈ [0, T ] with s ≤ t and all k > 0. Using a standard class monotone argument, the inequality in ( 5) can be generalized as follows

A rk Ñu 1 {Su≤k} du ≤ 1 2 A Ñu-dL k u (S), (6)
for all Borel set A and all k > 0. First, observe that in the two integrals in (6) the same function u → Ñuis integrated with respect to two dt-a.e. mutually singular2 measures (1/2)dL k u (S) and rk1 {Su≤k} du, which implies that such an inequality is verified only if N ≤ 0 dP ⊗ dt-a.e. on the set {(ω, t) : S t (ω) ≤ k}, for all k > 0. As a consequence, one has N ≤ 0 dP ⊗ dt-a.e. on Ω × [0, T ] and, since N is càdlàg, one has also that N T -≤ 0 a.s.. Finally, notice that Assumption 1 consequence F T = F T -implies that no martingale can jump at T (see Protter [START_REF] Protter | Stochastic integration and differential equations[END_REF], p. 191, for details), so that one has also that N T ≤ 0 a.s.. To end the proof, it suffices to recall that N is martingale with N 0 = 0, so that E[N T ] = 0 and, N T being a.s. negative, it must be a.s. zero.

Remark 7 A careful inspection of our proof reveals that Theorem 6 holds true even if the spot interest rate r is not necessarily constant but a positive and bounded deterministic function of time, more precisely: S 0 t = exp( t 0 r(u)du) where r(u) is measurable positive function defined on [0, T ] and such that S 0 T is bounded above by some constant.

Remark 8 Note that if one considers contingent claims depending on some randomness source different from S, then the previous completeness result breaks down. Indeed, just consider a model whose price processes are identically equal to one, i.e. S 0 ≡ S ≡ 1, so that the natural filtration of S is trivial and the collection of all American put option pay-offs {(k -1) + : k > 0} coincides with [0, ∞). Thus, R a = R. Then, consider a sufficiently large filtration (F t ) t∈[0,T ] such that F T contains at least one non-degenerate square-integrable positive random variable f . It is now clear that, even if we allow the strategies to be F-adapted, in such a market it is not possible to hedge f as in Theorem 6.

A counter-example to hedging with European call options

In this section, we will exhibit a financial market in discrete time with a finite horizon T ∈ N and defined on a finite probability space, where it is not possible to hedge all contingent claims by trading dynamically in a given underlying and statically in all European call options of every strike price k > 0 and every maturity before T . Let (Ω, F, P) be a finite probability space supporting a martingale S = (S t ) T t=0 modelling the price evolution of a stock. This space is assumed to be equipped with the filtration (F t ) T t=0 naturally generated by S and for which F = F T . As usual, we denote S 0 the price process for a riskless asset and assume S 0 ≡ 1. For a given discrete-time process X, we set ∆X t := X t -X t-1 , t = 1, . . . , T .

Consider the linear space R e spanned by all random variables f of the form

f = x + T t=1 θ t ∆S t + n i=1 α i (C T (T i , k i ) -C 0 (T i , k i )), (7)
where x ∈ R is an initial endowment, θ is any predictable process modelling the dynamic strategy in S, α = (α 1 , . . . , α n ) ∈ R n is any static strategy in n ≥ 1 European call options with maturities T i ≤ T and strike prices k i for 1 ≤ i ≤ n, whose no-arbitrage prices are denoted by

C t (T i , k i ) := E[(S Ti -k i ) + |F t ]
. n ≥ 1 is an arbitrary positive integer.

Our aim is to construct a process S such that R e is not dense in the set of all positive F Tmeasurable random variables L 0 + , equipped with the usual scalar product (f, g) = E[f g]. To do so, we use the following consequence of Theorem 3 in Campi [START_REF] Campi | Arbitrage and completeness in financial markets with given N -dimensional distributions[END_REF]. We provide its short proof for reader's convenience.

Lemma 9 Assume that R e ∩ L 0 + is dense in L 0 + . Then the set of all P-equivalent martingale measures Q under which S has the same marginals as under P reduces to a singleton.

Proof. Let Q be an equivalent martingale measure under which S has same marginals as under P. In this case, for all postive random variable f as in [START_REF] Corcuera | Completion of a Lévy market by power-jump assets[END_REF] we have E[f ] = x = E Q [f ] and, since the family of those random variables is assumed to be dense in L 0 + , we can conclude that Q = P on F T . In the light of that result, it suffices now to exhibit a process S admitting two different equivalent martingale measures P and Q under which S has the same marginals. Here it is: T = 2, S 0 = 3/2, the marginals at time t = 1 are given by P[S 1 = 1] = Q[S 1 = 1] = 1/2 and P[S 1 = 2] = Q[S 1 = 2] = 1/2, and S 2 takes the values 0, 1, 2, 3 each one with probability 1/4 under both P and Q. To complete the description of P and Q, we only need to assign the transition probabilities between t = 1 and t = 2. This can be done in many ways to get P, Q ∈ M and nonetheless keep them different. For instance, set p ij := P[S 2 = j|S 1 = i] and q ij := Q[S 2 = j|S 1 = i] for i, j ∈ {0, 1, 2, 3}, and consider p 23 = p 10 = 0.4, p 22 = p 11 = 0.3, p 21 = p 12 = 0.2 and p 20 = p 13 = 0.1 for the measure P, and q 23 = q 10 = 0.39, q 22 = q 11 = 0.33, q 21 = q 12 = 0.17 and q 20 = q 13 = 0.11 for the measure Q. It can be easily verified that this example is exactly what we were looking for.

Remark 10 As it is formulated, this example does not satisfy the assumption S T = S T -. Indeed, S T -= S 2-= S 1 = S 2 . Nonetheless, it can be easily embedded in the framework of the previous section, where the price process S is assumed to be left-continuous at T as in Assumption 1, by simply adding a date T + 1 and setting S T +1 = S T .

Remark 11 By the call-put parity, trading in all European call options as in ( 7) is equivalent to trading in S and in all European put options. As a consequence, our example also shows that it is in general not possible to replicate each square-integrable positive contingent claim by trading dynamically in the underlying and statically in all European put options.

We recall that, in this H

setting, weak orthogonality between two martingales M and N in H 2 (P) is equivalent to strong orthogonality, i.e. M N ∈ H 1 (P). See, e.g., Lemma 2, Section IV, in Protter's book[START_REF] Protter | Stochastic integration and differential equations[END_REF].

Indeed, the support of dL k u (S) is {u : Su = S u-= k} (see, e.g., Protter's book[START_REF] Protter | Stochastic integration and differential equations[END_REF], Theorem 69, p. 217) while that of 1 {Su≤k} du is {u : Su ≤ k}. Thus, their intersection is contained in {u : Su = k} which is at most countable and so it has zero Lebesgue measure.


Etienne Côme
email: etienne.come@ifsttar.fr

Marie Cottrell
email: marie.cottrell@univ-paris1.fr

Michel Verleysen
email: michel.verleysen@uclouvain.be

Jérôme Lacaille
email: jerome.lacaille@snecma.fr

Aircraft engine fleet monitoring using Self-Organizing Maps and Edit Distance

Aircraft engines are designed to be used during several tens of years. Ensuring a proper operation of engines over their lifetime is therefore an important and difficult task. The maintenance can be improved if efficient procedures for the understanding of data flows produced by sensors for monitoring purposes are implemented. This paper details such a procedure aiming at visualizing in a meaningful way successive data measured on aircraft engines and finding for every possible request sequence of data measurement similar behaviour already observed in the past which may help to anticipate failures. The core of the procedure is based on Self-Organizing Maps (SOM) which are used to visualize the evolution of the data measured on the engines. Rough measurements can not be directly used as inputs, because they are influenced by external conditions. A preprocessing procedure is set up to extract meaningful information and remove uninteresting variations due to change of environmental conditions. The proposed procedure contains four main modules to tackle these difficulties: environmental conditions normalization (ECN), change detection and adaptive signal modeling (CD), visualization with Self-Organizing Maps (SOM) and finally minimal Edit Distance search (SEARCH). The architecture of the procedure and of its modules is described in this paper and results on real data are also supplied.

Introduction

During the flights, some on-board sensors measure many parameters related to the behavior (and therefore the health) of aircraft engines. These parameters are recorded and used at short and long terms for immediate action and alarm generation, respectively. In this work, we are interested in the long-term monitoring of aircraft engines and we want to use these measurements to detect any deviations from a "normal" behavior, to anticipate possible faults and to facilitate the maintenance of aircraft engines. This work presents a tool that can help experts, in addition to their traditional tools based on quantitative inspection of some relevant variables, to easily visualize the evolution of the engine health. This evolution will be characterized by a trajectory on a two-dimensional Self-Organizing Map. Abnormal aging and fault will result in deviations with respect to normal conditions. The choice of Self-Organizing Maps is motivated by several points:

-SOMs are useful tools for visualizing high-dimensional data onto a lowdimensional grid; -SOMs have already been applied with success for fault detection and prediction in plants and machines (see [START_REF] Svensson | Self-organizing maps for automatic fault detection in a vehicle cooling system[END_REF] for example).

This article follows another WSOM paper [START_REF] Cottrell | Fault prediction in aircraft engines using self-organizing maps[END_REF] but contains necessary material (and possibly redundant) to be self-contained. It is organized as follows : first, in Section 2, the data and the notations used throughout the paper are presented. The methodology and the global architecture of the proposed procedure are described in Section 3. Each step is defined and results on real data are given in Section 4.

Data

Measurements are collected on a set of I engines. On each engine i ∈ {1, . . . , I}, n i measurements are performed successively flight after flight; there is thus no guarantee that the time intervals between two measures are approximately equal. Each observation is denoted by Z ij , where i ∈ {1, . . . , I} is the engine number and j ∈ {1, . . . , n i } is the flight number.

Each vector Z ij contains two kinds of variables: those which are strictly related to the behavior of the engine (fuel consumption, static pressure, ...), and those which are related to the environment (temperature, altitude, ...). Let the p engine variables be denoted by Y 1 ij , . . . , Y p ij and the q environmental variables by X 1 ij , . . . , X q ij . Each observation is therefore a (p + q)-vector Z ij , where

Z ij = [Y ij , X ij ] = Y 1 ij , . . . , Y p ij , X 1 ij , . . . , X q ij .
The variables at disposal are listed in Table 1. There are p = 5 engine variables and q = 15 environmental variables. The dataset contains measurements for approximately one year of flights and I = 91 engines, that leads to a global dataset with 91 i=1 n i = 59407 (p + q)dimensional observations.

Methodology

The goal is to build the trajectories of all the engines, that is to project the successive observations of each engine on a Self-Organizing Map, in order to follow the evolution and to eventually detect some "abnormal" deviation. It is not valuable to use the rough engine measurements: they are inappropriate for direct analysis by Self-Organizing Maps, because they are strongly dependent on environment conditions and also on the characteristics of the engine (its past, its age, ...). The first idea is to use a linear regression for each engine variable: the environmental variables (real-valued variables) and the number of the engine (categorical variable) are the predictors and the residuals of these regressions can be used as standardized variables (see [START_REF] Cottrell | Fault prediction in aircraft engines using self-organizing maps[END_REF] for details). For each engine variable r = 1, . . . , p, the regression model can be written as:

Y r ij = µ r + α r i + λ r 1 X 1 ij + . . . + λ r q X q ij + ǫ r ij (1)
where α r i is the engine effect on the r th variable, λ r 1 , . . . , λ r q are the regression coefficients for the r th variable, µ r is the intercept and the error term ǫ r ij is the residual. Figure 1 presents for example the rough measurements of the corespeed feature as a function of time (for engine 6) and the residuals computed by model [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]. The rough measurements seem almost time-independent on this figure, whereas the residuals exhibit an abrupt change which is linked to a specific event in the 1).

life of this engine. This simple model is therefore sufficient to bring to light interesting aspects of the evolution of this engine. However, the signals may contain ruptures, making the use of a single regression model hazardous. The main idea of this work is to replace model (1) by a new procedure which deals with the temporal behavior of the signals. The goal is therefore to detect the ruptures and to use different models after each rupture. This new procedure is composed of two modules. The first module (Environmental Conditions Normalization, ECN) aims at removing the effects of the environmental variables to provide standardized variables, independent of the flight conditions. It is described in section 4.1. The second module uses an on-line change detection algorithm to find the above mentioned abrupt changes, and introduces a piecewise regression model. The detection of the change points is done in a multi-dimensional setting taking as input all the normalized engine variables supplied by the ECN module. The Change Detection (CD) module is presented in Section 4.2. As a result of these first two steps, the "cleaned" database can be used as input to a Self-Organizing Map with a "proper" distance for trajectories visualization. The third module (SOM) provides the "map" on which the trajectories will be drawn. Finally, engine trajectories on the map are gathered in a trajectory database which can be accessed through a SEARCH module, which use a dedicated Edit Distance to find similar trajectories. This four-steps procedure is summarized in Figure 2.

4 Description of the four modules

Environmental Conditions Normalization -ECN

The first module aims at removing the effects of the environmental variables.

For that purpose, one regression model has to be fitted for each of the p engine variables. As the relationship between environmental and engine variables is complex and definitively not linear, the environmental variables can be supplemented by some non-linear transformations of the latter, increasing the number of explanatory variables. Interactions (all the possible products between two environmental variables), squares, cubes and fourth powers of the non binary environmental variables are considered. The number q of predictors in the model is therefore a priori equal to (11 + 4) * (11 + 4 -1)/2 = 105 for the interactions variables and 11 * 4 + 4 = 48 for the power of the continuous variables and the binary variables leading to a total of q = 153 predictors. This number is certainly too large and some of them are clearly irrelevant due to the systematic procedure used to build the non-linear transforms of environmental variables. A LASSO criterion [START_REF] Efron | Least angle regression[END_REF] is therefore used to estimate the regression parameters and to select a subset of significant predictors. This criterion can be written using the notations from Section 2 for one engine variable Y r , r ∈ {1, . . . , p} as :

β r = arg min β r ∈R q I,ni i,j=1 Y r ij - q l=1 β r l X l ij 2 , q l=1 |β r l | < C r (2)
The regression coefficients are penalized by a condition which forces some of them to be null for a well chosen value of C r . The LARS algorithm [START_REF] Efron | Least angle regression[END_REF] is used to estimate all the solutions of the optimization problem (2) for all possible values of C r . The optimal value of C r with respect to the prediction error estimated by cross-validation (with 20 blocs) is finally selected. The number of selected predictors and the coefficient of determination R 2 are listed in Table 2 for all engine variables. Engine variables are well explained by the proposed models as attested by the high value of the coefficients of determination.

A qualitative inspection of the model results was also carried out with the help of engine experts. The regularization path plot (as shown in Figure 3) is very interesting from the point of view of the experts, because it can be compared with their previous knowledge. Such a curve clearly highlights which are the more relevant predictors and they appear to be in very good adequateness with the physical knowledge on the system. In summary, the first preprocessing module (ECN) provides 5 standardized engine variables denoted by S r = [S r ij , i ∈ {1, . . . , I}, j ∈ {1, . . . , n i }], with r ∈ {1, . . . , p}, which are the residuals of the selected regressions. They are independent of environmental conditions but still contain some significant aspects such as linear trends and abrupt changes at specific dates. We therefore propose to use an on-line Change Detection algorithm (CD) together with an adaptive linear model to fit the data.

Change Detection -CD

To take into account the two types of variation (linear trend and abrupt changes), we implement an algorithm based on the ideas from [START_REF] Gustafsson | Adaptative filtering and change detetction[END_REF] and [START_REF] Ross | Online annotation and prediction for regime switching data streams[END_REF]. The solution is based on the joint use of an on-line change detection algorithm to detect abrupt changes and of a bank of recursive least squares (RLS) algorithms to estimate the slow variations of the signals. The algorithm works on-line in order to allows projecting new measurements on the map as soon as new data are available.The method can be described as follows:

1) One RLS algorithm is used for each one of the p standardized engine variables to recursively fit a linear model. For each r ∈ {1, . . . , p}, for each engine i ∈ {1, . . . , I} and at each date l, one has to solve the following equation:

(α r il , β r il ) = arg min α∈R,β∈R l j=1 λ (l-i) (S r ij -(βj + α)) 2 , ( 3
)
where λ is a forget coefficient. The estimates α r il and β r il are respectively the intercept and the slope of the linear relationship. These estimates are then used to define the variables ε r il = S r il -(β r il l + α r il ), which do not contain anymore the slow variations of the signals.

2) These values are concatenated in a vector ε l = [ε 1 l , . . . , ε p l ], which is then used in a multi-dimensional Generalized Likelihood Ratio (GLR) algorithm [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF] to detect the abrupt changes of the signals. The GLR algorithm is a sequential test procedure based on the following model:

ε k ∼ N p (θ(k), Σ), ∀k > 0, with θ(k) = θ ∈ Θ 0 si k < t 0 , θ ∈ Θ 1 si k ≥ t 0 .
3) Finally, when an alarm is sent by the GLR algorithm, all the RLS algorithms are re-initialized. The results supplied by this algorithm are the following:

the alarm dates supplied by the multi-dimensional GLR algorithm; -cleaned signals estimated by the RLS algorithm; -slopes and intercepts estimated by the RLS algorithm.

Figure 4 presents the obtained results for two engines. One abrupt change was found for the first engine and 3 for the second one; all of them seem to be reasonable and a comparison between estimated alarm dates and recorded real events of the engine life have confirmed this fact. The estimated signals are also shown on these two figures. For more information on this aspect of the analysis process see [START_REF] Côme | Aircraft engine health monitoring using self-organizing maps[END_REF]. From now, the observations corresponding to each flight are

F il = [F 1 il , . . . , F p il ]
, where F r il = β r il l + α r il are the results of the transformations of the raw data performed by the first two modules (ECN and CD).

Self-Organizing Maps -SOM

The cleaned signals F il provided by the previous two modules are then used as input to a SOM for visualization purpose. To project the observations we use a [20 × 20] SOM implemented with the Matlab toolbox [START_REF] Vesanto | Som toolbox for matlab 5[END_REF]. Since the variables F r il are correlated, a Mahalanobis distance is used to whiten the data. A classical learning scheme is used to train the map. Figure 5 (a) presents one example of engine trajectories on the map, which clearly have different shapes. For the studied engine, available maintenance reports inform us that this engine suffers from an deterioration of its high pressure core. This fault is visible on the map at the end of the trajectory: the engine which was projected on the middle north of the map during a large part of its trajectory, suddenly moves towards the north-west corner of the map. This area of the map furthermore corresponds to abnormal values of the engine variables.

Similar trajectory matching -SEARCH

One of the final goal of the proposed tool concerns clustering and prediction of engine trajectories or pieces of engine trajectories. For this end, we have to define a proper distance between pieces of trajectories, which can be of different lengths.

Before projection on the map, pieces of trajectories were sequences of R pvectors, but as soon as measurements are projected, they can be described by sequences of integers corresponding to the units where measurements are projected. Such sequences will be denoted by T = [k 1 , . . . , k L ] and as they take their values in a finite set {1, . . . , U } (where U is the number of units), we will call them strings. The difficulty comes from the fact that the strings can have different lengths. Such a problem has been already investigated in other fields and one classical solution is to use Edit Distance, which is commonly used for approximate string matching [START_REF] Navarro | A guided tour to approximate string matching[END_REF].

To compare two strings, Edit Distance uses a cost function. This function gives individual cost for each unitary operation such as: suppression, addition or substitution. The cost of a sequence of operations O = [o 1 , o 2 , . . .] is simply equal to the sum of all the unitary costs, so: cost(O) = t cost(o t ).

Then, the Edit Distance between two strings de(T, T ′ ) is defined as the minimal cost among all sequences of operations that fulfill the constraint O(T ) = T ′ .

Such a distance can be tuned to our end by carefully choosing unitary costs. We may in particular use the map topology to define meaningful substitution costs, by setting the cost of the substitution k ↔ k ′ to the distance between unit k and unit k ′ on the map. With such a choice, we will take benefit of the fact that close units on the map can be exchanged with a small cost. Suppression and insertion costs are equal to the average of all the pairwise distances between units of the map.

With such a distance one can build classes of pieces of trajectories using hierarchical clustering. But this distance can also be used to supply clues on the possible future evolution of one trajectory. To perform such a task, the following method is proposed. Let T be a piece of engine trajectory:

1. compute the Edit Distance between T and all the pieces of engine trajectories recorded in the fleet database [T 1 , T 2 , . . .] (all these distances can be computed efficiently using dynamic programming [START_REF] Navarro | A guided tour to approximate string matching[END_REF]); 2. search for matching pieces T x such that de(T x , T ) < η, where η is a given threshold; Note that these pieces are parts of already observed engine trajectories, which were recorded in the fleet database and their evolutions after the matching points are therefore known, the third step uses this property. 3. look at the pieces that are just after the matching pieces. That gives an idea about the possible future of T and enables the computation of probabilities of different types of maintenance events if the fleet database is connected to the maintenance database which recorded all the failures and maintenance operations performed on the fleet. We hope that it will be a useful tool to anticipate failures.

Figure 5 (b) presents preliminary results obtained using such an approach, T was build using the end of an engine trajectory of length 50. During this time period this engine stay in the same unit. We present the pieces of trajectories that follow the 3 best matching points found in the fleet database using the proposed Edit Distance. These possible futures for T seem reasonable. Further works concern the connection with the mainteance database to perform a quantitative analysis of the results.

Conclusion

The method proposed in this paper is an interesting tool to summarize and represent the temporal evolution of an aircraft engine health flight after flight. The regression approach used to deal with the problem of environmental condition normalization (ECN) seems to be effective. The joint use of an adaptive algorithm to estimate signal evolution (RLS) and of a change points detection method (GLR) is also an interesting solution to deal with the non-stationary of the signals and to clean them (GLR module). Finally, Self-Organizing Maps (SOM) can be used to show the engine health evolution in a synthetic manner and to provide codes for synthetic representation of trajectories, that enables the development of predictive analysis tools (SEARCH module).

Fig. 1 .

 1 Fig. 1. (a) Rough measurements of the corespeed variable as a function of time for engine 6, (b) residuals of the same variable and for the same engine using a simple linear model with the environmental variables and the engine indicator as predictors (see Table1).

Fig. 2 .

 2 Fig. 2. Global architecture of the health monitoring tools.

Table 2 .Fig. 3 .

 23 Fig. 3. Regularization path for the fuelflow variable: regression coefficients evolution with respect to C r . The more significant explanatory variables are given and the best solution with respect to cross-validation is depicted by a vertical line.

Fig. 4 .

 4 Fig. 4. Change detection results for two engines, (a) engine 2, (b) engine 41. Alarms are depicted by vertical lines, input signals are shown in light gray and signal estimates F r using RLS are depicted by a black line. One figure (with respect to egt) is bigger than the others to present more clearly the RLS estimate of the signal.

Fig. 5 .

 5 Fig. 5. Trajectories of engine 22 on the map. The sizes of the dots are proportional to the measurement date: smallest dots correspond to recent measurements, larger dots to older measurements.

Table 1 .

 1 Variables names, descriptions and type.

		Name	Description	Type	Binary
		aid	aircraft id
		eid	engine id
		fdt	flight date
	X 1 ij	temp	temperature	environment
	X 2 ij	nacelletemp nacelle temperature	environment
	X 3 ij	altitude	aircraft altitude	environment
	X 4 ij	wingaice	wings anti-ice	environment
	X 5 ij	nacelleaice	nacelle anti-ice	environment
	X 6 ij	bleedvalve	bleed valve position	environment
	X 7 ij	isolationleft valve position	environment
	X 8 ij	vbv	variable bleed valve position	environment
	X 9 ij	vsv	variable stator valve position	environment
	X ij	hptclear	high pressure turbine setpoint	environment
	X ij	lptclear	low pressure turbine setpoint	environment
	X ij	rotorclear	rotor setpoint	environment
	X ij	ecs	air cooling system	environment
	X ij	fanspeedi	N1	environment
	X ij	mach	aircraft speed	environment
	Y 1 ij	corespeed	N2	engine
	Y 2 ij	fuelflow	fuel consumption	engine
	Y 3 ij	ps3	static pressure	engine
	Y 4 ij	t3	temperature plan 3	engine
	Y 5 ij	egt	exhaust gas temperature	engine