Sabine Bottin-Rousseau

Mikaël Perrut

Christian Picard

Silvère Akamatsu

Gabriel Faivre

An experimental method for the in-situ observation of eutectic growth patterns in bulk samples of transparent alloys

We present an experimental method for the in-situ observation of directional-solidification fronts in bulk samples of transparent eutectic alloys. The growth front is observed obliquely in dark field through the liquid and a glass wall of the container with a long-distance microscope. We show that a focused image of the whole growth front can be obtained at a certain tilt angle of the microscope. At this tilt angle, eutectic fibers of about 3.5µm in diameter can be clearly seen over the whole growth front in 400-µm thick samples.

I. INTRODUCTION

Directionally solidified eutectic alloys form extended growth patterns of spacing in the 1-10µm range for solidification rates on the order of 0.1 µms -1 . The dynamics of these patterns has been a question of fundamental and practical interest for decades. Progress in this field requires experimental methods permitting a spatiotemporal follow-up of growth patterns during solidification. Transparent nonfaceted eutectic alloys, such as CBr 4 -C 2 Cl 6 or Succinonitrile-(d)Camphor (SCN-DC), have been introduced for this purpose a long time ago [1]. Usually, these alloys are directionally solidified in very thin (10-µm thick, typically) samples. This confers a strongly one-dimensional (1D) character to the growth patterns, which can then be observed in side view, i.e. along the direction normal to the sample plane, with a standard optical microscope. This thin-sample method has proved quite accurate [2,3,4], but can naturally be used only in the study of 1D eutectic growth patterns. A method of observation in situ applicable to 2D eutectic growth patterns has long been lacking. Recently, we published observations obtained with a new method permitting a spatiotemporal follow-up of growth patterns in transparent samples of thicknesses ranging from 100 to 500µm [5,6]. In this article, we present this experimental method in detail.

Figure 1 shows the principle of the method. A directional-solidification sample of rectangular crosssection receives light through a wall of the crucible, and is observed with a tilted long-distance (LD) microscope through the opposite wall. The incident light beam and the LD microscope are oriented in such a way that the microscope collects only the light rays that emerge from one of the solid phases, of which the growth front is composed. A dark-field image of, essentially, the sole growth front is thus obtained. An example is given in Figure 2, which shows a fibrous eutectic growth pattern in a 340- * Electronic address: sabine.bottin-rousseau@insp.jussieu.fr µm thick sample of SCN-DC. It can be seen that fibers about 10µm in diameter are well resolved over the whole growth front. The purpose of this article is to explain how such an image can be obtained in spite of the large tilt angle of the LD microscope with respect to the growth direction.

The article is organized as follows. A technical description of the experimental setup is given in Section II. Some data about the transparent alloys used are reported in Section III. The basic imaging properties of the setup are explained in Section IV. The nature of dark-field images is discussed in Section V. Section VI is devoted to the important question of the elimination of thermal bias, i.e. the fine control of the orthogonality between the isotherms and the growth axis.

II. TECHNICAL DESCRIPTION OF THE EXPERIMENTAL SETUP

The experimental setup is composed of three distinct parts (Fig. 3): a directional-solidification (DS) bench, a light source, and a LD microscope (Questar QM100) equipped with a numerical black-and-white video camera (Scion Corporation). The images are stored in a personal computer for further numerical treatment (contrast enhancement, and "rescaling", i.e. stretching of the captured image in the direction y so as to restore its width to the actual value w; see Fig. 2). The DS bench is put vertically, with the direction z pointing upwards. The LD microscope and the light source are mounted on rotation stages of axes parallel to x. The different parts of the setup are mounted on independent translation stages, which are fixed on an anti-vibrating table.

FIG. 3: Technical sketch of the experimental setup. DS: directional-solidification bench (see Fig. 5). V : dc motor. XYZ: motorized xyz translation stages. LDM: long-distance microscope (the video camera is not shown). M1 to M4: step motors. M1: control of the working distance. M2: rotation of the microscope. M3 and M4: automatic adjustment of focus. AVT: anti-vibrating table. L: linear fiberoptic array (the halogen lamp is not shown). The reference frame is the same as in Fig. 1 A. Directional-solidification bench

The DS bench has been designed for low melting point transparent alloys. In this study, both the CBr 4 -C 2 Cl 6 and SCN-DC eutectics were utilized. The crucibles are made of two 300-µm thick glass plates (a "base" and a "cover" plate) separated by plastic (PET) spacers (Fig. 4). The spacers are cut to a special funnel shape in order to form a grain selector. The cover plate and the spacers are fixed to the base plate with a UV polymerisable glue (Norland NOA81). A pressure is applied onto the cover plate during the gluing process. The thickness w of the inner volume depends on both the spacer thickness and the applied pressure. We determine w by measuring the outer thickness of the sample at the end of the process. The crucibles are heated above the liquidus temperature of the alloy, filled with liquid alloy by capillary action under a controlled Ar atmosphere (≈ 10 mbar), and sealed after a rapid cooling to room temperature. The alloy is then entirely solid. The base glass plate fits in a rigid holder, which is inserted into the DS bench and connected to the translation motor. The DS bench is composed of two ovens fixed on an insulating soil, and a translation motor. Each oven is made of a fixed (base), and a removable (cover) copper block, between which the sample is inserted (Fig. 5). The blocks are machined with V-shaped edges to permit oblique observation and illumination through the gap separating the ovens. The base glass plate of the crucible is in contact with the base blocks of the ovens. A pressure is exerted onto the cover blocks by means of a screw-and-spring system (not shown) in order to ensure a better thermal contact and an more accurate guiding of the sample translation. The temperature field is established by heat diffusion along the sample. The temperature gradient G at the growth front can be varied from 40 to 110 ±10% Kcm -1 by varying the width of the observation window. Temperature regulation is ensured by water circulation for the cold oven, and by resistive sheets connected to an electronic PID controller for the hot oven. The temperatures of the base and cover hot block are regulated independently of each other to permit the adjustment of the y component of the thermal gradient (see Section VI). A dc motor (Physik Instrumente) ensures the translation of the sample at velocities ranging from 0.01 to 10 µm.s -1 with an accuracy of ±1%.

B. Optical system

The Questar QM100 system is a Maksutov-Cassegrain catadioptric telescope utilizable at working distances ranging from 15 to 35 cm; at the working distance of 15 cm, the numerical aperture of the microscope is of 0.14, the diffraction-limited resolution is of 1.1µm and, and the depth of field is of 28µm(information from the constructor). We usually set the working distance at 15 cm. The field of view of the complete instrument (LD microscope plus video camera) is 440-µm wide in the direction x. The whole growth front can be scanned along x using the motorized translation stage of the DS bench (Fig. 6).

The tilt angle θ t of the microscope axis, called "observation angle", is measured from the horizontal (sideview) position. It can be varied continuously from 0 to 60 o . An automatized mechanical device comprising two translation motors keeps the field of view and the focus unchanged as θ t varies. An additional translation stage serves to vary the working distance, i.e. to zoom or unzoom the image. An iris diaphragm is situated in a plane close to the entry of the LD microscope. Incidentally, a slight misalignment of the microscope axis with respect to the plane yz forced us to decenter this diaphragm in the direction x in order to obtain a satisfactory dark-field contrast.

The light source is a 30 × 2mm linear fiberoptic array aligned on the direction x illuminated by a halogen lamp. It can be rotated around x, and displaced in the directions y and z. For a horizontal distance between the source and the object of 5 cm, the angular aperture of the incidence beam is of about 2 o in the plane yz, and 30 o in the direction x . We adjust the incidence angle θ o by changing the vertical position of the source. Simultaneously, we rotate the source in order to keep the intensity of the incident light at its maximum.

III. TRANSPARENT NONFACETED EUTECTIC ALLOYS

The two eutectic solid phases of the CBr 4 -C 2 Cl 6 system, denoted α and β, are solid solutions of fcc and bcc structures, respectively. The two eutectic solid phases of the SCN-DC system are practically pure bcc SCN , and hexagonal DC, respectively, and will therefore be denoted "SCN" and "DC". The concentrations and temperatures at the eutectic points in these two systems are given in Table I. More details about the phase diagrams can be found elsewhere [7,8]. Here we are interested in the refractive indices of the eutectic phases, especially, the differences in refractive index between eutectic solids and liquid, in these systems. These refractive indices have never been measured, to our knowledge, so we attempted to calculate them using the Lorentz-Lorenz formula [START_REF] Hartshorne | Practical Optical Crystallography[END_REF]. The ingredients of the calculation are the molecular refractivity of the components, and the density of the phases of interest. We calculated the molecular refractivities of the components from the measured values of the refractive indices of pure CBr 4 , C 2 Cl 6 , SCN and DC available in the literature [START_REF]Handbook of Chemistry and Physics[END_REF]. The densities of the pure substances at the temperatures of interest could also be calculated, or estimated with good accuracy from data available in the literature. Fortunately, all the phases of interest except one (the β phase of CBr 4 -C 2 Cl 6 ) have a low concentration in one of the components. We estimated their density (and hence their refractive index) assuming that the molar volumes in these phases were the same as in the pure substances (Table I). This assumption is correct in the case of the SCN and DC phases, reasonable in the case of the two liquids and the α phase, but questionable in the case of the β phase. The calculated difference in refractive index is -.005 for the SCN-liquid pair, 0.087 for the DC-liquid pair. This strong difference in order of magnitude is supported by experimental observations, which show that, in thin samples, the contrast of the DC-liquid interfaces much stronger than that of the SCN-liquid interfaces [4]. Similarly, experimental observations in thin samples of CBr 4 -C 2 Cl 6 show a strong contrast of the α-liquid interfaces and a weak contrast of the β-liquid interfaces [3], indicating that the refractive index of β in Table I is most probably largely overestimated.

IV. IMAGING PROPERTIES A. Optimal observation angle

Two basic properties of the experimental setup are the existence of an optimal observation angle, at which an infocus image of the whole growth front can be obtained, and a relatively small astigmatism aberration. To explain these properties, it will suffice to consider the propagation of the light rays in the plane yz (Figure 7a). We consider a rectilinear object (a planar growth front) im-

Material TE Concentration density Refractive θ ⊥ f ⊥ R ⊥ (C) (mol%) (g/cm
The equation for the light ray in the external medium reads

z = z o + py, (4.2)
where

z o = L 1 (p 1 + p b ) + L 2 p 2 , (4.3)
L 1 is the distance of m from the cover plate, and p b is the slope of the growth front due to thermal bias. The intersection of the light rays of slopes p and p+δp emitted by m is given by

z = z o (L 1 , p) + py = z o (L 1 , p + δp) + (p + δp)y. (4.4)
The image m i = (y i , z i ) of m corresponds to the limit δp → 0. One has

y i = -∂z o /∂p | p=pt , (4.5)
and z i = z o + p t y i , where p t = tanθ t . The image of the growth front is a curve of slope Since p i is independent of L 1 , the image of a planar front is a plane. A sharply focused image of the whole front can be obtained only when the image plane is normal to the axis of the microscope, i.e. when p i p t = -1. This condition defines the optimal observation angle θ ⊥ . By setting p i = -p -1 t in eq. 4.6, we obtain the following equation The contraction factor for this value of θ is I.

p i = ∂z i /∂y i | p=pt = p t - p 1t + p b ∂p 1 /∂p t . (4.6)
f ⊥ = n -1 1 sin 2 θ ⊥ . ( 4

B. Astigmatism

It is a well-known fact that the astigmatism aberration of a planar diopter is an increasing function of the distance of the object from the diopter and the aperture of the optical setup. This is the reason for which we took a thin cover glass plate and introduced an additional entry diaphragm. Indeed, the images obtained with the sole LD microscope (ǫ ≈ 8 o at a working distance of 15 cm) had a blurry aspect, which disappeared when ǫ was reduced by means of the entry diaphragm (Fig. 8). Reasonably sharp images were obtained for ǫ = 4.4 o and w ≈ 400µm (see Fig. 2) indicating that the astigmatism aberration of the setup was comparable to the diffraction-limited resolution of the LD at microscope (2.2 µm) at this value of ǫ.

We tested the validity of these empirical conclusions using eq. 4.5. According to this equation, the shift from the sharp-focus point due to a deviation δp from the chief ray is δy i = -1/2 ∂ 2 z o /∂p 2 δp along the axis y, and δy i × δp/ cos θ along the axis Y . Thus the broadening ∆Y of the image of m due to astigmatism is

∆Y ≈ | ∂ 2 z o /∂p 2 | 2 cos 3 θ ǫ 2 . (4.11) Assuming that n 2 ≈ n 1 , we obtain ∆Y ≈ R ⊥ (L 1 + L 2 )ǫ 2 . (4.12)
where R ⊥ = 1.5 n -1 1 sin 2 θ ⊥ . For n 1 = 1.451, L 1 + L 2 = 750 µm, and ∆Y = 2.2 µm, we obtain ǫ = 4.0 o , in reasonable agreement with the above conclusions. On the other hand, observations showed that, still for ǫ = 4.4 o and L 1 < 400µm, the image quality remained tolerable for deviations of the observation angle from the optimum value as large as 10 o . For this deviation and L 1 = 400µm, the maximum distance between a point of the image and Y is of about 70µm, which is comparable to the diffraction-limited depth of field of the LD microscope (90µm) at the value of ǫ under consideration .

V. DARK-FIELD IMAGES

When the solid and the liquid are homogeneous, only the light rays that intersect the solid-liquid interface undergo a deviation. This deviation is a function of the difference in refractive index between solid and liquid, and the slope angle ψ of the interface. At fixed observation angle, a dark-field image is obtained by adjusting the incidence angle so as to capture only the light rays that suffer a given deviation as they pass through the sample. Thus, dark-field images essentially reflect the composition distribution, and the topography of the solid-liquid interface. An example of a purely topographic contrast is provided by Figure 9, which shows a dark-field image of a single-phased solidification front. In this experiment, θ o was adjusted for the microscope to capture the light emitted by an interface normal to z. Non-deformed regions of the front appear bright, and distorted regions (grooves associated with grainboundaries) dark. Satisfactorily, the image width of large-angle grainboundary grooves is comparable to the known physical width (≈ 12µm) of such grooves in CBr 4 -C 2 Cl 6 [START_REF] Bottin-Rousseau | [END_REF]. The contrast as-sociated with low-angle grainboundaries is weak in accordance with the small amplitude of the corresponding grooves.

In the case of eutectic patterns, both the differences in refractive index between the different phases, and the strong local deformation of the interface contribute to the image contrast. Let us consider specifically quasistationary fibrous eutectic patterns (Fig. 7b). We first note that, when a pattern is stationary, the solid-solid interfaces left behind in the solid are parallel to z so that refraction at these interfaces does not affect the angle between the emergent light rays and z. Consequently, the solid behind the interface appears totally dark in spite of multiple refraction by solid-solid interfaces. Second we remark that, in the type of dark-field images illustrated in Fig. 2, each DC fiber gives rise to a diffuse light patch ("bright spot"), which does certainly not coincide with the interface between the fiber and the liquid ("fiber tip"). Let us discuss the nature of these bright spots in some detail.

In the case of SCN-DC, an inversion of the contrast between fibers and matrix is observed as a function of θ o at fixed θ t (Figure 10). The maximum of brightness is at about θ o = θ t -2 o for the SCN matrix, and θ o = θ t + 1 o for the DC fibers. Figure 11 shows the calculated value of θ o -θ t as a function of ψ at fixed θ t for the two solid phases. The calculated values are in reasonable agreement with the experimental ones for ψ ≈ 60 o , thus for a solid-liquid interface roughly normal to the light rays. On the other hand, we compared directly the area of the bright spots with the theoretical value A DC of the fiber cross-section. In a regular hexagonal pattern,

A DC = η √ 3λ 2 /2
, where η ≈ 0.23 is the equilibrium volume fraction of DC in the solid, and λ is the spacing. Figure 12 shows enlarged views of a series of SCN-DC fibrous patterns at average spacings ranging from 22 to 6 µm, on which white disks of areas A DC (λ) have been superposed. It can be seen that the area of the bright spots, although it is ill defined, is clearly smaller than A DC , at least when λ is sufficiently large. We conclude that the bright spots are caustics, which are generated by the rounded tips of the fibers, and are concentrated around a direction making a small angle with the direction of the incident beam. Fig. 12 also shows that the ratio of the diameter of the bright spots to the actual fiber diameter increases as the fiber diameter decreases. This is due to the diffractioninduced broadening of the bright spots. The value of the fiber diameter, at which the ratio is unity, can be considered as the empirical resolution power of the instrument. From Fig. 12, we estimate this resolution power to be of about 3.5µm, which is somewhat larger than the diffraction-limited resolution of the LD microscope (2.2 µm). Previous studies demonstrated a similar empirical resolution power for CBr 4 -C 2 Cl 6 lamellar eutectic patterns [5,6].

VI. THERMAL-BIAS CORRECTION

In a perfect directional-solidification experiment, the isotherms would be planar and perpendicular to the solidification axis. In any real experiment, the isotherms are both curved and, on average, tilted with respect to the solidification axis (thermal bias). The fact that the dynamics of solidification patterns is very sensitive to these instrumental imperfections has not always received sufficient attention. A thermal bias progressively drives the structure along the direction of the bias, while a curvature causes a progressive blow-up (when the isotherms are bulging outward from the solid, which is our case) of the structure over time. This will be the object of a detailed report in the near future [12]. Here we content ourselves with stressing the importance of these effects [13].

We determined the curvature of the isotherms in our setup by following the dynamics of fibrous eutectic patterns over long periods of time. We found that the curvature was negligible in the direction x, but not in the direction y. The radius of curvature in the direction y ranged from 2 to 6 mm in 400-µm thick samples, depending on the experiment, which corresponds to a difference in z between two points located at y = 0 and y = w/2, respectively, ranging from 3 to 10µm. There is apparently no way to correct this curvature, which mostly arises from the discontinuities in heat conductivity at the solid-liquid, solid-wall and liquid-wall interfaces.

We determine the thermal bias at the onset of an experiment by observing the solid-liquid interface at rest (V = 0) in side view, and correct it by changing the temperature of the hot cover block in the DS bench. The details of the procedure are as follows. Previous to a solidification run, a sample undergoes a partial directional melting, after which it is held at rest for equilibration. When the annealing is completed, the concentration C L in the liquid is homogeneous, and the shape of the solidliquid interface is that of the isotherm at the liquidus temperature T L of the alloy. (In the case of a eutectic alloy, T L is the liquidus temperature for the primary eutectic phase). At this stage, we generally observe a noticeable transverse thermal bias, which we measure using focal-series reconstructions (Fig. 13). The thermal bias is given by ψ b = arctan(∆z/w), where ∆z is the difference in height between the intersections of the interface with the two glass plates. The residual thermal bias after correction does not exceed 2 o . The validity of the procedure can be crosschecked in oblique view, by comparing the measured value of the contraction factor to the theoretical value (eq. 4.10).

VII. CONCLUSION

We have tested, and validated, a new method of observation of transparent eutectic growth fronts in bulk samples. The method makes use of three basic ingredients: a microscope with a very long working distance and a narrow aperture, an oblique direction of observation, and a dark-field contrast relying on specific features of eutectic patterns, namely, the difference in optical index between the different solid phases, and the strong local curvature of the interface. With a prototype instrument based on this method, we obtained satisfactory images of eutectic patterns of spacing no larger than 3.5µm in 400-µm thick samples. Much larger sample thicknesses could probably be utilized with more sophisticated instruments based on the same method. More importantly, the principle of a long-distance oblique observation (without the dark-field specification) could be applied to other objects than eutectic growth patterns, for instance, grain boundaries intersecting the single-phased growth fronts, as shown above, and thus become a standard method in solidification science.

FIG. 1 :

 1 FIG.1: Oblique observation of a transparent eutectic directional-solidification front. z: direction of the thermal gradient. The sample is pulled in the direction -z at a rate V . y: normal to the sample walls. The optic axis of the microscope lies in the plane yz.

FIG. 2 :

 2 FIG. 2: Upper panel: Oblique view of a fibrous eutectic pattern in a 340-µm thick sample of eutectic SCN-DC. V = 0.035µms -1 . Bright spots correspond to DC fibers, while the SCN matrix appears dark. The tilt angle of the LD microscope is about 40 o . Lower panel: rescaled image. Horizontal dimension: 1300µm.

  FIG. 4: a) Sketch of a crucible (glue seals are not shown). b) Cross-section of the grain selector.

FIG. 5 :

 5 FIG.5: Directional-solidification bench. Each oven is made of a base (dark grey) and a cover (soft grey) copper block. A sample fixed in its holder (mid grey) is represented. The connection to the motor is not shown. Dash-dotted lines: rotation axes of the cover blocks.

  FIG.6: Whole-front mosaic micrograph of the same growth front as in Fig.2. Note the lamellar growth patterns near the edges of the sample. Horizontal dimension: 4800µm.

FIG. 7 :

 7 FIG.7: Sketches illustrating the processes under study. a) Large-scale image of a planar growth front. (The direction of y is opposite to that defined in Fig.1). b) Origin of dark-field contrast in the case of a fibrous SCN-DC pattern.

  for θ ⊥ : p(p 1 + p b ) = (1 + p 2 )∂p 1 /∂p. (4.7) For an arbitrary value of θ, the contraction factor of the image in the direction Y normal to x and the microscope axis is f = ∂z o /∂L 1 cos θ cos ψ. (4.8) From eq. 4.3 and the definition of p, one obtains f ∝ (p 1 + p b )(1 + p 2 ) -1/2 , and hence ∂f /∂p ∝ ∂p 1 /∂p(1 + p 2 ) -(p 1 + p b )p. Comparing this formula with eq. 4.7 shows that f is maximum at the observation angle θ ⊥ . It is useful to note that, when the growth front is normal to z (p b = 0), eq. 4.7 has the simple form tan θ ⊥ = (1 -n -2 1 ) -1/4 . (4.9)

. 10 )

 10 These functions vary slowly with n 1 , except when n 1 is close to unity. Their estimated values at the eutectic points of the CBr 4 -C 2 Cl 6 and SCN-DC alloys are displayed in Table

FIG. 8 :

 8 FIG.8:A fibrous pattern in a 350-µm thick CBr4-C2Cl6 sample viewed without (upper panel) and with (lower panel) an additional diaphragm. V = 0.25µms -1 . Rescaled images. Horizontal dimension: 750µm.

FIG. 9 :

 9 FIG. 9: Solidification front of a α crystal in a 400-µm thick sample of a dilute CBr4-C2Cl6 alloy. V = 2 µms -1 . Nonrescaled image. A dark zigzag line runs along the intersection of a grainboundary with the solid-liquid interface. Angular points on this line are junctions with subboundaries. Dark spots are liquid droplets migrating in the solid. The residual brightness appearing under the interface is due to multiple diffusion by these droplets. Horizontal dimension: 425 µms -1 .

FIG. 11 :

 11 FIG. 11: Calculated deviation of a light beam intersecting a SCN-, or a DC-liquid interface as a function of the slope angle of the interface for θt = 40 o .

FIG. 13 :

 13 FIG. 13: Side view of a 400-µm thick sample of a dilute CBr4-C2Cl6 alloy at rest. Upper panel: A focal series before thermal-bias correction. Horizontal dimension of a micrograph: 130µm. Left to right: the focus changes from the base, to the cover glass plate. In this example, ∆z = 83µm and ψ b = 12 o . Lower panel: Same interface after thermal-bias correction. The focus is on the cover glass plate. Horizontal dimension: 420µm.

TABLE

VIII. ACKNOWLEDGMENTS

We thank V.T. Witusiewicz, L. Sturz, and S. Rex from ACCESS (Aachen, Germany) for kindly providing us with purified Succinonitrile-(d)Camphor alloy. This work was supported by the Centre National d'Etudes Spatiales, France.


Seungwoo Kang

Sandrine Selosse

Nadia Maïzi

WHAT DEVELOPMENT FOR BIOENERGY IN ASIA: A LONG-TERM ANALYSIS OF THE EFFECTS OF POLICY INSTRUMENTS USING TIAM-FR MODEL

Keywords: bioenergy, climate change, action plan, strategy, policies

Asian countries have sharply increased their energy consumption with economic growth, inducing an extremely high dependency on energy fossil due to a low level of energy resources. In 2011, Asia's share in global energy consumption reached 40%, of which China, India, Japan and South Korea represented 70%. Furthermore, these countries are ranked among the world's largest CO 2 emitters. Hence, bioenergy is currently being highlighted to diversify the energy mix and to mitigate climate change in these countries. Along with the direct bioenergy promotion policies, a variety of economic policy instruments are introduced as emission trading scheme, carbon and energy taxes. In this study, we estimated the evolution of bioenergy deployment from current energy and climate policies through scenarios analysis using bottom-up energy system optimization model, TIAM-FR, a TIMES family model from ETSAP/IEA. Different scenarios were developed with basis on political instruments in these countries:

(1) BAU, (2) Global climate change (50% of world GHG emissions reduction), (3) National GHG mitigation targets from INDCs, (4) National GHG mitigation targets + Bioenergy policy, (5) Carbon tax system. According to the results, the higher is climate constraints, the higher is bioenergy development. Hence, carbon tax system derived the most bioenergy consumption by 2030 followed by global climate change scenario. Current INDCs target is evaluated as not sufficient to promote bioenergy development without renewable energy targets and biofuel blending mandates.

INTRODUCTION

During the last decades, Asia countries have drastically increased their energy consumption due to fast economic developments and growing population. Currently, Asia is highly dependent on fossil energy and it led this region to be ranked at top Greenhouse gas (GHG) emitter in the world. Then, growing concerns about climate change, bioenergy attracts more and more attention as one of the best alternative sources of energy to fossil energy, as currently shown in the number of ambitious bioenergy promotion policies implemented in Asia. However, the relation between bioenergy deployment and several political instruments in harmonized energy system still needs to be addressed. Therefore, this study aims to observe how the different political tools encourage or discourage bioenergy development and especially for the largest contributors to GHG emissions in Asia, i.e. China, India, Japan and South Korea.

METHODOLOGY

The evaluation of bioenergy's future is performed through scenario analysis. Several different scenarios have been analyzed and compared to current bioenergy policies of four Asian countries, China, India, Japan and South Korea, for the period 2005-2050. This long-term analysis is performed with a bottom-up energy system optimization model, TIAM-FR, developed by the MINES ParisTech's Centre for Applied Mathematics [START_REF] Selosse | A regional perspective to achieve the future climate regime : a long-term analysis with the TIAM-FR model[END_REF].

Presentation of TIAM-FR

TIAM-FR is the French version of ETSAP-TIAM, the world TIMES model developed under IEA's ETSAP (Energy Technology Systems Analysis Program). This model is based on a bottom-up approach and a technology-rich representation of the energy system to estimate its changes and evolutions for the long-term.

The model covers the entire world in 15 regions, of which China, India, Japan and South Korea are separately represented. And the model has 5 energy service sectors (agriculture, industry, commercial, residential, transport) whose demands are based on socio-economic assumptions (GDP, household, population, sectorial growth, etc.). Technico-economic characteristics are integrated into the model to describe several technologies in all sectors of each stage of the energy chain (extraction, conversion, processing, transmission, and end-uses) in each region (see figure 1). Where NPV is the net present value of the total cost for all regions over the projected period; ANNCOST (r,y) is the total annual cost in region r and year y; dr,y is the discount rate; REFYR is the reference year for discounting; YEARS is the set of years and R is the set of regions (15 regions).

Results of the optimization are the structure of the energy system for each region, i.e. type and capacity of the energy technologies, energy consumption by fuel, development of emissions, energy trade flows between the regions as well as the therefore required transport capacities, and detailed energy system costs plus information on the marginal costs of environmental measures, etc.

Energy-climate Scenarios

In this study, several different climate scenarios were explored on four countries (China, India, Japan and South Korea). A business as usual (BAU) scenario is firstly developed and calculated. This baseline scenario without any emission constraint outlines some key patterns in the evolution of the energy system, and serves as the starting point for comparing other different avenues of carbon mitigation over the period 2005-2050.

Secondly, Global_Factor2_GHG scenario, which consists in reducing 50% of GHG emissions by the year 2050 compared to the year 2005.

From the recent decision at COP 20 (Conference of Parties) in Lima, Peru, each country is invited to communicate its GHG mitigation action plan by reporting Intended Nationally Determined Contributions (INDCs) in advanced of the COP 21 in Paris. Currently, China, India, Japan and South Korea have submitted their INDCs report to the UNFCCC secretariat [START_REF][END_REF]. According to these reports, it comprises national GHG emission reduction target until 2030 and action plans as renewable energy and biofuel supply targets. Based on INDCs reports [START_REF]Japan, Intended Nationally Determined Contribution[END_REF][START_REF]India, Intended Nationally Determined Contribution (INDC) -India[END_REF][START_REF]China, Intended Nationally Determined Contribution[END_REF][START_REF]Republic of Korea, Intended Nationally Determined Contribution (INDC) -Republic of Korea, Unfccc[END_REF], two scenarios are developed in this study.

First of them, called "INDC", applies only the GHG emission reduction targets until 2030 to investigate the bioenergy evolution without other mandates. Second scenario, called "INDC+RNW+BIOFUEL", combines GHG emission targets with renewable and bioenergy supply targets. The table 1 presents the different GHG emission targets applied to this scenario:

Table I: GHG emission reduction targest applied to scenario ( recalculated INDCs target by author)

In case of Japan, announced GHG emission reduction target is a reduction of 26% by 2030 compared to 2005 level. However, the Japanese government plans to use 50 -100 million tCO2eq of JCM (Japanese Crediting Mechanism), which achieves 4-8% of their target. In addition, we applied the target of GHG emission related to energy consumption, which is published as 24% of reduction by 2030 in INDCs. Hence, we applied recalculated target of 16-20% of reduction by 2030 compared to 2005 level to our scenario.

South Korea announced their target to reduce 36% of BAU GHG emissions by 2030. However, the ministry of environment of South Korea communicates in official press that 11.3% emission reduction will be achieved by international carbon crediting system. Hence, we applied only domestic GHG emission targets of 25.7% reduction to our scenario. [START_REF] Choi | Korea confirmed 2030 greenhouse gas reduction target of 37%[END_REF] Different renewable energy targets and biofuel blending rates are explained in table 2. Particularly, in case of Japan, maximum blending rates of bioethanol and biodiesel are implemented for avoiding possible degradation of vehicle performance.

Table II: Renewable energy targets and biofuel blending mandates

In addition to these direct GHG emission constraints, economic instrument to reduce GHG emission, a carbon tax system, was analyzed in this study. We developed two different levels of carbon tax scenarios described in table 3. For low level of carbon tax, 80$/tCO2 was applied by 2020 and 100$/tCO2 by 2030. Then the carbon tax level progresses to 270$/tCO2 until 2050. On the other hand, under high level of carbon tax scenario, tax rate is doubled by 2030 as 200$/tCO2 compared to low level scenario. 3 RESULTS

Final energy consumption: bioenergy

In this study, we calculated each scenarios until 2030 for each countries (China, India, Japan and South Korea).

The first option to response to theses climate constraints in these countries was transition from coal to gas and other less GHG emitting energies. Furthermore, the results state that stricterclimate constraints derive more bioenergy consumption share in final energy consumption. In all four studied countries, high level of CO2 tax system, which is the most restrictive climate scenario, promotes more bioenergy consumption than other scenarios. In case of India, biofuel consumption start to appear since 2020 under INDCs scenario combined with 20% of biofuel blending mandate. However, India's voluntary GHG target can not enoughly drive bioenergy consumption than BAU scenario. For example, INDCs scenario derives between 14% ~ 15% by 2020 while 15% of bioenergy consumption shows up under BAU scenario by 2020. By strengthening climate constraiants by 2030, bioenergy consumption becomes more sigfinifcant under climate scenarios than BAU scenario. While BAU scenario shows 11% of bioenergy consumption, intensified climate scenarios promote it up to 22% from 12%. In case of fossil energy consumption, India's INDCs pledges allow still significant quantity of coal consumption, about 30% of total final consumption, due to the low restriction of GHG emissions (see figure 3). In Japan, the share of bioenergy in final energy consumption is insignificant during all projected years, which represents only between 0.4%~2.5% (see figure 4). In fact, Japan is the only country, who set a maximum blending rate of biofuels for the reason of vehicle performance and security issues. Hence, this regulation limits promotion of biofuels in transport sector. Bioenergy evolution in South Korea is not dissimilar to Japan's situation. Bioenergy is not the first promising option to reduce GHG emissions. However, tight climate scenarios since 2030 brings up bioethanol and biodiesel consumption more than required minimum blending rate. Share of biofuel in final energy consumption increases until 3% (alcohol 2.3% and biodiesel 0.7%) at most by 2030 (see figure 5).

Power supply sector

According to scenario results (see figure 6), energy transition from coal to gas appeared as the first considered option. Then, coal CCS technology dominates the electricity generation with tighter climatic constraints as "high level of carbon tax system" and "Global Factor2 GHG" scenarios. By 2030, coal CCS technology in power supply sector represents about 57% for China, 58% for India and 39% for Japan under higher carbon tax scenario. Meanwhile, coal CCS technology does not come up with South Korea but only the first phase of energy transition from coal to gas occurs until 2030 and fossil energy still dominates energy mix for power supply. However, higher climate constraints lead to develop gas CCS technology for South Korea by 2050.

In energy mix of power supply sector, renewable energies and CCS technologies are in competition. And we observed that CCS technologies become more representatives with higher climatic constraints and the share of renewable energies, including bioenergy, decreases. For example, the share of bioenergy drops to almost 0% with higher carbon tax system in all four studied countries by 2030 while "INDC+RNW+BIOFUEL" scenario produces 22% of electricity from biomass in India, 7.3% in China, 5% in Japan and 1.6% in South Korea.

GHG emission reduction: INDCs

In this section, we studied GHG emission reduction of each climate scenario with a focus on announced targets of GHG emission reduction in INDCs reports. In case of China and South Korea, the announced GHG targets are situated near the emission levels of "carbon tax system" and "global Factor2 GHG" scenarios. In 2030, China's GHG emissions increases to 7.5 GtCO 2eq with INDCs target as well as 7.1GtCO 2eq with global GHG reduction scenario by 50% (see figure 7). On the other hand, the targets of GHG emission reduction by 2030 in Japan and India derive the GHG emissions, far away from global GHG emissions reduction target, and remains on BAU level. In 2030, Japan's GHG emissions with INDCs target (1.04 GtCO 2eq ) exceeds about 37% over Global Factor2 scenario emissions (0.75 GtCO 2eq ) and achieve only 9% of reduction from BAU level (1.14GtCO 2eq ) (see figure 9).

Figure 9: GHG emission projection for Japan

This non conformity with global Factor2 GHG scenario intensifies in India. Current carbon intensity reduction pledge of India brings GHG emissions up to 3.2 GtCO2eq by 2030, which exceeds about 62% over allowed GHG emissions under global Factor2 GHG scenario (1.9 GtCO2eq). And India achieves only 10% of GHG reduction from their BAU GHG emissions level (3.6GtCO 2eq ) (see figure 10). In all four studied countries, bioenergy consumption increased with different climate constraints. As a result, the tighter are climate constraints imposed, the more is bioenergy developed. However, as we identified in power supply sector, CCS technology becomes dominant from a certain level of climate constraint. In this case, biomass' share decreases in total energy mix by replacing with more economic and less GHG emitting technologies as Coal and Gas CCS technologies. Hence, the results show that higher climate constraints do not always guarantee more bioenergy development and it is required to investigate optimal level of climate constraints that may maximize bioenergy development.

In terms of GHG emissions, INDCs GHG emissions targets of China and South Korea approached closely to global GHG reduction target of factor 2. However, current INDCs GHG target of India and Japan are evaluated quite generous to achieve global GHG mitigation target based on out assumptions, for example, the projected evolution of GDP, population, technology progress and many other different factors.

In conclusion, all calculated scenarios do not give an important advantage in bioenergy development apart from mandated minimum bioenergy share in energy mix either of biofuel in transport or of biomass in power supply. However, we identified that enough level of climate constraint increases the bioenergy share in energy mix without other supporting tools for bioenergy before the transition to CCS technology commences. So, for further promotion of bioenergy in this region, it is recommended to investigate in finding out a proper level of climate constraints and effective combination with direct imposing policy on bioenergy.

Figure 1 :

 1 Figure 1: Energy system structure in TIAM-FR This model is based on the optimization theory, which aims to minimize the discounted global energy system cost over the entire model horizon until 2050. The model calculates the net present values of total cost for each region with the objective function as follows:

Figure 2 :

 2 Figure 2: Finale energy consumption until 2030 : ChinaIn China, bioenergy consumption share increases up to 22.9% by 2030 with high level of carbon tax scenario although it remains at 11.4% with BAU scenario. On the contrary, coal consumption disappears in the first reaction to climate-energy constraints by being replaced with gas (see figure2).

Figure 3 :

 3 Figure 3: Finale energy consumption until 2030 : India

Figure 4 :

 4 Figure 4: Finale energy consumption until 2030 : Japan

Figure 5 :

 5 Figure 5: Finale energy consumption until 2030 : South Korea

Figure 6 :

 6 Figure 6: Power supply energy mix in 4 countries by 2030

Figure 7 :

 7 Figure 7: GHG emission projection for ChinaAlso, South Korea's GHG emissions reach 0.45 GtCO 2eq with INDCs target, which is similar to global GHG Factor2 scenario's estimate GHG emissions of 0.42 GtCO 2eq (see figure8).

Figure 8 :

 8 Figure 8: GHG emission projection for South Korea

Figure 10 :

 10 Figure 10: GHG emission projection for India

Table III :

 III Carbon tax levelsIn our study, carbon tax scenarios are independently represented, which are not combined with other scenarios. The stand-alone carbon tax scenarios are calculated to estimate the effectiveness on GHG emission reduction and bioenergy development.

			Country		Target
				15% of non-fossil in TPES by 2020 and
				20% by 2030
			China	10% minimum blending rate of bioethanol
				in gasoline by 2015
				40% of power supply capacity from
			India	renewables in 2030 20% minimum blending rate of biofuel
				Solar 7%, Wind 1.7%, Geo 1.1%, Hydro
				9.2%, Biomass 4.6% by 2030 [3]
			Japan	maximum blending rate (3% bioethanol,
				5% biodiesel) [8]
				4.6% renewables by 2029 [9]
			S. Korea	3% minimum blending rate of biodiesel by
				2020; 5% by 2030 [10]
	Country	Target
		Reduce CO 2 emissions per unit of GDP by
	China	60-65% in 2030 compared to 2005 level
		Reduce GHG emissions per unit of GDP
	India	by 33-35% in 2030 compared to 2005 level	Unit : $/tCO2eq	2020	2030	2050
			Low Tax		80	100	270
	Japan	Reduce GHG emissions level 16-20% in 2030 compared to 2005 level	High Tax		150	200	530
	S. Korea	Reduce BAU GHG emissions level by 25.7 in 2030

International Bioenergy (Shanghai) Exhibition and Asian Bioenergy Conference 2015, 21-23 October 2015, Shanghai, P.R. China

ACKNOWLEDGEMENTS

This research was supported by the Chair Modeling for sustainable development, driven by MINES ParisTech, Ecole des Ponts ParisTech, and Agro ParisTech, supported by ADEME, EDF, GRTgaz and SCHNEIDER ELECTRIC.


Pierre Raphaël Bertrand

Doha Hadouni

Change point detection by Filtered Derivative with p-Value : Choice of the extra-parameters

Keywords: Change points detection, Filtered Derivative with p-Value, Filtered Derivative extra-parameters, Mean Integrated Square Error (MISE), Impact on MISE

Change point detection by Filtered Derivative with p-Value : Choice of the extra-parameters

Pierre, Raphael Bertrand,

Introduction

Change point detection is an important problem in various applications: signal processing [START_REF] Boutoille | Systèmes de fusion pour la segmentation hors-ligne de signaux gpd multiporteuses[END_REF], global warming [START_REF] Wang | Testing for increasing weather risk[END_REF], magnetospheric dynamics [START_REF] Wanliss | Space storm as phase transition[END_REF], neuro-physiological studies [START_REF] Schneider | Messages of oscillatory correlograms : A spike train model[END_REF][START_REF]Unitary events in multiple single neuron spiking activity: I. detection and significance[END_REF][START_REF] Grun | Unitary events in multiple single neuron spiking activity. II. non-stationary data[END_REF], motion of chemical or physical particles [START_REF] Lim | Modeling single-file diffusion by step fractional brownian motion and generalized fractional langevin equation[END_REF], finance [START_REF] Bianchi | Efficient markets and behacioral finance : a comprehensive multifractal model[END_REF][START_REF] Xiao | Parameter identification for drift fractional brownian motions with application to the chinese stock markets[END_REF], health [START_REF] Khalfa | Heart rate regulation processed through wavelet analysis and change detection., Some case studies[END_REF]. . . Most of the previous examples concern detection of change on the mean of series derived from the original one, as the series of energy calculated by the wavelet analysis [START_REF] Khalfa | Heart rate regulation processed through wavelet analysis and change detection., Some case studies[END_REF] and the series of Hurst index [START_REF] Bianchi | Efficient markets and behacioral finance : a comprehensive multifractal model[END_REF][START_REF] Wanliss | Space storm as phase transition[END_REF][START_REF] Xiao | Parameter identification for drift fractional brownian motions with application to the chinese stock markets[END_REF]. However, in all those cases, we still detect change on the mean of the derived series, that is change on the mean value of the Hurst series [START_REF] Bianchi | Efficient markets and behacioral finance : a comprehensive multifractal model[END_REF][START_REF] Wanliss | Space storm as phase transition[END_REF][START_REF] Xiao | Parameter identification for drift fractional brownian motions with application to the chinese stock markets[END_REF] or change of the mean value of the wavelet transform for the series of energy [START_REF] Khalfa | Heart rate regulation processed through wavelet analysis and change detection., Some case studies[END_REF]. To sum up, change point detection on the mean is a relevant question in many applications. On the other hand, in statistics, the change point analysis field has been studied for more than forty years [START_REF] Csörgö | Strong approximations in probability and statistics[END_REF][START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF][START_REF] Brodsky | Nonparametric methods in change-point problems[END_REF] or [START_REF] Hušková | Change point analysis based on the empirical characteristic functions[END_REF][START_REF] Gombay | Monitoring parameter change in ar(p) time series models[END_REF][START_REF] Rigaill | Learning sparse penalties for changepoint detection using max margin interval regression[END_REF] for an updated overview. Depending on the method of data acquisition, we distinguish two kinds of change point detection :

• We observe the whole time series and we want to detect all the change point a posteriori or offline, see e.g. [START_REF] Bertrand | Off-line detection of multiple change points by the filtered derivative with p-value method[END_REF][START_REF] Boutoille | A hybrid fusion system applied to off-line detection and change-points estimation[END_REF].

• We observe the time series and we want to detect a change point as soon as possible. It is the online change point detection, see e.g [START_REF] Fearnhead | On-line inference for multiple change point problems[END_REF][START_REF] Bertrand | Detecting small shift on the mean by finite moving average[END_REF].

In this work, we only consider the 'a posteriori' detection which is called change point analysis in the statistical literature. We describe our framework with a toy model in section 1.

At the beginning of 21st century, the method used for this kind of problem was the Penalized Least Square Criterion. This algorithm is based on the minimisation of the contrast function when the number of change points is known [START_REF] Bai | Estimating and testing linear models with multiple structural changes[END_REF]. When the number of change point is unknown, many authors use the penalized version of the contrast function [START_REF] Lavielle | Least-squares estimation of an unknown number of shifts in a time series[END_REF]. From a computational point of view, the PLS method uses dynamic programming algorithms and requires matrix operations. Therefore, the time and memory complexity of PLS algorithm is of order O(n 2 ), where n denotes the size of the dataset. Due to the data deluge, the sizes of datasets become larger and larger, to the point where the computational complexity of this statistical method has become a challenge, see e.g [START_REF] Ji | Multi-scale internet traffic analysis using piecewise self-similar processes[END_REF] for internet traffic, [START_REF] Bai | Estimating and testing linear models with multiple structural changes[END_REF][START_REF] Bertrand | Modelling NASDAQ series by sparse multifractional brownian motion[END_REF][START_REF] Bianchi | Efficient markets and behacioral finance : a comprehensive multifractal model[END_REF][START_REF] Frezza | Modeling the time-changing dependence in stock markets[END_REF] for economics or High Frequency finance, [START_REF] Khalfa | Heart rate regulation processed through wavelet analysis and change detection., Some case studies[END_REF][START_REF] Azzaoui | Classifying heart-rate by change detection and wavelet methods for emergency physicians[END_REF] for heartbeat series and health. Among the different methods for a posteriori change detection, the use of a Filtered Derivative function has been introduced by [START_REF] Benveniste | Detection of abrupt changes in signals and dynamical systems: some statistical aspects, Analysis and optimization of systems[END_REF][START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. The advantage of the Filtered Derivative method is the time and memory complexity, both of order O(n) [START_REF] Bertrand | Off-line detection of multiple change points by the filtered derivative with p-value method[END_REF][START_REF] Messer | The shark function -asymptotoc behavior of the filetered derivative for point process in case of change points[END_REF][START_REF] Soh | High-dimensional change-point estimation: Combining filtering with convex optimization[END_REF][START_REF] Messer | A multiple filter test for the detection of rate changes in renewal processes with varying variance[END_REF]. On the other hand, Filtered Derivative method leads to many false discoveries of change points. Recently, a new method called Filtered Derivative with p-value (FDpV) has been introduced [START_REF] Bertrand | Off-line detection of multiple change points by the filtered derivative with p-value method[END_REF]. FDpV is a two-step procedure: the first step is based on the Filtered Derivative function and detects the potential change points. In the second step we calculate their p-value to eliminate the false alarms. In [START_REF] Messer | A multiple filter test for the detection of rate changes in renewal processes with varying variance[END_REF], the first step is still based on Filtered Derivative, but the second step consists on increasing the window size A in order to find the true positives.

Yet, the problem of the false discoveries with the Filtered Derivative function, in the first step, was not resolved, even if in the second step, the number of false alarms and non-detections drop substantially. Indeed, most of the false discoveries at the first step will be discarded during the second step by calculating their p-values. However, these calculations still increase the computational time. This problem led us to think of a way to minimise the number of false alarms and also the number of undetected change points in the first step. Furthermore, we investigate the impact of the false positives and the undetected change points on the Mean Integrated Square Error (MISE). The rest of this article is organized as follows. In Section 1, we describe the problem of change point analysis with a toy model and we give some comparison criterion. In Section 2, we recall the method of the Penalized Least Square and the Filtered derivative with p-Value in order to analyse the problem of change points. In section 3, we expand on the Filtered Derivative with p-Value method by providing a method to choose the extra parameters of Step 1 and we show which impact is more important on MISE. All the technical proofs are postponed in appendices.

Change point analysis

In this section, we describe the problem of the change point analysis in a toy model that will be used throughout the sequel of this work. Then, we give some comparison criterion.

Toy model

Let X = (X 1 , X 2 , . . . , X n ) be a series indexed by the time t = 1, 2, . . . , n. We assume that a segmentation τ = (τ 1 , . . . , τ K ) exists such that :

• X t is a family of independent identically distributed (iid) random variables for t ∈ (τ k , τ k+1 ],

• k = 0, . . . , K, where by convention τ 0 = 0 and τ K+1 = n.

The most simple model is X ∼ N (µ(•), σ 2 ) a sequence of independent standard Gaussian variables such that X t ∈ N (µ(t), 1), where N (µ, σ 2 ) denotes the Gaussian law with mean µ and variance σ 2 . The function of time t -→ µ(t) is piecewise constant that is to say µ(t) = µ k for all t ∈ (τ k , τ k+1 ], see eg. Fig. 3 and Fig. 4. To sum up, we have :

• a configuration of K change points τ = (τ 1 , . . . , τ K ) enlarged, by convention, by adding τ 0 = 0 and τ K+1 = n,

• associated to the configuration of mean values µ = (µ 0 , . . . , µ K ),

• X t ∈ N (µ k , σ), for t ∈ (τ k , τ k+1
] and for all k = 0, . . . , K.

• For notational convenience, we define the configuration of shifts δ = (δ 1 , . . . , δ K ) where

δ k = µ k -µ k-1 , for k = 1, . . . , K.
• The minimal distance between two consecutive change points is defined by

L 0 = inf{|τ k+1 -τ k |, f or k = 0, . . . , K}.
• The minimal absolute value of the shifts is

δ 0 = inf{|δ k |, k = 1, . . . , K}. (1.1)
Let us also recall the definition of the cumulative distribution function for standard Gaussian law

Φ(x) = 1 √ 2π x -∞ e -u 2 2 du and Ψ(x) = 1 -Φ(x). (1.2)

The Comparison Criterion

We have to estimate the configuration τ = (τ 1 , . . . , τ K ) and the values of the mean µ = (µ 0 , µ 1 , . . . , µ K ). We denote the corresponding estimates by τ = (τ 1 , . . . , τ K ) and μ = (μ 0 , μ1 , . . . , μ K ). Stress that in real life situations the number of change points is also unknown and is estimated by K. In this frame, the comparison criterion concerning the different methods for change point analysis are :

1. The quality of estimation. For one sample, this quality can be measured by :

• The number of estimated change points. More precisely, the absolute value of the difference between the number of estimated and the number of true change points and | K -K|.

• The accuracy of the estimation of the change point. It is the distance between the true change points and the estimated change points, as defined by

d 2 (τ, τ ) = K k=1 |τ k -τ j(k) | 2 , (1.3)
where τ j(k) , for each k = 1, . . . , K, denotes the potential change point which is the closer to the right change point τ k .

• The integrated square error (ISE). Actually, we can reformulate the problem as a problem of estimation of a noisy signal, see eg. [START_REF] Arlot | Segmentation of the mean of heteroscedastic data via crossvalidation[END_REF][START_REF] Birgé | Minimal penalties for gaussian model selection[END_REF]. The signal is

s(t) = K k=0 µ k × 1 (τ k ,τ k+1 ] (t)
where we have set by convention τ 0 = 0 and τ K+1 = n. Therefore the estimated signal is

s(t) = K k=0 µ k × 1 ( τ k , τ k+1 ] (t)
and the integrated square error (ISE) is defined by

ISE = n t=1 ( s(t) -s(t)) 2
2. The mean value of estimations: a result obtained for just one simulation can be hazardous. So, we have to do M simulations, with e.g. M = 1, 000. Then, we calculate the mean integrated square error (MISE) and the histogram of K with the percentage of the true changes, or the mean and standard deviation of the misestimation of the number of change point : ( K -K).

3. The time complexity and the memory complexity: it is the mean CPU (Central Processing Unit) time for estimating s and the amount of memory is used.

2 Some methods for change point analysis

Penalized Least Square method (PLS)

Set S K = {τ = (τ 1 , . . . , τ K ) such that card(τ ) = K}, where card(τ ) denotes the dimension of the change points configuration τ .

1 st case : the number of change points K is known

For each configuration of change τ ∈ S K , we can define

µ k = mean(X, [τ k + 1, τ k+1 ]) f or k = 0, . . . , K (2.2)
where mean(X, Box) denotes the mean of the family X t for the indices t ∈ Box as defined by (2.1). Next, we search the configuration of change points τ K ∈ S K which minimises the square error Q(τ ) defined by

Q(τ ) = K k=0 τ k+1 t=τ k +1 |X t -µ k | 2 . (2.3)
That is to say such that τ K = argmin

τ K ∈S K Q(τ ).
2 nd case : the number of change points K is unknown

We remark that minimising the function Q(τ ) with an unknown number of changes will lead to consider the trivial configuration of changes τ = (1, 2, . . . , n) as optimal. To avoid this drawback, we add a penalty term proportional to the length of the change point configuration. Eventually, we want to minimise

pen(K) = Q( τ K ) + β × K for K = 0, . . . , n.
The parameter β adjusts the trade-off between minimising the square error Q( τ K ) and minimising the dimension of the change point configuration card( τ K ) = K. Indeed, a large value of β would allow to detect only the most significant change points, while a low value of β produces a high number of changes, with many false detections. Thus, different choices of the penalty coefficient β are possible. According to the criterion of information AIC and the Schwarz criterion, [START_REF] Lavielle | Least-squares estimation of an unknown number of shifts in a time series[END_REF] suggests to use a positive β n which converges to 0 when the series size n converges to infinity. For this model, [START_REF] Yao | Estimating the number of change-points via schwarzs criterion. statistics and probability letters[END_REF] has proved the consistency of the Schwarz criterion, with

β n = 2σ 2 (ln n) n .
The same choice is proposed in [START_REF] Lavielle | Least-squares estimation of an unknown number of shifts in a time series[END_REF][START_REF] Lavielle | Detection of multiple change points in multivariate time series[END_REF]. In [START_REF] Birgé | Minimal penalties for gaussian model selection[END_REF], the proposed choice is

β n = σ 2 n × 2 + 5 × ln( n K )
where σ 2 is the variance assumed to be constant and known and n the size of the series. In Fig. 1 below, we have plotted the contrast function and the penalized contrast function [START_REF] Lavielle | Least-squares estimation of an unknown number of shifts in a time series[END_REF][START_REF] Lavielle | Detection of multiple change points in multivariate time series[END_REF]. We clearly see in the figure 2.1 that the penalized contrast is almost horizontal. Thus, the minimal value fluctuates largely depending on the choice of the parameter β.

Filtered Derivative with p-Value method (FDpV)

In this subsection, we describe the FDpV method which is based on two procedures. In the first procedure, we use the function of Filtered Derivative to select a set of potential change points, whereas in the second one, we calculate the p-value for each change point in order to keep only the true positive. Precisely, the method is defined as follows:

Step 1 : Filtered Derivative

The first step (FD selection) depends on two parameters: the window size A and the threshold C 1 .

Computation of the Filtered Derivative function :

The Filtered Derivative function is defined as the difference between the estimators of the mean computed in two sliding windows respectively to the right and to the left of the time t, both of size A, with the following formula :

F D(t, A) = µ(X, [t + 1, t + A]) -µ(X, [t -A + 1, t]), (2.4) f or A < t < n -A,
where mean(X, Box) denotes the mean of the family X t for the indices t ∈ Box as defined by (2.1). This method consists on filtering data by computing the estimators of the parameter µ before applying a discrete derivation. So this construction explains the name of the algorithm, so called Filtered Derivative method [START_REF] Benveniste | Detection of abrupt changes in signals and dynamical systems: some statistical aspects, Analysis and optimization of systems[END_REF][START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. Next, remark that quantities A × F D(t, A) can be inductively calculated by using

A × F D(t + 1, A) = A × F D(t, A) + X(t + 1 + A) -2X(t + 1) + X(t -A + 1) (2.5)
Thus, the computation of the whole function t -→ F D(t) for t ∈ [A, n -A] requires O(n) operations and the storage of n real numbers.

Determination of the potential change points

Let us point that the absolute value of the Filtered Derivative |F D| presents hats at the vicinity of the change points as seen on the figure 2 below. In the case of a true change point, an error of estimation on the location of the change exists. Thus, we have to cancel a small vicinity of size ε k around each point τ * k , see [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF][START_REF] Bertrand | Off-line detection of multiple change points by the filtered derivative with p-value method[END_REF]. Then, for each segment, we calculate an estimation of the mean

µ k := µ(X, [τ k + ε k , τ k+1 -ε k+1 ]), (2.6)
where µ(X, Box) is defined by (2.1), and, as in [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF],

ε k = 5 × σ δ k 2
where x denotes the ceiling function of the real number x.

Remark 2.1 From definition (1.1), δ 0 is the lower bound of the shifts, thus we can deduce the following upper bound

ε k ≤ ε 0 := 5 × σ δ 0 2 . (2.7)
The standard deviation σ can be empirically estimated. Next, we can use this bound in Formula (2.6) and set

µ k := µ(X, [τ k + ε 0 , τ k+1 -ε 0 ]). (2.8)
After that, we eliminate the false detections in order to keep (as possible) only the true change points. In [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF], we apply the following hypothesis testing, for all 1 ≤ k ≤ K :

(H 0,k ) : µ k = µ k+1 versus (H 1,k ) : µ k = µ k+1
where the terms µ k are defined by (2.8). By using this second single hypothesis test, we calculate the p-values p * 1 , . . . , p * K associated to each potential change point τ 1 , . . . , τ K .

p-value computation

We choose the statistic Student T. Indeed, under the null hypothesis, t k has a Student distribution of degree

d = N k + N k-1 -2 such that N k = (τ * 0 -ε 0 ) -(τ * k + ε 0 ) , where t * k = µ k -µ k-1 S 2 k-1 N k-1 + S 2 k N k , (2.9)
and the sample standard deviation is :

S k = 1 N k τ k+1 -ε 0 t=τ k +ε 0 X 2 t -µ 2 k , (2.10)
By construction, d > 2A -4ε 0 , thus for A > 30 the distribution of t * k is approximatively Gaussian and we can set

p k = 2 × 1 -St d (|t k |) 2 × 1 -Φ(|t k |) (2.11)
where St d denotes the cumulative distribution function of a Student law of degree d and Φ the cumulative distribution function of the zero mean standard Gaussian law as given by (1.2). In [START_REF] Bertrand | Off-line detection of multiple change points by the filtered derivative with p-value method[END_REF], we only keep the change points corresponding to a p-value smaller than a fixed threshold p * 2 . Consequently, Step 2 is much more selective and allows us to deduce an estimator of the piecewise constant map t -→ µ(t).

3 How to choose the extra-parameters of Step 1 (Filtered Derivative)?

All the change points methods depend on extra-parameters which must be well chosen. The filtered derivative method depends on two extra-parameters, namely the window size A and the threshold C 1 .

Criterion for choosing the extra-parameters of FD

In the case of "At Most one Change point " (AMOC), the usual criteria are :

-Error of type I which corresponds to the Probability of False Alarm (PFA).

-Error of type II which corresponds to the Probability of Non Detection (PND).

see eg. [START_REF] Csörgo | Limit theorem in change-point analysis[END_REF][START_REF] Kolmogorov | Probabilist-statistical methods of detecting spontaneously occuring effects[END_REF]. The PND is well suited for AMOC. But for detecting more than one change point, we have to impose that one detected change point has to be at the vicinity of each real change point. Following [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF], for each real change point τ k , we define the local PND as

P N D local (τ k ) = I P (B k ) ,
where

B k = ∀t ∈ [τ k -A, τ k + A], |F D(t, A)| < C 1 .
With these notations, we can define the global PND by

P N D global = I P K k=1 B k . (3.1)
Next, we can define the Number of Non Detection (N N D) in the vicinity of the right change points as

N N D(ω) = K k=1 1 B k (ω). (3.2)
Furthermore, in case of no change points, the probability of false alarm is defined as :

α(A, C 1 ) = I P τ (C 1 , A) ≤ n -A , where τ (C 1 , A) is the first hitting time of C 1 , that is τ (C 1 , A) := inf{t ≥ A ; |F D(t, A)| ≥ C 1 }. (3.3)
However, type I error is the probability of at least one false alarm and thus appears as a rough criterion see [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF]. Following [START_REF] Benjanimi | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF], the number of false alarm (NFA) is a more relevant criterion.

Definition 3.1 The number of false alarms is defined as follows iii) The impacts of false alarm and non detection on MISE are mainly described by the two following propositions (3.3 and 3.4). Stress that, the potential change points can be different from the true change points.

N F A = K -K + N N D. ( 3
Definition 3.2 (MISE) Let τ 1 , τ 2 be two change times, set

s(t) = K k=0 µ k × 1 (τ k ,τ k+1 ] (t)
the true signal, and denote by

s(t) = K k=0 µ k × 1 ( τ k , τ k+1 ] (t)
the estimated signal.

The mean integrated square error between the times τ 1 and τ 2 is then defined by

M ISE(τ 1 , τ 2 ) := E τ 2 t=τ 1 | s(t) -s(t)| 2 .
Proposition 3.3 (False alarm) Let τ 1 , τ 2 be two successive change points and τ 1 and τ 2 the potential change points such that τ1 = τ 1 + ε 1 , τ2 = τ 2 -ε 2 . Furthermore, we assume that

ε i ∞ ≤ M ε i , with i ∈ {1, 2} where
M ε i is a finite constant.
i) Without false alarm : Assume that τ1 , τ2 are the two successive potential change points obtained after Step 1. Then

M ISE(τ 1 , τ 2 ) = σ 2 + r 3.5 (3.5) whith |r 3.5 | ≤ M ε 1 σ 2 τ 1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 + . . . • • • + M ε 2 σ 2 τ 3 -τ2 + τ 3 -τ 2 τ 3 -τ2 2 (µ 1 -µ 2 ) 2 .
ii) With false alarm : Assume that τ1 , τ2 , and τ3 are the three successive potential change points obtained after Step 1 such as τ3 is the false alarm, and that

τ 1 < τ1 < τ3 < τ2 < τ 2 . Then M ISE(τ 1 , τ 2 ) = 2σ 2 + r 3.6 (3.6)
with :

|r 3.6 | ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 + . . . • • • + M ε 2 σ 2 τ2 -τ 4 + τ 2 -τ 4 τ2 -τ 4 2 (µ 1 -µ 3 ) 2 .
Proof. See Appendix A.

Proposition 3.4 (Non detection) Let τ 1 , τ 2 , and τ 3 be three successive right change points Furthermore, we assume that ε i ∞ ≤ M ε i , with i ∈ {1, 2, 3} where M ε i is a finite constant and the quantity ε i are precisely defined below.

i) Without non detection : Assume that τ1 , τ2 , τ3 are the three successive potential change points obtained after Step 1 such that τ1 =

τ 1 + ε 1 , τ2 = τ 2 -ε 2 , τ3 = τ 3 -ε 3 and τ 1 < τ1 < τ 2 < τ2 < τ3 < τ 3 . M ISE(τ 1 , τ 3 ) = σ 2 + r 3.7 (3.7)
with :

|r 3.7 | ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 + . . . • • • + M ε 2 σ 2 τ3 -τ2 + τ3 -τ 2 τ3 -τ2 2 (µ 1 -µ 2 ) 2 + . . . • • • + M ε 3 σ 2 τ 4 -τ3 + τ 4 -τ 3 τ 4 -τ3 2 (µ 1 -µ 4 ) 2 .
ii) With non detection : Assume that τ 2 is the undetected change point. For notational convenience, let us still denote by τ1 and τ3 the two successive potential change points obtained after Step 1 such that τ1 = τ

1 + ε 1 , τ3 = τ 3 -ε 3 and τ 1 < τ1 < τ 2 < τ3 < τ 3 . Then M ISE(τ 1 , τ 3 ) = 2σ 2 + r 3.8 (3.8)
with :

|r 3.8 | ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 + . . . • • • + M ε 3 σ 2 τ 4 -τ3 + τ 4 -τ 3 τ 4 -τ3 2 (µ 1 -µ 3 ) 2 .
Proof. See Appendix B.

In order to measure which kind of alarm (error) impacts more the MISE, we study the difference between the MISE in the case of the non detection and the MISE in the case of the false alarm. To begin with, we restricted ourselves to the particular case where r 1,i -→ 0 with i ∈ {1, 2, 3, 4}, which means that ε j = 0 with j ∈ {1, 2, 3} and the potential change points are the true change points. Moreover, we obtain the same results at the first order in the general case when we put ε i = 0. ii) Let τ1 , τ2 , and τ3 be the three successive potential change points obtained after Step 1 such that

τ 1 = τ1 < τ3 < τ2 = τ 2 . Then M ISE(τ 1 , τ 2 ) = 2σ 2 .
iii) The difference between the MISE in the case with false alarm and the case without false alarm is as follows :

∆M ISE F A := M ISE withF A -M ISE withoutF A (3.9) = σ 2
τ 1 = τ1 < τ 2 < τ3 = τ 3 . Then M ISE(τ 1 , τ 3 ) = σ 2 + (τ 2 -τ 1 )(τ 3 -τ 2 ) τ 3 -τ 1 (δµ) 2
with δµ = (µ 2 -µ 1 ).

iii) The difference between the MISE in the case with non detection and the case without non detection is as follows :

∆M ISE N D := M ISE withN D -M ISE withoutN D (3.10) = (τ 2 -τ 1 )(τ 3 -τ 2 ) τ 3 -τ 1 (δµ) 2 -σ 2
Proposition 3.7 Let us assume that Propositions 3.3 and 3.4 are satisfied. Then

• if δµ σ 2 > 2, then ∆M ISE N D > ∆M ISE F A
, namely the impact of the non detection is more important than the impact of the false alarm .

• On the other side, if δµ σ

2 < 8 τ 3 -τ 1
, then ∆M ISE N D < ∆M ISE F A , namely the impact of the false alarm is more important than the impact of the non detection.

Proof. We want to show that :

δµ σ 2 > 2 =⇒ ∆M ISE N D > ∆M ISE F A (3.11)
We have :

∆M ISE N D -∆M ISE F A = σ 2 (δµ) σ 2 (τ 2 -τ 1 )(τ 3 -τ 2 ) τ 3 -τ 1 - 2
Thus,

∆M ISE N D -∆M ISE F A > 0 =⇒ ∆M ISE N D > ∆M ISE F A =⇒ δµ σ > 2 τ 3 -τ 1 (τ 2 -τ 1 )(τ 3 -τ 2 ) Set L = τ 3 -τ 1 , λ = τ 2 -τ 1 where L > 0 and λ ∈ [1, L -1] Thus, we have : f (λ) = 2L λ(L-λ)
The minimum of f (λ) is reached for λ = L/2, and the maximum is reached for λ = 1. We deduce then :

8 L ≤ f (λ) ≤ 2L L -1 < √ 2 Therefore, δµ σ > √ 2 implies ∆M ISE N D > ∆M ISE F A .
On the other hand, when δµ σ > 8 L we have ∆M ISE N D < ∆M ISE F A .

Bound on the error of type II for Filtered Derivative

Proposition 3.8 Let τ = (τ 1 , τ 2 , . . . , τ K ) be a configuration of K change points, with means µ = (µ 0 , . . . , µ K ) and shifts δ = (δ 1 , . . . , δ K ) as described in Subsection 1.1. Then

P N D global ≤ K × β * (C 1 , A),
where P N D global is defined by (3.1),

β * (C 1 , A) := Ψ δ 0 -C 1 σ A 2 × Φ C 1 -δ 0 /3 σ A 2 2 , (3.12)
δ 0 is defined by (1.1), and Φ and Ψ are given by (1.2).

Proof. Actually, the proposition Prop.

B k = ∀t ∈ [τ k -A, τ k + A], F D(t, A) < C 1 .
Thus, I P B k ≤ I P B k . On the other hand, following [7, Prop. 3.2, p 222], we have,

I P (B k ) ≤ Ψ |δ k | -C 1 σ A 2 × Φ C 1 -|δ k |/3 σ A 2 2 . ( 3
I P (B k ) ≤ β * (C 1 , A) := Ψ δ 0 -C 1 σ A 2 × Φ C 1 -δ 0 /3 σ A 2 2 . (3.14)
On the other hand, we obviously have

P N D global ≤ K k=1 I P (B k )
which combined with (3.14) provides the bound (3.8). This finishes the proof of Proposition 3.8.

Control of the number of false alarms

In this subsection, we want to control not only the Probability of False Alarm (PFA) but also the Number of False Alarms (NFA) (see Definition 3.1). Let us denote by K the number of change point select in step 1 of the FDpV method (FD), then the number of false alarms is ( K -K). Moreover, in order to control the number of false alarms, we need to choose the extraparameters A and C 1 .

The choice of parameter A

From the subsection 3.3, we can get the feeling that the larger the window size A is, the smaller type I and type II errors will be. This reasoning holds true as long as

2 × A < L 0 := inf{|τ k+1 -τ k |, k = 1, . . . , K}. (3.15)
Thus, we have to choose a parameter A < L 0 /2, even if we do not exactly know the quantity L 0 .

The choice of parameter C 1

In [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF] [START_REF] Bertrand | A local method for estimating change points: the "hat-function[END_REF], we have C 1 < δ 0 with δ 0 = inf{|δ k |, k = 1, ..., K} where δ k is the size of the average of µ k . With different Monte-Carlo simulations (see the subsection 3.5), we note that the best values of C 1 are between 0.1 and 0.2.

Monte-Carlo simulation

This Monte-Carlo simulation is done for M = 1000. Let (X j 1 , X j 2 , . . . , X j n ) be a sequence of simulated Gaussian random variables where n = 5000 or n = 50, 000 and j = 1, . . . , M , with a variance σ 2 = 1 and a mean µ t = f (t) where f is a piecewise constant function with a specified number of change points at different times τ with different means µ. On each sample, we apply the FDpV method for different values of the extra-parameters A and C 1 . We vary the parameter A between 20 and 220 with a step of 10 while we vary the parameter C 1 between 0.1 and 1 with a step of 0.05.

Example 3

In this example, we keep the variables n and τ of the example 2 and we take shifts on the mean µ of the change point smaller than those in the example 1 and 2, that is to say δµ = (-2, -1.

Conclusion

Our analysis suggest a comprehensive way to optimize the Filtered Derivative function extraparameters in the first step of the FDpV method. By giving the precise values to choose for the window size A and the threshold C 1 where we have to choose A > 70 and C 1 ∈ [0.15, 0.22], we obtain fewer false alarms and undetected change points. Thus, in the second step, we calculate less p-values, gaining computational time and memory. Another method is used in [START_REF] Messer | A multiple filter test for the detection of rate changes in renewal processes with varying variance[END_REF]. The first step is still based on the Filtered Derivative function with the two extra-parameters: the window size A and the threshold C 1 . Yet, in the second step, they increase the window size A such that for each value of A they detect a real change point. This algorithm may detect all the true change points. But the problem of the false discoveries is not resolved, and can be enhanced by a good initial choice of the parameters A and C 1 . On the other hand, by fixing a value for the window size A and the threshold C 1 in Step 1 of the FDpV method, we detect the real change points and, at the same time, we minimise the number of the false discoveries. Furthermore, in Step 2, we eliminate more false positives. Finally, we studied the impact of the false alarms and the undetected change points on the Mean Integrated Square Error in order to determine which one has a more significant impact. To sum up, for the FDpV method, we have specified the values of the window size A and the threshold C 1 in Step 1. In forthcoming studies, we need to find a way to estimate the value of the threshold p * 2 in Step 2 to discard as many false positives as possible.

A Proof of Proposition 3.3

Proof.

i) Without false alarm :

ISE(τ 1 , τ 2 ) = ISE(τ 1 , τ1 ) + ISE(τ 1 , τ2 ) + ISE(τ 2 , τ 2 ) -ISE(τ 1 , τ1 ) = ε 1 (µ 1 -μ0 ) 2
For τ 0 = τ0 et τ 3 = τ3 , we have :

μ0 = σ √ τ1 -τ 0 × U 0 + (τ 1 -τ 1 )µ 1 +(τ 1 -τ 0 )µ 0 τ1 -τ 0 µ 1 -μ0 = µ 1 - σ √ τ1 -τ 0 × U 0 - (τ 1 -τ 1 )µ 1 -(τ 1 -τ 0 )µ 0 τ1 -τ 0 = (τ 1 -τ 0 )µ 1 -(τ 1 -τ 1 )µ 1 -(τ 1 -τ 0 )µ 0 τ1 -τ 0 - σ √ τ1 -τ 0 × U 0 = τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 M ISE(τ 1 , τ1 ) = E ε 1 (µ 1 -μ0 ) 2 = E ε 1 τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ≤ ε 1 ∞ × E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ε 1 ∞ = sup|ε 1 k | ≤ M ε 1 E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) -σ √ τ1 -τ 0 × U 0 2 = σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2
Thus :

M ISE(τ 1 , τ1 ) ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 ISE(τ 1 , τ2 ) = (τ 2 -τ1 )(µ 1 -μ1 ) 2
We have :

μ1 = µ 1 + σ √ τ2 -τ1 × U 1 ⇒ µ 1 -μ1 = -σ √ τ2 -τ1 × U 1 Thus : M ISE(τ 1 , τ2 ) = E (τ 2 -τ1 )(µ 1 -μ1 ) 2 = (τ 2 -τ1 )E(µ 1 -μ1 2 ) = (τ 2 -τ1 ) σ 2 τ2 -τ1 × E(U 2 1 ) = σ 2 ISE( τ2 -τ 2 ) = ε 2 (µ 1 -μ2 ) 2 we have : μ2 = σ √ τ3 -τ2 × U 2 + (τ 2 -τ 2 )µ 1 +(τ 3 -τ 2 )µ 2 τ3 -τ 2 µ 1 -μ2 = µ 1 - σ √ τ 3 -τ2 × U 2 - (τ 2 -τ2 )µ 1 -(τ 3 -τ 2 )µ 2 τ 3 -τ2 = (τ 3 -τ 2 )µ 1 -(τ 2 -τ2 )µ 1 -(τ 3 -τ 2 )µ 2 τ 3 -τ2 - σ √ τ 3 -τ2 × U 2 = τ 3 -τ 2 τ 3 -τ2 (µ 1 -µ 2 ) - σ √ τ 3 -τ2 × U 2 M ISE( τ2 , τ 2 ) = E ε 2 (µ 1 -μ2 ) 2 = E ε 2 τ 3 -τ 2 τ 3 -τ2 (µ 1 -µ 2 ) - σ √ τ 3 -τ2 × U 2 2 ≤ ε 2 ∞ × E τ 3 -τ 2 τ 3 -τ2 (µ 1 -µ 2 ) - σ √ τ 3 -τ2 × U 2 2 ε 2 ∞ = sup|ε 2 k | ≤ M ε 2 E τ 3 -τ 2 τ 3 -τ 2 (µ 1 -µ 2 ) -σ √ τ 3 -τ 2 × U 2 2 = σ 2 τ 3 -τ 2 + τ 3 -τ 2 τ 3 -τ 2 2 (µ 1 -µ 2 ) 2
Thus :

M ISE(τ 1 , τ1 ) ≤ M ε 2 σ 2 τ 3 -τ2 + τ 3 -τ 2 τ 3 -τ2 2 (µ 1 -µ 2 ) 2
Thus, we can deduce that :

M ISE(τ 1 , τ 2 ) ≤ σ 2 + ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 (A.1) + ε 2 σ 2 τ 3 -τ2 + τ 3 -τ 2 τ 3 -τ2 2 (µ 1 -µ 2 ) 2
ii) With false alarm :

ISE(τ 1 , τ 2 ) = ISE(τ 1 , τ1 ) + ISE(τ 1 , τ3 ) + ISE(τ 3 , τ2 ) + ISE(τ 2 , τ 2 ) ISE(τ 1 , τ1 ) = ε 1 (µ 1 -μ0 ) 2 we have : μ0 = σ √ τ1 -τ 0 × U 0 + (τ 1 -τ 1 )µ 1 +(τ 1 -τ 0 )µ 0 τ1 -τ 0 µ 1 -μ0 = µ 1 - σ √ τ1 -τ 0 × U 0 - (τ 1 -τ 1 )µ 1 + (τ 1 -τ 0 )µ 0 τ1 -τ 0 = τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 M ISE(τ 1 , τ1 ) = E ε 1 (µ 1 -μ0 ) 2 = E ε 1 τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ≤ ε 1 ∞ × E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ε 1 ∞ = sup|ε 1 k | ≤ M ε 1 E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) -σ √ τ1 -τ 0 × U 0 2 = σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2
Thus :

M ISE(τ 1 , τ1 ) ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 ISE(τ 1 , τ3 ) = (τ 3 -τ1 )(µ 1 -μ1 ) 2
We have :

μ1 = µ 1 + σ √ τ3 -τ1 × U 1 ⇒ µ 1 -μ1 = -σ √ τ3 -τ1 × U 1 Thus : M ISE(τ 1 , τ3 ) = E (τ 3 -τ1 )(µ 1 -μ1 ) 2 = (τ 3 -τ1 )E(µ 1 -μ1 ) 2 = (τ 3 -τ1 ) σ 2 τ3 -τ1 × E(U 2 1 ) = σ 2 ISE(τ 3 , τ2 ) = (τ 2 -τ3 )(µ 1 -μ2 2 ) We have : μ1 = µ 2 + σ √ τ2 -τ 3 × U 2 ⇒ µ 1 -μ2 = -σ √ τ2 -τ 3 × U 2 Thus : M ISE(τ 3 , τ2 ) = E (τ 2 -τ3 )(µ 1 -μ1 ) 2 = (τ 2 -τ3 )E(µ 1 -μ2 ) 2 = (τ 2 -τ3 ) σ 2 τ2 -τ3 × E(U 2 2 ) = σ 2 ISE(τ 2 , τ 2 ) = ε 2 (µ 1 -μ3 ) 2 we have : μ3 = σ √ τ2 -τ 4 × U 2 + (τ 2 -τ 2 )µ 1 +(τ 2 -τ 4 )µ 3 τ2 -τ 4 µ 1 -μ3 = µ 1 - σ √ τ2 -τ 4 × U 2 - (τ 2 -τ 2 )µ 1 + (τ 2 -τ 4 )µ 3 τ2 -τ 4 = τ 2 -τ 4 τ2 -τ 4 (µ 1 -µ 3 ) - σ √ τ2 -τ 4 × U 2 M ISE(τ 2 , τ 2 ) = E ε 2 (µ 1 -μ3 ) 2 = E ε 2 τ 2 -τ 4 τ2 -τ 4 (µ 1 -µ 3 ) - σ √ τ2 -τ 4 × U 2 2 ≤ ε 2 ∞ × E τ 2 -τ 4 τ2 -τ 4 (µ 1 -µ 3 ) - σ √ τ2 -τ 4 × U 2 2 ε 2 ∞ = sup|ε 2 k | ≤ M ε 2 E τ 2 -τ 4 τ2 -τ 4 (µ 1 -µ 3 ) -σ √ τ2 -τ 4 × U 2 2 = σ 2 τ2 -τ 4 + τ 2 -τ 4 τ2 -τ 4 2 (µ 1 -µ 3 ) 2
Thus :

M ISE(τ 2 , τ2 ) ≤ M ε 2 σ 2 τ2 -τ 4 + τ 2 -τ 4 τ2 -τ 4 2 (µ 1 -µ 3 ) 2
Then :

M ISE(τ 1 , τ 2 ) ≤ 2σ 2 + M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 (A.2) + M ε 2 σ 2 τ2 -τ 4 + τ 2 -τ 4 τ2 -τ 4 2 (µ 1 -µ 3 ) 2 B Proof of Proposition 3.4 Proof. i) Without non detection : ISE(τ 1 , τ 3 ) = ISE(τ 1 , τ1 ) + ISE(τ 1 , τ2 ) + ISE(τ 2 , τ 2 ) + ISE(τ 2 , τ3 ) + ISE(τ 3 , τ 3 ) -ISE(τ 1 , τ1 ) = ε 1 (µ 1 -μ0 ) 2
For τ 0 = τ0 et τ 4 = τ4 , we have :

μ0 = σ √ τ1 -τ 0 × U 0 + (τ 1 -τ 1 )µ 1 +(τ 1 -τ 0 )µ 0 τ1 -τ 0 µ 1 -μ0 = µ 1 - σ √ τ1 -τ 0 × U 0 - (τ 1 -τ 1 )µ 1 -(τ 1 -τ 0 )µ 0 τ1 -τ 0 = (τ 1 -τ 0 )µ 1 -(τ 1 -τ 1 )µ 1 -(τ 1 -τ 0 )µ 0 τ1 -τ 0 - σ √ τ1 -τ 0 × U 0 = τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 M ISE(τ 1 , τ1 ) = E ε 1 (µ 1 -μ0 ) 2 = E ε 1 τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ≤ ε 1 ∞ × E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ε 1 ∞ = sup|ε 1 k | ≤ M ε 1 E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) -σ √ τ1 -τ 0 × U 0 2 = σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2
Thus :

M ISE(τ 1 , τ1 ) ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 ISE(τ 1 , τ2 ) = (τ 1 , τ2 )(µ 1 -μ1 ) 2
We have :

μ1 = µ 1 + σ √ τ2 -τ1 × U 1 ⇒ µ 1 -μ1 = -σ √ τ2 -τ1 × U 1 Thus : M ISE(τ 1 , τ2 ) = E (τ 2 -τ1 )(µ 1 -μ1 ) 2 = (τ 2 -τ1 )E(µ 1 -μ1 ) 2 = (τ 2 -τ1 ) σ 2 τ2 -τ1 × E(U 2 1 ) = σ 2 ISE( τ2 -τ 2 ) = ε 2 (µ 1 -μ2 ) 2 we have : μ2 = σ √ τ3 -τ2 × U 2 + (τ 2 -τ 2 )µ 1 +(τ 3 -τ 2 )µ 2 τ3 -τ 2 µ 1 -μ2 = µ 1 - σ √ τ 3 -τ2 × U 2 - (τ 2 -τ2 )µ 1 -(τ 3 -τ 2 )µ 2 τ 3 -τ2 = (τ 3 -τ 2 )µ 1 -(τ 2 -τ2 )µ 1 -(τ 3 -τ 2 )µ 2 τ 3 -τ2 - σ √ τ 3 -τ2 × U 2 = τ 3 -τ 2 τ 3 -τ2 (µ 1 -µ 2 ) - σ √ τ 3 -τ2 × U 2 M ISE( τ2 , τ 2 ) = E ε 2 (µ 1 -μ2 ) 2 = E ε 2 τ 3 -τ 2 τ 3 -τ2 (µ 1 -µ 2 ) - σ √ τ 3 -τ2 × U 2 2 ≤ ε 2 ∞ × E τ 3 -τ 2 τ 3 -τ2 (µ 1 -µ 2 ) - σ √ τ 3 -τ2 × U 2 ε 2 ∞ = sup|ε 2 k | ≤ M ε 2 E τ 3 -τ 2 τ 3 -τ 2 (µ 1 -µ 2 ) -σ √ τ 3 -τ 2 × U 2 2 = σ 2 τ 3 -τ 2 + τ 3 -τ 2 τ 3 -τ 2 2 (µ 1 -µ 2 ) 2
Thus :

M ISE(τ 1 , τ1 ) ≤ ε 2 σ 2 τ 3 -τ2 + τ 3 -τ 2 τ 3 -τ2 2 (µ 1 -µ 2 ) 2 ISE(τ 2 , τ3 ) = (τ 3 -τ 2 )(µ 1 -μ2 ) 2
We have :

μ2 = µ 1 + σ √ τ3 -τ 2 × U 3 ⇒ µ 1 -μ2 = -σ √ τ2 -τ1 × U 3 Thus : M ISE(τ 1 , τ2 ) = E (τ 2 -τ1 )(µ 1 -μ2 ) 2 = (τ 2 -τ1 )E(µ 1 -μ1 2 ) = (τ 2 -τ1 ) σ 2 τ2 -τ1 × E(U 2 3 ) = σ 2 ISE( τ3 , τ 3 ) = ε 3 (µ 1 -μ4 ) 2
we have : μ4 = σ √ τ 4 -τ3 × U 4 + (τ 3 -τ 3 )µ 1 +(τ 4 -τ 3 )µ 4 τ 4 -τ3

µ 1 -μ4 = µ 1 - σ √ τ 4 -τ3 × U 4 - (τ 3 -τ3 )µ 1 + (τ 4 -τ 3 )µ 4 τ 4 -τ 3 = τ 4 -τ 3 τ 4 -τ3 (µ 1 -µ 4 ) - σ √ τ 4 -τ3 × U 4 M ISE( τ3 , τ 3 ) = E ε 3 (µ 1 -μ4 ) 2 = E ε 3 τ 4 -τ 3 τ 4 -τ3 (µ 1 -µ 4 ) - σ √ τ 4 -τ3 × U 4 2 ≤ ε 3 ∞ × E τ 4 -τ 3 τ 4 -τ3 (µ 1 -µ 4 ) - σ √ τ 4 -τ3 × U 4 2 ε 3 ∞ = sup|ε 3 k | ≤ M ε 3 E τ 4 -τ 3 τ 4 -τ 3 (µ 1 -µ 4 ) -σ √ τ 4 -τ 3 × U 4 2 = σ 2 τ 4 -τ 3 + τ 4 -τ 3 τ 4 -τ 3 2 (µ 1 -µ 4 ) 2
Thus :

M ISE(τ 3 , τ 3 ) ≤ M ε 3 σ 2 τ 4 -τ3 + τ 4 -τ 3 τ 4 -τ3 2 (µ 1 -µ 4 ) 2
Thus, we can deduce that :

M ISE(τ 1 , τ 3 ) ≤ 2σ 2 + M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 (B.1) + M ε 2 σ 2 τ 3 -τ2 + τ 3 -τ 2 τ 3 -τ2 2 (µ 1 -µ 2 ) 2 + M ε 3 σ 2 τ 4 -τ3 + τ 4 -τ 3 τ 4 -τ3 2 (µ 1 -µ 4 ) 2
ii) With non detection : ISE(τ 1 , τ 3 ) = ISE(τ 1 , τ1 )+ISE(τ 1 , τ 2 )+ISE(τ 2 , τ3 )+ISE(τ 3 , τ 3 ) ISE(τ 1 , τ1 ) = ε 1 (µ 1 -μ0 ) 2 For τ 0 = τ0 et τ 4 = τ4 , we have :

μ0 = σ √ τ1 -τ 0 × U 0 + (τ 1 -τ 1 )µ 1 +(τ 1 -τ 0 )µ 0 τ1 -τ 0 µ 1 -μ0 = µ 1 - σ √ τ1 -τ 0 × U 0 - (τ 1 -τ 1 )µ 1 -(τ 1 -τ 0 )µ 0 τ1 -τ 0 = (τ 1 -τ 0 )µ 1 -(τ 1 -τ 1 )µ 1 -(τ 1 -τ 0 )µ 0 τ1 -τ 0 - σ √ τ1 -τ 0 × U 0 = τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 M ISE(τ 1 , τ1 ) = E ε 1 (µ 1 -μ0 ) 2 = E ε 1 τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ≤ ε 1 ∞ × E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) - σ √ τ1 -τ 0 × U 0 2 ε 1 ∞ = sup|ε 1 k | ≤ M ε 1 E τ 1 -τ 0 τ1 -τ 0 (µ 1 -µ 0 ) -σ √ τ1 -τ 0 × U 0 2 = σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2
Thus :

M ISE(τ 1 , τ1 ) ≤ M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 ISE(τ 1 , τ2 ) = (τ 1 , τ2 )(µ 1 -μ1 ) 2
We have :

μ1 = µ 1 + σ √ τ2 -τ1 × U 1 ⇒ µ 1 -μ1 = -σ √ τ2 -τ1 × U 1 Thus : M ISE(τ 1 , τ 2 ) = E (τ 2 -τ1 )(µ 1 -μ1 ) 2 = (τ 2 -τ1 )E(µ 1 -μ1 ) 2 = (τ 2 -τ1 ) σ 2 τ 2 -τ1 × E(U 2 1 )
= σ 2 ISE(τ 2 , τ3 ) = (τ 3 -τ 2 )(µ 1 -μ2 ) 2 We have :

μ2 = µ 1 + σ √ τ3 -τ 2 × U 2 ⇒ µ 1 -μ2 = -σ √ τ2 -τ1 × U 2 Thus : M ISE(τ 1 , τ2 ) = E (τ 2 -τ1 )(µ 1 -μ2 ) 2 = (τ 2 -τ1 )E(µ 1 -μ1 2 ) = (τ 2 -τ1 ) σ 2 τ2 -τ1 × E(U 2 2 ) = σ 2 ISE( τ3 , τ 3 ) = ε 3 (µ 1 -μ3 ) 2
we have : μ3 = σ √ τ 4 -τ3 × U 3 + (τ 3 -τ 3 )µ 1 +(τ 4 -τ 3 )µ 3 τ 4 -τ3

µ 1 -μ3 = µ 1 - σ √ τ 4 -τ3 × U 3 - (τ 3 -τ3 )µ 1 + (τ 4 -τ 3 )µ 3 τ 4 -τ3 = τ 4 -τ 3 τ 4 -τ3 (µ 1 -µ 3 ) - σ √ τ 4 -τ3 × U 3 M ISE( τ3 , τ 3 ) = E ε 3 (µ 1 -μ3 ) 2 = E ε 3 τ 4 -τ 3 τ 4 -τ3 (µ 1 -µ 3 ) - σ √ τ 4 -τ3 × U 3 2 ≤ ε 3 ∞ × E τ 4 -τ 3 τ 4 -τ3 (µ 1 -µ 3 ) - σ √ τ 4 -τ3 × U 3 2 ε 3 ∞ = sup|ε 3 k | ≤ M ε 3 E τ 4 -τ 3 τ 4 -τ 3 (µ 1 -µ 3 ) -σ √ τ 4 -τ 3 × U 3 2 = σ 2 τ 4 -τ 3 + τ 4 -τ 3 τ 4 -τ 3 2 (µ 1 -µ 3 ) 2
Thus :

M ISE(τ 1 , τ 3 ) ≤ M ε 3 σ 2 τ 4 -τ3 + τ 4 -τ 3 τ 4 -τ3 2 (µ 1 -µ 3 ) 2
Thus, we can deduce that :

M ISE(τ 1 , τ 3 ) ≤ 2σ 2 + M ε 1 σ 2 τ1 -τ 0 + τ 1 -τ 0 τ1 -τ 0 2 (µ 1 -µ 0 ) 2 (B.
2)

+ M ε 3 σ 2 τ 4 -τ3 + τ 4 -τ 3 τ 4 -τ3 2 (µ 1 -µ 3 ) 2

Figure 1 :

 1 Figure 1: Blue with red crosses : the contrast function Q( τ K ); green : the penalized contrasted pen(K).

Figure 2 :

 2 Figure 2: A graph of the Filtered Derivative function

Figure 3 :

 3 Figure 3: Filtered Derivative function without noise (σ = 0).

Figure 4 :

 4 Figure 4: Filtered Derivative function with noise (σ = 1). The blue circles correspond to false alarms.

. 4 ) 3 . 2

 432 Impact of non detection and false alarm on the MISE i) From (3.4), we get | K -K| = |N F A -N N D|. Thus non detection and false alarm impact in the same way the criterion | K -K|.ii) Clearly, false alarm does not impact the quantity d 2 (τ, τ ) defined by (1.3), whereas non detection at the right place increases the quantity d 2 (τ, τ ).

Corollary 3 . 5 (

 35 False alarm) i) Let τ1 , τ2 be the two successive potential change points obtained after Step 1. Then M ISE(τ 1 , τ 2 ) = σ 2

Corollary 3 . 6 (

 36 Non detection) i) Let τ1 , τ2 and τ3 be the three successive potential change points obtained after Step 1 such that τ1 = τ 1 , τ2 = τ 2 and τ3 = τ 3 . Then M ISE(τ 1 , τ 3 ) = 2σ 2 ii) Next, assume that τ 2 is the undetected change point. For mathematical convenience, let us still denote by τ1 and τ3 the two successive potential change points obtained after Step 1 such that

. 13 )

 13 Next, by remarking that the right side of (3.13) is a decreasing function of |δ k | and recalling that δ 0 = inf k=1,...,K |δ k |, we can deduce that

Example 1

 1 The figures 5 and 6 show respectively the variation of the mean number of non detected change points and the variation of the mean number of the false alarms in function of the extra-parameters A and C 1 where n = 5000 ; τ = (1000, 1250, 1500, 2000, 3500, 4000, 4500) and δµ = (-0.5, -1, -1.5, 0.5, 1, 1.5) with δ 0 = 0.5.

Figure 5 :Figure 6 :

 56 Figure 5: The mean of the number of undetected change points

Figure 7 : 2

 72 Figure 7: The contours of the figures 5 and 6 and their overlap

Figure 8 :

 8 Figure 8: The mean of the number of undetected change points

Figure 9 :

 9 Figure 9: The mean of the number of false alarm

Figure 10 :

 10 Figure 10: The contours of the figures 8 and 9 and their overlap

  5, -1, -0.25, 0.25, 1, 1.5, 2) with δ 0 = 0.25. The figures 11 and 12 show the variation of the mean number of the undetected change points and those of the mean number of the false alarms for different values of the extra-parameters A and C 1 .

Figure 11 :Figure 12 :

 1112 Figure 11: The mean of the number of undetected change points

Figure 13 :

 13 Figure 13: The contours of the figures 11 and 12

Research supported by grant ANR-12-BS01-0016-01 entitled "Do Well B."