Hédi Hamdane Ifsttar-Lma

Thierry Serre
email: thierry.serre@ifsttar.fr

Catherine Masson
email: catherine.masson@ifsttar.fr

Robert Anderson
email: robert@casr.adelaide.edu.au

Issues and challenges for pedestrian active safety systems based on real world accidents 1 Issues and challenges for pedestrian active safety systems based on real world accidents

Keywords: AEB systems, Pedestrian detection, Safety performance, Accident reconstruction

The purpose of this study was to analyze real crashes involving pedestrians in order to evaluate the potential effectiveness of Autonomous Emergency Braking systems (AEB) in pedestrian protection. A sample of 100 real accident cases were reconstructed providing a comprehensive set of data describing the interaction between the vehicle, the environment and the pedestrian all along the scenario of the accident. A generic AEB system based on a camera sensor for pedestrian detection was modeled in order to identify the functionality of its different attributes in the timeline of each crash scenario. These attributes were assessed to determine their impact on pedestrian safety. The influence of the detection and the activation of the AEB system were explored by varying the Field Of View (FOV) of the sensor and the level of deceleration. A FOV of 35° was estimated to be required to detect and react to the majority of crash scenarios. For the reaction of a system (from hazard detection to triggering the brakes), between 0,5 and 1s appears necessary.

Introduction

Every year, more than 1.2 million people in the world die in traffic accidents and among these casualties 22% are pedestrians [START_REF] Who | WHO | Global status report on road safety[END_REF]. This proportion of fatalities that are pedestrians is 14% in Australia and 16% in the EU-15 (OECD/ITF, 2013). Worldwide organizations such as the OECD (Organization for Economic Co-operation and Development) and the WHO (World Health Organization) have outlined a set of objectives and actions to enhance pedestrian safety. Among these measures, the development of new safety-based technologies in the vehicle is promoted. Passive safety measures such as energy absorbing bumpers and hoods are designed to reduce the injury outcomes for pedestrians. These devices are evaluated under regulatory and non-regulatory tests (EuroNCAP). Along with these passive systems, active safety systems are being developed and introduced to prevent crashes. These active systems use sensors to monitor the forward path of the vehicle in order to detect a pedestrian. Once a hazard is detected, these systems trigger various countermeasures to avoid or mitigate collisions. These measures may include autonomous emergency braking or autonomous steering (e.g. [START_REF] Broggi | A New Approach to Urban Pedestrian Detection for Automatic Braking[END_REF][START_REF] Hayashi | Autonomous collision avoidance system by combined control of steering and braking using geometrically optimised vehicular trajectory[END_REF]Keller et al., 2011a).

Active safety systems are mainly composed of the three following components: sensors for detection, a unit for processing and actuators for triggering an emergency maneuver. Concerning the first component, in order to detect various obstacles, cameras operating in visible light or infrared radiation (Near, Mid, Far) as well as RADARs and Laser Scanners are used. These different sensors have complementary functions as described by [START_REF] Gandhi | Pedestrian Protection Systems: Issues, Survey, and Challenges[END_REF]. So the combination of multiple sensors allows mare accurate detection. This is achieved by merging and filtering the data collected from the environment in order to distinguish pedestrians from other background obstacles. As soon as pedestrians are detected, they are tracked in order to predict any collision. If a crash is imminent, the system applies the appropriate countermeasure.

Several methods have been developed to assess the effectiveness of these systems and estimate their safety impact on real world crashes. Approaches based on numerical simulation have been explored to assess the effect of systems in various accident scenarios. These scenarios are provided from in-depth crash investigations and are simulated using simple models [START_REF] Rosén | Pedestrian injury mitigation by autonomous braking[END_REF][START_REF] Seiniger | An open simulation approach to identify chances and limitations for vulnerable road user (VRU) active safety[END_REF]. More complex methods (probabilistic methods as Monte Carlo) expand the scenarios by slightly varying the initial pre-crash conditions in order to cover a wider range of crash configurations [START_REF] Lindman | A Method for Estimating the Benefit of Autonomous Braking Systems Using Traffic Accident Data[END_REF]. Besides computational simulations, researchers are currently attempting to develop standardized on-track tests to assess the performance of these systems. However, these track tests are limited in numbers of test scenarios even if they are focusing on reference scenarios defined through cluster and statistical analyses of real world crashes [START_REF] Ebner | Identifying and Analyzing Reference Scenarios for the Development and Evaluation of Active Safety: Application to Preventive Pedestrian Safety[END_REF][START_REF] Eckert | An integrated ADAS solution for pedestrian collision avoidance[END_REF][START_REF] Lenard | Typical pedestrian accident scenarios for the development of autonomous emergency braking test protocols[END_REF][START_REF] Wisch | European project AsPeCSS -interim result: Development of test scenarios based on identified accident scenarios[END_REF].

Several active systems for pedestrian safety have already been introduced on the market such as the CWAB-PD (Collision Warning with Full Auto Brake and Pedestrian Detection) firstly used in the Volvo S60 MY2011 and the Eyesight system available in Subaru's 2013 models. The regulations and protocols are being developed [START_REF] Schram | Implementation of autonomous emergency braking (AEB), the next step in Euro NCAP's safety assessment[END_REF] and some studies on the effectiveness of these active systems have been already presented. For example, [START_REF] Lindman | Benefit estimation model for pedestrian auto brake functionality[END_REF] used a simulation based method to conclude that the CWAB-PD system may lead to a reduction of about 30% of all front-end accidents and 24% of pedestrian fatalities. [START_REF] Ando | An evaluation protocol for collision avoidance and mitigation systems and its application to safety estimation[END_REF] assessed two systems using "standard tests" (vehicle driving towards a pedestrian dummy positioned on the course). The results attest that the functionality of the systems depends on the vehicle speed to avoid collisions but it is limited at a certain speed (40 km/h). Although these systems are intended to effectively reduce pedestrian injury outcomes, the collision avoidance performance of these systems remains limited.

The objective of this study is to understand the functionality of these active safety systems to prevent pedestrian crashes in terms of detecting the hazard and triggering Autonomous Emergency Braking (AEB). These issues explored in this paper are focused on the implication of using sensors with a wide field of view to detect pedestrians and the compromise between high effectiveness and a low false activation rate (C. G. Keller et al., 2011). This paper will address the problem of defining criteria for AEB activation in terms of:

 The trigger time for emergency maneuvers, and  The remaining distance before impact. Indeed, these two parameters encompass the most critical aspect of the performance and effects of pedestrian active safety systems because they affect the amount that the speed of the vehicle can be reduced prior to impact. In this study, the trigger times (relative to the time of impact) and the distance before impact at which a system is able to trigger are analyzed through the simulation of real world crash scenarios based on in-depth accident investigation.

Materials and Methods

In summary, the study uses a sample of crashes from in-depth accident databases. Each crash in the sample is reconstructed and modeled numerically. In parallel, a model representing the functioning of an AEB system has been established from a bibliographic research. Finally the system's model is coupled to the kinematic of the vehicle preceding the impact in order to evaluate the characteristics of the triggering of an AEB system in the sample of crashes.

Accident data

A sample of one hundred crashes was selected from two in-depth at-scene crash investigation databases: 60 cases were provided by CASR (Centre for Automotive Safety Research, University of Adelaide, Australia) and 40 cases were compiled from the IFSTTAR-LMA crash database (the laboratory of accident mechanism analysis of the French institute of science and technology for transport, development and networks, France). Both of these centres proceed in a similar way to perform in-depth investigations [START_REF] Girard | In-depth investigation of accidents: the experience of INRETS at Salon de Provence[END_REF][START_REF] Woolley | The Adelaide metropolitan indepth crash investigation study 2002-2005[END_REF]. In particular, these in-depth studies consist of investigations by a multidisciplinary team, composed of a psychologist and a technician, and are made on the scene of the accident, at the same time as the intervention of the rescue services. Those performing the survey were asked to collect a maximum amount of data. Generally, this data consisted of the final positions of the vehicles, the marks left on the ground (tyres, fluids, debris, etc.), the point of impact on the vehicle (bonnet, windscreen, etc.), the direction of the impact, the weather conditions, the obstructions and also to collect statements of drivers, witnesses and make a record of any injuries on the basis of the medical report, etc. All this data was then pooled and compared in order to make an initial reconstruction of the accident and to make hypotheses regarding the process involved: direction of travel of the pedestrian, speed of the vehicle, etc. When the reconstruction proposed is in agreement with all indications available, it is adopted as being the most probable scenario. To summarize the data of the reconstruction, a global synthesis of the accident is drawn up by the investigators relating in details the story of the accident.

IFSTTAR-LMA crash cases include investigated accidents occurred between 1995 and 2011 near the townships of Salon-de-Provence and Aix-en-Provence. CASR accidents were investigated in the Adelaide Metropolitan Area in the period April 2002 to October 2005.

Data collected during each investigation included:

-Photographs and videos of the crash scene and vehicles involved; -Statements of people involved in the crash, witnesses, and police; -Details of the road environment, involved vehicles and pedestrians; -Details of injuries from medical records; -A site diagram of the accident drawing to scale including the marks observed on the scene (skid, debris, blood, etc.), the final position of the vehicle and the pedestrian, the estimated impact location and the estimated trajectories of the different subjects involved in the crash. Some examples of these in-depth investigations have already been published [START_REF] Hamdane | Accident simulation and reconstruction for enhancing pedestrian safety: issues and challenges[END_REF].

The inclusion criteria of the cases used in this study were that the impact location on the roadway was known and the speed of the vehicle could be reliably assessed from standard crash reconstruction techniques.

The accidents were systematically clustered into four distinct groups according to the predominant features of the pre-crash sequence, as they might relate to the performance of the sensing system. These groups were: obstacles which hid the pedestrian, curved vehicle trajectory, pedestrian crossing from the sidewalk and pedestrian crossing from the far side of the roadway.

Accident modeling

The first step in each was to graphically reconstruct it. The approximate trajectories of the vehicle and the pedestrian were extracted from the scaled accident diagram provided from the in-depth investigation. Obstacles that mask pedestrians were located also using the diagram. Then, a temporal reconstruction was set up to emulate the kinematics of both the vehicle and the pedestrian from a pre-defined initial state until the impact. The initial conditions of the crash-involved pedestrian and vehicle were defined. As the pedestrian walking speed is often missing in the in-depth accident databases, it was estimated from the speed of the 50 th percentile based on the age of the pedestrians [START_REF] Huang | Evaluation of remote pedestrian sensor system based on the analysis of car-pedestrian accident scenarios[END_REF][START_REF] Montufar | Pedestrians' Normal Walking Speed and Speed When Crossing a Street[END_REF][START_REF] Zhang | Modeling pedestrian walking speed at signalized crosswalks considering crosswalk length and signal timing[END_REF] and considered as constant (Table 1). In particular, the speed change of the pedestrian when he/she starts walking has not been taken into account since this speed variation concerns less than the first half meter of his moving as shown by [START_REF] Zębala | Pedestrian acceleration and speeds[END_REF]. Similarly, the vehicle travel speed is also considered as constant over the chain of the precrash events in order to "rewind" and retrieve the position of the vehicle at the beginning of the accident scenario. This vehicle speed is given in the in-depth database and it is calculated or estimated using crash reconstruction methods [START_REF] Mclean | Vehicle travel speeds and the incidence of fatal pedestrian collisions[END_REF][START_REF] Serre | Detailed investigation and reconstructions of real accidents involving vulnerable road users[END_REF]. In cases where the driver did react by a hard braking, a constant braking deceleration was estimated for each accident case according to the road surface condition provided by the indepth database (for example, the deceleration is assumed to be -8m/s² for dry conditions and -6m/s² for bad conditions). Cases where the driver did brake gently were rejected from this study. Some other cases were added where the vehicle was stopped at an intersection before moving forward and colliding a pedestrian. In these cases, the acceleration of the vehicle was considered constant. When completed, the reconstruction includes the trajectories of both the vehicle and the pedestrian to locate in space and time the pedestrian relative to the vehicle for several seconds prior to the collision.

2.3. Active safety system modeling Most of the systems found on the market use both a camera vision mounted beside the rear view mirror and radar located in the front end of the vehicle behind the grille [START_REF] Coelingh | Collision Warning with Full Auto Brake and Pedestrian Detection -a practical example of Automatic Emergency Braking[END_REF][START_REF] Hayashi | Development of pre-crash safety system with pedestrian collision avoidance assist[END_REF][START_REF] Meinecke | SAVE-U: First Experiences with a Pre-Crash System for Enhancing Pedestrian Safety[END_REF][START_REF] Scheunert | Multi sensor based tracking of pedestrians: a survey of suitable movement models[END_REF]. Data flows from the on-board sensors are merged during the vision processing. Sensors described as "time-offlight sensors" are designed to determine the distance between the vehicle and any obstacle. Imaging sensors are mainly used for detecting and tracking pedestrians [START_REF] Gandhi | Pedestrian Protection Systems: Issues, Survey, and Challenges[END_REF]. Pedestrians are accurately detected by the cameras when they are located less than 40m from the vehicle [START_REF] Meinecke | Deliverable D6: Strategies in Terms of Vulnerable Road User Protection[END_REF]. Hence, in the detection algorithm, it is the camera sensor which is the most relevant because it validates the pedestrian presence and the possibility of a collision. Therefore, in this study, it was decided that the detection system could be modeled solely by the camera placed behind the rear view mirror of the vehicle and without taking into account the effect of others sensors such as radar or scanner. It has been also assumed in this study that the camera is not affected by bad weather conditions. This model is also acceptable for stereo cameras located near the rear-vision mirror of the vehicle since these cameras analyze the forward scene by overlapping the view of both lenses.

Camera sensors have been modeled in this study by their range and Field Of View (FOV). Different FOVs were considered from 20° to 45° in order to evaluate its influence on detection. The range was fixed to 40 m according to the research conducted in the European project SAVE-U [START_REF] Meinecke | Deliverable D6: Strategies in Terms of Vulnerable Road User Protection[END_REF]. The other characteristics of the camera like image processing or system lag are ignored since it is very difficult to get the required characteristics in the literature and consequently to model it numerically. Figure 1 illustrates a configuration of one camera as it has been considered in the accident modeling with a range and a FOV. Concerning the actuators for triggering an emergency maneuver, they may include autonomous emergency braking and/or steering capabilities. But because the purpose of this paper is focusing on the detection of the pedestrian and the evaluation of the remaining distance before impact, the actuator considered was a generic AEB system. The generic AEB system was modeled with an appropriate level of average deceleration of -8m/s² [START_REF] Brach | Vehicle accident analysis and reconstruction methods[END_REF][START_REF] Byatt | Manual of road accident investigation[END_REF][START_REF] Lechner | Analysis and reconstruction of accident sequences[END_REF]. This braking deceleration can be associated with "ideal" conditions: dry surface, good tires, efficient brake system and the like. The brake force was assumed to have a step response excluding the existing system lag and the transient state. Figure 2 displays this brake model and compares it to a conventional one.

Simulation process

Once the modeling had been set up, all the required components of an accident (crash environment, vehicle, pedestrian, and sensor) were simulated computationally. These components interacted in a customized virtual environment similar to the real world crash scenario (see Figure 1). The simulation emulated the trajectories of both the vehicle and the pedestrian starting from initial conditions until the crash using predefined time-steps. The program accounted for blind spots due to obstacles that mask the pedestrian and computed the instant when the pedestrian was detected by the sensor (that is to say located on the sensor's field of view but not hidden).

From the simulation, a set of data was extracted describing the characteristics of the accident at different pre-defined events. Data included the location of the pedestrian relative to the vehicle (the forward and lateral position), the speed of the vehicle (taking into account when the drivers accelerate or brake) and the Time-To-Collision (TTC). This last value is computed from the remaining travel distance before impact for the vehicle and its speed [START_REF] Horst | Time-to-collision and collision avoidance systems[END_REF].

In this paper, there is a focus on two specific events occurring in the crash sequence: the first instant when the pedestrian is visible by the sensor (considered as detection) and the last moment when the brakes need to be applied to avoid the crash. For the detection, it is a question of determining the moment when the pedestrian is entirely inside the camera's Fields Of View and was not masked by any obstacle. The ability of an AEB system to assist in complete collision avoidance is evaluated with reference to a time in the crash sequence that is the last time-to-brake (LTTB). It corresponds to the time when the vehicle is located at a distance d stop before the impact defined by the following equation:

𝑑 𝑠𝑡𝑜𝑝 = 𝑉 2 2 * |𝑎| + 𝑑 𝑜𝑓𝑓𝑠𝑒𝑡 (1)
where V is the vehicle travelling speed (m/s), a is the deceleration that fluctuates depending on the road conditions of the reconstructed accidents (m/s²), d offset is the vehicle longitudinal clearance from pedestrian set at 0.3 m. At this sequence of the crash, an estimation of the avoidance rate of collisions is established for each camera FOV. It is calculated by verifying if the pedestrian is visible at LTTB even if the pedestrian is not yet on the roadway; i.e. it is possible for a system to avoid the crash (depending on the duration that the pedestrian is visible preceding the LTTB which is not measured here). Otherwise, if the pedestrian is out of the FOV or inside but masked by an obstacle, the case is considered mitigated (according to the moment of visibility before the collision) or unavoidable.

Results

Accident reconstruction results

Considering the environment of the crashes, specific factors are challenging the detection of the pedestrians. These factors are clustered into 4 scenarios.

-Scenario 1: pedestrian obscured by an obstacle (traffic, parked vehicle, tree, bin…); -Scenario 2: pedestrian crossing through or after a turn; -Scenario 3: pedestrian starting crossing from the sidewalk or nearby (from the right of the vehicle in France, from the left in Australia); -Scenario 4: pedestrian crossing from the far side of the roadway (from the left of the vehicle in France, from the right in Australia). Figure 3 illustrates each of these scenarios and gives their percentage distribution in the accident database. In order to better understand the relative trajectories between the vehicle and the pedestrian all along the accident scenario, the kinematic parameters were located in space and time. Figure 4 gives for the 100 accident reconstructions, the pedestrian position relatively to the vehicle position at different times which are respectively 2.5s, 2s, 1.5s, 1s and 0.5s before the impact. This figure highlights in particular that until 1.5s, the positions of the pedestrians relatively to the vehicle are still scattered, and will largely not invoke any response from the AEB system.

Then between 1.5s and 0.5s, their positions are gathered mainly in front of the vehicle. This period between 1.5 s and 0.5 s TTC appears to be critical. Moreover, these results show that the lateral position of the pedestrian a short time before collision appears to be very important because at 1s before the impact, it can be observed that most of the pedestrians are located no more than 3m of the side of the vehicle. Concerning their longitudinal position relative to the vehicle, the pedestrians are mainly to be found less than 20 m far from the front of the vehicle within 1 s TTC. The travel speed of the vehicles ranged from 20 to 60 km/h with an average value of 40 km/h (S.D. 20.2 km/h). In 33% of cases, drivers perceive the hazard and then brake before the impact; the average impact speed was 32 km/h. More specifically, 23% of drivers applied the brakes at a TTC greater than 1 s after perceiving the pedestrian located at a lateral distance greater than 0.5 m from the side of the vehicle.

The average walking speed of an adult is about 1.4 m/s (SD = 0.39). Yet, there is a nonnegligible number of running pedestrians representing 24 cases; among that group, 62% are children which raise the running speed average to 3.5 m/s (SD = 0.74). The specific scenario of these crashes is running children masked by another vehicle. Half of these running children are coming from the sidewalk. Three cases concern children hit at the beginning of their travel.

Parameters relevant to pedestrian detection

The detection of the pedestrians in each case is characterized by the Time-To-Collision (TTC) and the pedestrian location relative to the vehicle (the longitudinal and lateral position) at the first instant of detection. These detection parameters were evaluated for different Fields-Of-View (FOV) of the sensor: 20°, 25°, 30°, 35°, 40° and 45°.

Figure 5 gives the complementary cumulative frequencies for each of the three kinematic parameters according to the different FOV. From a general point of view, all the curves never reach 100% since there are about 10% of pedestrians that remain undetected till the crash. These undetected pedestrians are mainly due to a pedestrian location outside the sight of view of the sensor. These cases are collisions on the side of the vehicle. Concerning the time remaining before the impact, about 60% of pedestrians are visible by sensors with a FOV over 35° at 2.5 s before impact (first time of the simulation). Because the number of visible pedestrian remains quite similar with upper FOV, it appears that the optimum FOV for the camera is above 35°. As important as it is, this optimized FOV is expected to detect beyond 80% of hazards 1s before impact.

About the remaining distance before collision, it can be highlight that for all FOV more than half of the pedestrians are visible 20m before the impact and so can be detected. The variation of the FOV has a main role only during the last 20m with a better visibility according to a wider FOV. An angle of 35° seems to be optimal again. It appears also that 10% of the pedestrians are never visible even with a FOV upper than 35°. 40m before the impact, only 10% of pedestrians are visible.

Considering the lateral distance, about 50% of pedestrians are visible with a FOV of 20° from a distance in lateral of 2m. With a FOV upper than 35°, this rate is reached when pedestrians are located above 3m. For a FOV below of 35°, approximately no pedestrians are visible from a lateral distance of 6m. If it is considered a pedestrian situated on the far side of the road (i.e. at about 4.5m in lateral), a FOV lower than 30° allows a visible rate of 10% while for upper FOV, this rate reaches 20%. This rate takes into account all visible pedestrians even if they are not yet on the roadway.

Parameters relevant to AEB performance

The kinematic parameters (TTC, longitudinal and lateral position of the pedestrian) are analyzed at the LTTB. This defines the requirements of an AEB system that can avoid the collision. The AEB performance refers to the crash avoidance rate obtained from the analysis of the 100 reconstructed crashes.

Figure 6 gives the cumulative frequency of the crash avoidance rate as a function of the kinematic parameters according to different sensor's FOV. It appears that 50% of accidents could be avoided if systems are able to be triggered 1s before the impact with a FOV upper than 35°. With a 20° FOV, the avoidance rate decreases to 40%. A threshold is reached at a LTTB equal to about 1.5 s; beyond this LTTB value, there is little improvement in the rate at which crashes can be avoided. Most challenges concerning crash avoidance occur between 0.5 s and 1.5 s (corresponding to the sharp slope of the curve). With a FOV of 35° or more, this threshold reaches an avoidance rate of 80%. Once more, it seems that this FOV is an optimal value as it includes pedestrians not yet on the roadway.

Regarding the longitudinal distance, if the system could be triggered more than 20 m before the impact, about 60% of accidents in this sample would have been avoided with a FOV of 20° (or 75% with a FOV of 35° or more). Beyond 20 m there seems to be no more gain in the rate of avoidance.

For the lateral distance, similar patterns are observed except that the threshold corresponds to a distance of approximately 3 m. The avoidance rate is affected by the FOV only at a distance beyond 2 m. At this lateral distance, about 50% of accidents can be avoided. In order to assess the available time for the system to react, the elapsed time between the instant when the pedestrian is visible (t visible ) and the LTTB was evaluated. Indeed this elapsed time corresponds to the duration available to detect the pedestrian and to trigger an AEB. This duration has been defined by Keller et al. (2011a) as the TTB (Time To Brake). Figure 7 shows the evolution of the complementary cumulative frequencies of avoided accidents according to this elapsed time and for the different FOV. These curves highlight that if the system is ideal, i.e. needs 0 second to react, the avoidance rate is comprised between 70% with a FOV of 20° and 83% with a FOV of 35°. With a FOV of 35° or more, if the system requires 0.5s for triggering (like the systems described by [START_REF] Edwards | Estimate of Potential Benefit for Europe of Fitting Autonomous Emergency Braking (AEB) Systems for Pedestrian Protection to Passenger Cars[END_REF], the avoidance rate decreases to 75%. If the reaction time is 1s (respectively 2s) the avoidance rate is 64% (respectively ~8%). It is then noted that effectiveness rapidly declines for reaction times above one second.

For fields-of-view above 35°, it is confirmed that the benefits of a system does not improve with increasing fields of view.

Discussion

This research is based on the understanding of the functioning of pedestrian active safety systems through real world crash reconstructions. A sample of 100 accidents involving pedestrians and passenger cars was gathered in order to illustrate the diversity of world crash configurations. Indeed, the sample is not only covering frontal impact configurations but also including side impacts along the fender of the vehicle.

Modelling the crash scenario depends of the availability of crash data required for the reconstruction. There is inevitably some fuzzy and even missing data to complete the reconstruction of a crash. This matter led to establish assumptions and define a procedure for modelling the different components of an accident in order to realise reliable computational simulations of the crash scenarios. To remove any doubt concerning the robustness of using data from in-depth investigation, a sensitivity study has been made on the assessments of two critical physical parameters: the speed of both pedestrian and vehicle. Results show that a variation of +/-10% of these speeds can be neglected on the accident reconstruction. For example, for an accident with a car speed of 11,12 m/s (40km/h) and a pedestrian walking speed of 1,4 m/s, a car speed variation of +/-5km/h and a pedestrian walking speed variation of +/-0,4 m/s induce a position of the pedestrian at LLTB of +/-0,4m on his trajectory. Similar analysis was conducted on the variance of the car deceleration and results show again a marginal effect.

The modeling of the detection function of an active safety system and its implementation in the computational simulation of crash scenarios allows the relevant parameters in detecting pedestrians to be highlighted. It has to be noted that the effect of lighting/weather conditions was not considered in the simulation of the crash scenarios. This constraint concerns more specifically the algorithms of detection and the limits of the technology. Additionally, no sensitive analysis was conducted in order to assess the false positive rate. This issue is also considered to be related to the algorithms of detection and probably differs from one active safety system to another. To include an assessment of the false positives to the simulation process, it suggests that the analysis should be conducted for a determined system. In the set of 100 real accidents selected for this research, it appears that not all the pedestrians are visible prior to the crash. The remaining set of not visible pedestrians (about 10%) corresponds mostly with scenarios where vehicles are turning (scenario S2) and obviously cases where pedestrians are masked by obstacles (scenario S1). These scenarios end up with side and front corner impact configurations.

Different FOV were evaluated for the sensor model. The evaluation highlighted that the rate of visible pedestrians is increasing with a wider FOV for the camera. This rate is about 79% for a 20° cone angle while it reaches 92% for a 45° cone angle. Additionally, from a 35° FOV, a threshold in the visibility rate is observed. For example, at 1s before impact about 80% of pedestrians could be detected with this FOV. Beyond it, the pattern of the visibility rate is similar for any kinematic parameter: the Time-To-Collision, the longitudinal and lateral pedestrian position relative to the vehicle. These results are complementary and in accordance with those of [START_REF] Rosén | Pedestrian injury mitigation by autonomous braking[END_REF] which show a slight reduction of the severely injured (as well as fatality) for camera sensors with a FOV from 40° to 180°. Thus, it can be considered that a FOV of 35° is relevant for pedestrian detection. An expansion of this work could be interesting to study the influence of a 35° FOV on the AEB for vehicle-to-vehicle collisions. Regarding the triggering issue of an active safety system, the braking maneuvers were considered in this research but not the steering maneuver because this last emergency procedure must be handled with caution. Indeed, it is restricted by a considerable number of factors such as the traffic situation and an available clearance enabling a safe evasive action. Furthermore, drivers' behavior and driving situation in a bend are additional constraints that can cancel the steering maneuver. Moreover, if a front-end collision can be avoided by steering, it can possibly lead to a side collision when the pedestrian keeps moving [START_REF] Hayashi | Autonomous collision avoidance system by combined control of steering and braking using geometrically optimised vehicular trajectory[END_REF].

The deployment of an AEB system (Autonomous Emergency Braking system) was assessed by computing the LTTB (Last Time To Brake). The system lag and the building rate of a full braking (i.e. the transient phase of an AEB system) were not considered in this research because not enough information was available. However, this delay can be taking into account in the processing time (Figure 7). For example, adding the braking system delay of 0.5s to a processing time of 0.5s leads to increase the elapsed time from the visibility of the pedestrian to the LTTB to 1s. So, according to Figure 7, for a FOV of 35°, it would decrease the avoidance rate from 75% to 63%. This rate could drop off considerably if the sum of the processing time and the delay is above 1s. Future systems will have to reduce these aforementioned parameters to improve their effectiveness.

Once the LTTB is calculated, it is possible to determine the effect of the FOV of a camerabased system in terms of avoidance or mitigated cases. As for detection, the avoidance rate increase with the FOV until reaching a threshold at 35°. This FOV can avoid approximately 50% of crashes if the system trigger at 1s before impact. This rate is a little overestimated compared to the literature. [START_REF] Lindman | Benefit estimation model for pedestrian auto brake functionality[END_REF] presented in a Case Study based on accident data from GIDAS the potential effectiveness of an active safety system developed by Volvo Cars (the CWAB-PD). It was estimated that CWAB-PD autonomous braking could avoid about 30% of all pedestrian accidents.

Concerning the reaction time of an active safety system (from detection of the hazard to the decision making and deployment of the emergency maneuver), the elapsed time between the instant when the pedestrian is visible and the LTTB was studied. During this elapsed time, it could be expected to alert the driver by a warning system (visual and/or audible alarm). Since a driver needs in average 1s to react [START_REF] Lee | Collision warning timing, driver distraction, and driver response to imminent rear-end collisions in a high-fidelity driving simulator[END_REF], only cases with an available time more than 1s can be considered. Such cases represent 50% of our database. So, among this rate, it is possible to avoid some accidents by prompting a response from the driver. Finally, the results of this research can be used to establish specifications for an active safety system. Indeed, according to the Figure 7, it can be observed that the objective of 75% of avoidance rate require 0,5s of reaction time for a 35° FOV.

Conclusion

The purpose of this paper was to highlight issues and challenges in pedestrian active safety through reconstruction of 100 real accident cases involving a vehicle and a pedestrian. The reconstructed crashes provided information concerning the crash configurations, the travel and impact speeds of vehicles, the velocity of pedestrians and their relative position in time and space, etc. In particular, the time between 0.5 and 1s appeared as critical for the pedestrian safety.

The functionality of AEB systems were analyzed according to the reconstructed crash cases. The analysis considered different attribute of a system: the detection, reaction and triggering of the brakes. According to these attribute, a system based on a camera with a FOV of 35° appeared relevant in terms of efficiency: high rates in detection (~90%) and avoidance (~75%). Depending on the direction taken in using the graphs established in this research, it is possible to evaluate the efficiency of active safety systems as well as to define their specifications.

Next work will concern the evaluation of active safety systems in terms of crash mitigation (i.e. their potential effects on casualty reduction). The decrease of the impact speed induced by these systems will be assessed in order to identify the potential benefit of these systems on casualty reduction. This analysis will be carried using risk curves derived by logistic regression of the accident data. A coupled approach joining active and passive safety systems can be added in the future to establish assessed risk curves.

Figure 1 :

 1 Figure 1: Scheme illustrating a crash representation including the active system

Figure 2 :

 2 Figure 2: Comparison between the current and conventional brake model

  Figure 3: Description of the 4 crash scenarios

Figure 4 :

 4 Figure 4: Pedestrian positions relatively to the vehicle at 2.5s, 2s, 1.5s, 1s and 0.5s before impact

Figure 5 :

 5 Figure 5: Rate of visible pedestrians for each kinematic parameter according to different FOVs

Figure 6 :

 6 Figure 6: Pedestrian avoidance rate in function of kinematic parameters according to the camera FOV and the vehicle braking

Figure 7 :

 7 Figure 7: Complementary cumulative frequencies of the avoided accidents in function of the elapsed time from the visibility of the pedestrian to the LTTB for different FOV

Table 1 Pedestrian speed according to the age (from [Huang et al, 2008])

 1

	Age	50% speed (m/s)
		Walking	Running
	5-9	1.83	3.94
	10-14	1.68	4.20
	15-19	1.65	4.20
	20-29	1.62	3.54
	30-39	1.62	3.35
	40-44	1.62	2.90
	45-49	1.52	2.90
	50-54	1.52	2.83
	55-59	1.46	2.83
	60-64	1.46	2.47
	65+	1.28	2.47

Acknowledgements

The authors would like to thank the both in-depth investigation teams of CASR (Centre for Automotive Safety Research, the University of Adelaide, Australia) and IFSTTAR-LMA (the Laboratory of accident mechanism analysis of the French institute of science and technology for transport, development and networks, Salon-de-Provence, France) for their help and cooperation in the reconstruction of the accidents.


P O Goffard

S Loisel

D Pommeret

Polynomial approximations for bivariate aggregate claims amount probability distributions

Keywords: Bivariate aggregate claims model, bivariate distribution, bivariate Laplace Transform, numerical inversion of Laplace transform, natural exponential families with quadratic variance functions, orthogonal polynomials

A numerical method to compute bivariate probability distributions from their Laplace transforms is presented. The method consists in an orthogonal projection of the probability density function with respect to a probability measure that belongs to a Natural Exponential Family with Quadratic Variance Function (NEF-QVF). A particular link to Lancaster probabilities is highlighted. The procedure allows a quick and accurate calculation of probabilities of interest and does not require strong coding skills. Numerical illustrations and comparisons with other methods are provided. This work is motivated by actuarial applications. We aim at recovering the joint distribution of two aggregate claims amounts associated with two insurance policy portfolios that are closely related, and at computing survival functions for reinsurance losses in presence of two non-proportional reinsurance treaties.

Introduction

In insurance and in reinsurance, some common events may cause simultaneous, correlated claims in two lines of business. For example, in third-party liability motor insurance, an accident may cause corporal damage and material damage losses. These parts of the same claim are then handled by different claim managers and reserving is most often also done separately, which makes it necessary to study the joint distribution and not only the marginals or the distribution of the sum (except when one loss dominates the other one, which makes correlation of secondary importance). The study of the sum may be reduced to a one-dimensional problem, see [START_REF] Ambagaspitiya | On the distribution of a sum of correlated aggregate claims[END_REF]. Similarly, a reinsurer who accepts two stop-loss treaties from two customers operating on the same market is exposed to the sum of two excesses of aggregate losses, which contain in general some independent part and some common shock part, arising from events that generate claims for the two insurance companies at the same time. To compute the Value-at-Risk of the reinsurer's exposure, or to compute the probability that the reinsurer looses more than certain risk limits, we need the bivariate distribution of the aggregate claims amount of the two insurers. In addition to classical couples of independent aggregated losses drawn from compound distributions, one needs to study the common shock part of the bivariate loss vector. We define a bivariate collective model with common shock by

S 1 S 2 = N j=1 U 1j U 2j + N 1 i=1 V i N 2 i=1 W i , (1.1)
= X 1 X 2 + Y 1 Y 2 , (1.2)
where S = (S 1 , S 2 ) denotes the aggregate claims amount of two non-life insurance portfolios over a given time period. An event may be associated to a unique claim incurred simultaneously in the portfolios. The number of such events is modeled through a counting random variable N . These events cause losses that are positively correlated and represented by a sequence of independent and identically distributed (i.i.d.) random vectors {(U 1j , U 2j )} j∈N * . We are interested in bivariate distributions that allow dependence between U 1j and U 2j .

An event may also cause several losses in the two portfolios. The number of claims are then correlated and modeled by two counting random variables N 1 and N 2 . The vector N = (N 1 , N 2 ) admits a bivariate distribution allowing dependency between N 1 and N 2 .

The claim sizes are sequences of i.i.d. random variables {V i } i∈N * and {W i } i∈N * . The model (1.1), proposed in the work of Jin et al. [START_REF] Jin | Recursions and fast fourier transforms for a new bivariate aggregate claims model[END_REF], enables one to include dependency among both the claim sizes and the claim frequencies. Aggregate claims are governed by compound distributions. These distributions, in the univariate case, are already difficult to study because closed formulae for the probability density function are available in only a few cases. The cumbersome computation time induced by Monte Carlo simulations to obtain accurate values of probabilities motivates the use of numerical techniques.

Recursion based methods have been designed to cope with bivariate aggregate claims distributions, they can be viewed as extension of the so-called Panjer's algorithm. The distribution of the vector X = (X 1 , X 2 ) is recovered through a recursion procedure in [START_REF] Ambagaspitiya | On the distributions of two classes of correlated aggregate claims[END_REF][START_REF] Sundt | On multivariate Panjer recursions[END_REF], whereas the distribution of the vector Y = (Y 1 , Y 2 ) is studied in [START_REF] Hesselager | Recursions for certain bivariate counting distributions and their compound distributions[END_REF][START_REF] Sundt | On multivariate Panjer recursions[END_REF][START_REF] Vernic | Recursive evaluation of some bivariate compound distibutions[END_REF]. Recursion based techniques rely on the existence of recurrence relationships between the probabilities of the claim frequencies which limits their use to certain kind of distributions. The Fast Fourier Transform is a procedure that permits one to inverse the discrete Fourier transform of the aggregate claim distribution. A multivariate extension is derived and applied to recover the distribution of Y and S in [START_REF] Jin | Recursions and fast fourier transforms for certain bivariate compound distributions[END_REF][START_REF] Jin | Recursions and fast fourier transforms for a new bivariate aggregate claims model[END_REF]. The authors provide a numerical comparison of the two procedure. The required arithmetization of the claim amounts distribution remains the main drawback of these methods. One may be interested in methods that do not start with a discretization of the claim size distribution.

In this paper, we develop an effective method to compute bivariate Probability Density Functions (BPDFs) from the knowledge of their bivariate Laplace transforms or equivalently from the moments of the distribution. First, we choose a reference distribution that belongs to a NEF-QVF. Then, we express the BPDF with respect to this reference measure as an expansion in terms of orthonormal polynomials. The polynomials are orthogonal with respect to the reference probability measure. When we aim at recovering a BPDF on the positive quadrant of the plane, as it is the case within the frame of aggregate claims amount distribution, our method looks like a well established method called the Laguerre method. The similarity is due to the fact that we choose a Gamma distribution as reference distribution associated with Laguerre polynomials to ensure the validity of the expansion. We refer to [START_REF] Abate | On the Laguerre method for numerically inverting Laplace transforms[END_REF] where an effective variant of the Laguerre method is presented. The same authors also work out a multivariate extension, see [START_REF] Abate | Numerical inversion of multidimensional Laplace transform by the Laguerre method[END_REF]. The main contribution of our work lies in the possibility of choosing the parameters of the reference distribution to enhance the performance of the method. A great body of the applied probability literature is dedicated to the numerical inversion of the Laplace transform of univariate PDF but to the best of our knowledge only a few attempts have been made in a multivariate context. A Fourier series based method is presented in [START_REF] Choudhury | Multidimensional transform inversion application to the transient M / G / 1 queue[END_REF] motivated by applications in queueing theory. The desired function is recovered from the Fourier series coefficients of a periodic version of the function constructed by aliasing. The inversion formula involves integral computations that are completed using a trapezoidal rule justified by Poisson summation formula. It also implies truncation but not a simple one: Euler summation formula for alternating series is employed consequently to the neat replacement of the original infinite serie by an equivalent alternating one. Interesting applications are allowed because of the possibility of recovering BPDFs of couples of random variables formed by a discrete and a continuous one. Recently, a moment-recovered approximation has been proposed in [START_REF] Mnatsakanov | Moment-recovered approximations of multivariate distributions: The Laplace transform inversion[END_REF]. The stable approximants for the BPDF and the joint survival function are derived. A very interesting feature of this method is that the approximation can be turned into a statistical estimation as there exist empirical counterparts of the aforementioned approximants. The expansion we propose is also well adapted to empirical estimation. The empirical estimation of BPDF when data are available will be at the center of a forthcoming paper. Finally, one key aspect of the method is that once the parameters have been chosen and the coefficients computed, approximations for the joint probability density and survival functions are available in an explicit form.

In Section 2, we introduce a BPDF expansion formula based on orthogonal projection within the frame of NEF-QVF. In Section 3, the expansion for bivariate aggregate claims amount distributions is derived along with conditions to ensure the validity of the polynomial expansion. Section 4 presents the link established between our expansion and Gamma-Gamma Lancaster probabilities. Section 5 is devoted to numerical illustrations and examples of application. A comparison to the moment-recovered method is also proposed.

Expression of the joint density

Orthogonal polynomials associated with NEF-QVF

Let F = {F θ , θ ∈ Θ} with Θ ⊂ R be a Natural Exponential Family (NEF), see [START_REF] Barndorff-Nielsen | Information and Exponential Families in Statistical Theory[END_REF], generated by a probability measure ν on R such that

F θ (X ∈ A) = A exp [xθ -κ(θ)] dν(x) = A f (x, θ)dν(x),
where A ⊂ R, κ(θ) = log R e θx dν(x) is the Cumulant Generating Function (CGF), and f (x, θ) is the density of P θ with respect to ν. Let X be a random variable

F θ distributed. The mean of X is µ = E θ (X) = xdF θ (x) = κ (θ),
and its variance is

V(µ) = V ar θ (X) = (x -µ) 2 dF θ (x) = κ (θ).
The application θ → κ (θ) is one to one. Its inverse function µ → h(µ) is defined on M = κ (Θ). With a slight change of notation, we can rewrite F = {F µ , µ ∈ M}, where F µ has mean µ and density g(x, µ) = exp {h(µ)x -κ [h(µ)]} with respect to ν. A NEF has a Quadratic Variance Function (QVF) if there exist reals v 0 , v 1 , v 2 such that

V(µ) = v 0 + v 1 µ + v 2 µ 2 . (2.1)
The Natural Exponential Families with Quadratic Variance Function (NEF-QVF) include the normal, Gamma, hyperbolic, Poisson, binomial and negative binomial distributions. Define

P n (x, µ) = V n (µ) ∂ n ∂µ n g(x, µ) /g(x, µ), (2.2)
for n ∈ N. Each P n (x, µ) is a polynomial of degree n in both µ and x. Moreover, if F is NEF-QVF, {P n } n∈N is a family of orthogonal polynomials with respect to P µ in the sense that

< P n , P m >= P n (x, µ)P m (x, µ)dF µ (x) = δ nm ||P n || 2 , m, n ∈ N,
where δ mn is the Kronecker symbol equal to 1 if n = m and 0 otherwise. For the sake of simplicity, we choose ν = F µ 0 . Then, f (x, µ 0 ) = 1 and we write

P n (x) = P n (x, µ 0 ) = V n (µ 0 ) ∂ n ∂µ n g(x, µ) µ=µ 0 . (2.3)
We also consider in the rest of the paper a normalized version of the polynomials defined in (2.3) with Q n (x) = Pn(x) ||Pn|| . For an exhaustive review regarding NEF-QVF and their properties, we refer to [START_REF] Morris | Natural Exponential Families with Quadratic Variance Functions[END_REF].

Bivariate probability measures and their Laplace transform

Let (X, Y ) be a couple of random variables, governed by a bivariate probability measure F X,Y . We also assume that (X, Y ) admits a joint probability density function with respect to a bivariate measure denoted by λ, such that dF X,Y (x, y) = f X,Y (x, y)dλ(x, y).

(2.4)

Its bivariate Laplace transform is defined as f X,Y (s, t) = e sx+ty f X,Y (x, y)dλ(x, y).

(2.5)

We consider the bivariate NEF-QVF generated by

dν(x, y) = dν 1 (x) × dν 2 (y), (2.6)
where ν 1 and ν 2 are two probability measures that belong to NEF-QVF absolutely continuous with respect to λ. With the notation of Section 2.1, the probability density function associated with ν is

f ν (x, y) = f ν 1 (x)f ν 2 (x), (2.7)
where

f ν i denotes the PDF of ν i for i = 1, 2. Let Q ν k,l (x, y) = Q ν 1 k (x)Q ν 2 l (y)
, k, l ∈ N be a bivariate orthonormal polynomial system in the sense that

< Q ν k,l , Q ν i,j > = Q ν k,l (x, y)Q ν i,j (x, y)dν(x, y) < Q ν k,l , Q ν i,j > = δ ki δ lj .
We denote by L 2 (ν) the set of square integrable functions with respect to ν.

Proposition 1. Assume that dF X,Y dν ∈ L 2 (ν), we have the following expansion dF X,Y dν (x, y) = +∞ k=0 +∞ l=0 dF X,Y dν , Q ν k,l Q ν k,l (x, y) (2.8) = +∞ k=0 +∞ l=0 E Q ν k,l (X, Y ) Q ν k,l (x, y).
(2.9)

Applying Proposition 1 yields an expansion of the probability density function

f X,Y (x, y) = +∞ k=0 +∞ l=0 f X,Y (x, y)Q ν k,l (x, y)dxdy × Q ν k,l (x, y)f ν (x, y).
(2.10)

We have

Q ν k,l (x, y) = k i=0 l j=0 q i,j (k, l)x i y j , (2.11)
where q i,j (k, l) is expressed in terms of the coefficents of

Q ν 1 k and Q ν 2 l . Reinjecting the expression of Q ν k,
X,Y f X,Y (x, y) = +∞ k=0 +∞ l=0 k i=0 l j=0 q i,j (k, l) x i y j f X,Y (x, y)dxdy × Q ν k,l (x, y)f ν (x, y) = +∞ k=0 +∞ l=0 k i=0 l j=0 q i,j (k, l) ∂ i+j ∂s i ∂t j f X,Y (s, t) s=0,t=0 Q ν k,l (x, y)f ν (x, y) . = +∞ k=0 +∞ l=0 a k,l Q ν k,l (x, y)f ν (x, y).
Remark 1. The probability density function takes the form of an expansion. The coefficients are expressed in terms of the moments and equivalently in terms of the derivative of the bivariate Laplace transform. The expression of the probability density function is therefore a Laplace transform inversion formula.

Approximations of the bivariate probability density function are obtained through truncations of the infinite series

f K,L X,Y (x, y) = K k=0 L l=0 a k,l Q ν k,l (x, y)f ν (x, y),
where K and L denote the order of truncation of the approximation. In view of future applications, we are interested in getting approximations of cumulative distribution functions, which are obtained through integration of the approximated density

F K,L X,Y (x, y) = x -∞ y -∞ f K,L X,Y (u, v)dλ(u, v) (2.12) = K k=0 L l=0 a k,l x -∞ y -∞ Q ν k,l (u, v)f (u, v)dλ(u, v). (2.13)
The goodness of the approximation lies in the decay of the sequence defined by the coefficients of the expansion. Note that from

dF X,Y dν ∈ L 2 (ν), a Parseval type relation holds f 2 X,Y (x, y)dxdy = +∞ k=0 +∞ l=0 a 2 k,l < +∞. (2.14)
Relation (2.14) implies that the sequence {a k,l } k,l∈N decreases towards 0 with respect to k and l. The choice of the NEF-QVF and its parameters depends on the problem, and are discussed in the next section.

3 Application to a bivariate aggregate claims amount distribution

General formula

Aggregate claims amount are modeled by non-negative random variables. The only continuous NEF-QVF having support on R + is generated by the Gamma distribution, according to the classification given in [START_REF] Morris | Natural Exponential Families with Quadratic Variance Functions[END_REF]. The two other continuous NEF-QVF distributions are the normal and the hyperbolic ones which have support on R. The Gamma distribution admits a probability density function with respect to the Lebesgue measure defined as

f (x) = x r-1 e -x/m Γ(r)m r . (3.1)
The orthogonal polynomials associated with the Gamma distribution are the generalized Laguerre polynomials

Q k (x) = (-1) k L r-1 k (x/m) k+r-1 k , (3.2)
where

L r-1 n (x) = n i=0 n+r-1 n-i (-x) i
i! , see [START_REF] Szegö | Orthogonal Polynomials[END_REF]. The random vector (S 1 , S 2 ), defined in (1.1), is divided into two parts that are considered separately. The first part,

X 1 X 2 = N j=1 U 1j U 2j , (3.3)
is governed by a bivariate compound distribution. The bivariate probability measure associated with (X 1 , X 2 ) is divided into a discrete part with an atom at (0, 0) with a probability mass equal to p 0 = P (N = 0) and a continuous part that is absolutely continuous with respect to the bivariate Lebesgue measure. Thus we have

dF X 1 ,X 2 (x 1 , x 2 ) = p 0 δ {x 1 =0,x 2 =0} (x 1 , x 1 ) + dG X 1 ,X 2 (x 1 , x 2 ), (3.4)
where

δ {x 1 =0,x 2 =0} (x 1 , x 2 )
is the Dirac measure equal to 1 at (0, 0) and 0 otherwise. We aim at expanding the continuous part using Proposition 1. The defective density probability function of dG X 1 ,X 2 is denoted by g X 1 ,X 2 . We define a bivariate NEF-QVF probability measure as in Section 2 by

ν(x 1 , x 2 ) = ν 1 (x 1 ) × ν 2 (x 2 ), (3.5)
where ν 1 and ν 2 are Gamma distributions with parameters (m 1 , r 1 ) and (m 2 , r 2 ) respectively. Orthogonal polynomials with respect to ν are defined as

Q ν l,k (x 1 , x 2 ) = Q ν 1 k (x 1 ) × Q ν 2 l (x 2 ), (3.6)
where {Q ν 1 k } k∈N and {Q ν 2 l } l∈N are orthogonal polynomial systems with respect to ν 1 and ν 2 respectively, and defined as in (3.2). If

dG X 1 ,X 2 dν ∈ L 2 (ν), applying Proposition 1 yields g X 1 ,X 2 (x 1 , x 2 ) = +∞ k=0 +∞ l=0 a k,l Q ν k,l (x 1 , x 2 )f ν (x 1 , x 2 ), (3.7)
where f is the BPDF associated to the probability measure ν.

Remark 2. The Laguerre series method in [START_REF] Abate | Numerical inversion of multidimensional Laplace transform by the Laguerre method[END_REF] is related to our approach with a particular parametrization that is

m 1 = m 2 = 2 and r 1 = r 2 = 1.
We know from the Parseval type relation (2.14) that the sequence of coefficients of the expansion is decreasing towards 0. But how fast does it actually decrease?

We address this problem using generating function theory just like what is done in [START_REF] Abate | Numerical inversion of multidimensional Laplace transform by the Laguerre method[END_REF]. The interesting fact is that the generating function of the coefficients of the expansion is linked to the bivariate Laplace transform of the expanded probability density function. We start from the expression given in (3.7). By taking the bivariate Laplace transform, we get

g X 1 ,X 2 (s 1 , s 2 ) = +∞ 0 +∞ 0 e s 1 x 1 +s 2 x g X 1 ,X 2 (x 1 , x)dx 1 dx 2 = +∞ k=0 +∞ l=0 a k,l k + r 1 -1 k l + r 2 -1 l × 1 1 -s 1 m 1 r 1 1 1 -s 2 m 2 r 2 s 1 m 1 1 -s 1 m 1 k s 2 m 2 1 -s 2 m 2 l = 1 1 -s 1 m 1 r 1 1 1 -s 2 m 2 r 2 B s 1 m 1 1 -s 1 m 1 , s 2 m 2 1 -s 2 m 2 , (3.8)
where

B(z 1 , z 2 ) = +∞ k=0 +∞ l=0 b k,l z k 1 z l 2 (3.9)
is the generating function of the two-dimensional sequence {b k,l } k,l∈N . The coefficients of the expansion follow from

a k,l = b k,l k+r 1 -1 k l+r 2 -1 l . (3.10)
The generating function B is derived from the bivariate Laplace transform after a simple change of variable in Equation (3.8):

B(z 1 , z 2 ) = (1 + z 1 ) -r 1 (1 + z 2 ) -r 2 g Y 1 ,Y 2 z 1 m 1 (1 + z 1 ) , z 2 m 2 (1 + z 2 ) . (3.11)
In [START_REF] Abate | Numerical inversion of multidimensional Laplace transform by the Laguerre method[END_REF], the authors compute the coefficients of the expansion through the derivative of their generating function. They use Cauchy contour integral to express the derivatives as integrals and approximate them via a trapezoidal rule. In some cases, when the Laplace transform expression is simple, the derivatives can be calculated directly with a computational software. The attractive feature of our method lies in the possibility of tuning the parameters of the Gamma distributions in order to modify the generating function. We illustrate this fact through the following example:

Example 1. Let (X, Y ) be a couple of independent random variables exponentially distributed with parameters δ 1 and δ 2 . The bivariate Laplace transform associated with this couple of random variables is

f X,Y (s 1 , s 2 ) = δ 1 δ 1 -s 1 × δ 2 δ 2 -s 2 . (3.12)
We inject the bivariate Laplace transform expression (3.12) in the generating function defined in (3.11) to get

B(z 1 , z 2 ) = (1 + z 1 )(1 + z 2 ) (1 + z 1 ) r 1 (1 + z 2 ) r 2 × δ 1 δ 2 [δ 1 + z 1 (δ 1 -1/m 1 )] [δ 2 + z 2 (δ 2 -1/m 2 )] . (3.13)
In this case, the parametrization is straightforward. We set

r 1 = r 2 = 1, m 1 = 1/δ 1 and m 2 = 1/δ 2 , so B(z 1 , z 2 ) = 1 which implies that a 0,0 = 1 and a k,l = 0 for all k, l ≥ 1.
The Laplace transforms of the distributions considered in the rest of the paper are more complicated and it is less straightforward to choose the parameters.

We now turn to the study of the second part of the random vector (S 1 , S 2 ). Recall that

Y 1 Y 2 = N 1 i=1 V i N 2 i=1 W i . (3.14)
The counting random variables N 1 and N 2 are assumed to be independent. Thus the components of the random vector (Y 1 , Y 2 ) are also independent which simplifies the model stated in the introduction. The probability measure of (Y 1 , Y 2 ) is given by

dF Y 1 ,Y 2 (y 1 , y 2 ) = P (N 1 = 0)P (N 2 = 0)δ {y 1 =0,y 2 =0} (y 1 , y 2 ) + dG Y 1 (y 1 ) × dG Y 2 (y 2 ) + P (N 1 = 0)dG Y 2 (y 2 ) + P (N 2 = 0)dG Y 1 (y 1 ).
The defective probability measures dG Y 1 and dG Y 2 are associated to univariate compound distribution. The expression of the probability density function are not always available in a closed form. In that case, we can use a univariate polynomial approximation described in [START_REF] Goffard | A polynomial expansion to approximate the ultimate ruin probability in the compound Poisson ruin model[END_REF][START_REF] Goffard | Approximations polynomiales de densités de probabilité et applications en assurance[END_REF] to recover the probability function.

Finally, the polynomial approximation of the probability measure of (S 1 , S 2 ) follows from the combination of the polynomial approximations of the BPDF of (X 1 , X 2 ) and (Y 1 , Y 2 ) through a two-dimensional convolution procedure. In the next Section, we show how to choose the parameters of the expansion in order to ensure the integrability condition mentionned in Proposition 1.

Integrability condition

The polynomial expansion is valid under the following condition:

+∞ 0 +∞ 0 dG X 1 ,X 2 dν (x 1 , x 2 ) 2 dν(x 1 , x 2 ) < +∞, (3.15)
which is equivalent to

+∞ 0 +∞ 0 g X 1 ,X 2 (x 1 , x 2 ) 2 x 1-r 1 1 x 1-r 2 2 e x/m 1 +y/m 2 dx 1 dx 2 < +∞. (3.16)
The parameters m 1 , m 2 , r 1 and r 2 of the reference distribution play a key role regarding the validity of the polynomial expansion.

Let (X, Y ) be a couple of non-negative random variables with a joint PDF f X,Y (x, y). We define the set

E = {(s, t) ∈ [0, +∞) × [0, +∞); f X,Y (s, t) < +∞}, (3.17)
and denote by E c its complementary. The following result gives us insight on how to choose the parameters of the reference distribution.

Theorem 1. Let (X, Y ) be a couple of non-negative random variables with a joint PDF f X,Y (x, y). We assume that E = {(0, 0)}. We assume additionally that there exist two real numbers a, b ≥ 0 such that for all (x, y) ∈ [a, +∞) × [b, +∞), the applications x → f X,Y (x, y) and y → f X,Y (x, y) are strictly decreasing. Then we have

f X,Y (x, y) ≤ A(s, t)e -sx-ty , ∀(x, y) ∈ [a, +∞) × [b, +∞), (3.18)
where A(s, t) is some real number independent of x and y.

Proof. Let (x, y) ∈ [a, +∞) × [b, +∞), we have f X,Y (s, t) ≥ x a y b e su+tv f X,Y (u, v)dudv = x a e su 1 t e tv f X,Y (u, y) -e tb f X,Y (u, b) dudv - x a e su y b e tv t ∂f X,Y (u, v) ∂v dudv ≥ x a e su 1 t e tv f X,Y (u, y) -e tb f X,Y (u, b) dudv ≥ e sx+ty st f X,Y (x, y) + e sa+tb st f X,Y (a, b) - 1 st e sa+ty f X,Y (a, y) + e sx+tb f X,Y (x, b) . (3.19)
In order to deal with the part (3.19), we need the following lemma.

Lemma 1. Let f be a continuously differentiable function on R + . Assume that there exist some a > 0 such that x → f (x) is strictly decreasing for x > a. Assume additionaly that there exist s > 0 such that f (s) < +∞. Then we have

f (x) ≤ A(s)e -sx ∀(x, y) ∈ [a, +∞), (3.20)
where A(s) is some real number independent of x.

Proof. For x ≥ a, we have e ty f X,Y (x, y) < +∞ for all x ∈ [0, +∞), and for x = a in particular. We can therefore apply Lemma 1 to the application y → f X,Y (a, y), which yields f X,Y (a, y) ≤ A 1 (t)e -ty .

f (s) ≥
(3.24)

Using the same arguments, we get that

f X,Y (x, b) ≤ A 2 (t)e -sx . (3.25)
Reinjecting inequalities (3.24) and (3.25)

in (3.19) yields f X,Y (x, y) ≤ st f X,Y (-s, -t) -e sa+tb f X,Y (a, b) + A 2 (t)e tb + A 1 (s)e sa e -sx-ty . (3.26)
Theorem 1 allows us to bound continuous joint PDF when the bivariate Laplace transform is well defined for some positive arguments. We aim at recovering the distribution of (X 1 , X 2 ) that admits a singular and a continuous part. We put aside the singular part and expand the continuous part represented through its defective bivariate defective probability density function g X 1 ,X 2 . Theorem 1 holds for g X 1 ,X 2 as for any common bivariate probability density function that satisfies the hypothesis of Theorem 1. We define E * = E ∩ E c , where E is the closure of E. One can note that for any (s, t) ∈ [0, +∞) × [0, +∞) if there exist (s * , t * ) such that s < s * and t < t * then (s, t) ∈ E. The next result follows from that last remark. Corollary 1. We consider a couple of random variables (X 1 , X 2 ) defined as in (3.3). Assume that g X 1 ,X 2 satisfies the hypothesis of Theorem 1. For all

(s * 1 , s * 2 ) ∈ E * , if we take {1/m 1 < 2s * 1 , 1/m 2 < 2s * 2 , r 1 = 1, r 2 = 1}, then the integrability condition (3.16) is satisfied.
Corollary 1 will help us greatly in Section 5. Theorem 1 and Corollary 1 show that the polynomial expansion is valid under the condition that the bivariate Laplace transform exists. If it is not the case, it is difficult to ensure the good behavior of the coefficients of the polynomial expansion. It is the case for instance when the claim sizes admit an heavy tailed distribution. In those cases, the method proposed in this paper might not applicable. Within the frame of insurance applications, a solution can be to consider truncated distributions for the claim sizes justified by the existence of a non proportional reinsurance treaty with a given retention level. The next section links our approach with Lancaster probabilities.

Downton Bivariate Exponential distribution and Lancaster probabilities

A classical problem, which has become known as the Lancaster problem, see [START_REF] Lancaster | The structure of bivariate distributions[END_REF], is to characterize bivariate distributions with given margins and orthogonal eigenfunctions. Let ν 1 and ν 2 be two probability measures with associated orthonormal polynomials sequences {Q 1 n } n∈N and {Q 2 n } n∈N respectively. We denote by f ν 1 and f ν 2 the PDF of ν 1 and ν 2 . Exchangeable Lancaster bivariate densities have the form

f X,Y (x, y) = +∞ k=0 a k Q ν 1 k (x)Q ν 2 k (y) f ν 1 (x)f ν 2 (y), (4.1)
where {a k } k∈N is a sequence of real numbers with a 0 = 1. The Lancaster problem can be reduced to finding a proper sequence {a k } k∈N such that f is nonnegative and therefore a BPDF. The characterization of the Lancaster probabilities constructed via NEF-QVF has been widely studied by Koudou, see [START_REF] Koudou | Problèmes de marges et familles exponentielles naturelles[END_REF][START_REF] Koudou | Probabilités de Lancaster[END_REF][START_REF] Koudou | Lancaster bivariate probability distributions with Poisson, negative binomial and gamma margins[END_REF]. He gave existence conditions and explicit forms of Lancaster sequences in various cases. In a recent publication, Diaconis et al. [START_REF] Diaconis | Exchangeables pairs of Bernouilli random variables, Krawtchouck polynomials, and Ehrenfest urns[END_REF] studied the binomial-binomial case which means that the marginal distributions are both binomial. Characterizations with a more "probabilistic flavor" than the classical one are given. Recall that the classical way to construct a Lancaster probability consists in considering the distribution of (X, Y ) = (U +W, V +W ), where U, V, W , are independent random variables governed by the same distribution. This "random element in common" way allows to construct Lancaster probabilities based on any NEF-QVF distribution. To connect Lancaster probabilities with our approach, we start by considering our representation of the BPDF for a couple of nonnegative random variables (X, Y )

f X,Y (x, y) = +∞ k=0 +∞ l=0 E Q ν 1 k (X)Q ν 2 l (Y ) Q ν 1 k (x)Q ν 2 l (y) f ν 1 (x)f ν 2 (y), (4.2)
where f ν i is a Gamma PDF and {Q i k } k∈N its orthonormal polynomials sequence for i = 1, 2. In order to get a Lancaster probability from the BPDF defined in (4.2), we need a biorthogonality condition, namely

E Q ν 1 k (X)Q ν 2 l (Y ) = δ kl a k . (4.3)
Let (X, Y ) be a couple of random variables governed by a Downton Bivariate Exponential distribution DBV E(µ 1 , µ 2 , ρ) introduced in [START_REF] Dowton | Bivariate exponential distributions in reliabilty theory[END_REF]. We apply our method to recover the joint PDF of this distribution. The bivariate Laplace transform associated with the DBVE distribution is

f X,Y (s 1 , s 2 ) = µ 1 µ 2 (µ 1 -s 1 )(µ 2 -s 2 ) -ρs 1 s 2 , (4.4)
where µ 1 , µ 2 ≥ 0 and 0 ≤ ρ ≤ 1. We inject the bivariate Laplace transform expression (5.5) in the generating function defined in (3.11) to get

B(z 1 , z 2 ) = (1 + z 1 ) 1+r 1 (1 -z 2 ) 1-r 2 µ 1 µ 2 [z 1 (1/m 1 -µ 1 ) + µ 1 ] [z 2 (1/m 2 -µ 2 ) + µ 2 ] -z 1 z 2 ρ m 1 m 2 . (4.5) We set r 1 = r 2 = 1, m 1 = 1 µ 1 and m 2 = 1 µ 2 , so B(z 1 , z 2 ) = 1 1-z 1 z 2 ρ which implies that a k . = a k,l = ρ k δ kl . (4.6)
The double sum turns into a simple one, leaving us with the same form as in (4.1) and yielding the following result:

Proposition 2. The Downton Bivariate Exponential distribution belongs to Lancaster probabilities and its PDF admits an infinite series representation of the form

f ν (x, y) = +∞ k=0 a k Q ν 1 k (x)Q ν 2 k (y) f ν 1 (x)f ν 2 (y), (4.7)
where

a k = ρ k δ kl , f ν 1 is the PDF of a Γ 1 µ 1 , 1 distribution and f ν 2 is the PDF of a Γ 1 µ 2 , 1 distribution.
Historically, the DBVE distribution has been used to capture the joint lifetimes of two components for reliability purposes. These two components are assumed to collapse after a random numbers of shocks with exponential inter-arrival times. We therefore have

X Y = N +1 j=1 U 1j U 2j , (4.8)
where N is a counting random variable governed by a geometric distribution with parameters ρ. The {U ij } i∈{1,2},j∈N are mutually independent and exponentially distributed with parameter µ i and independent of N . We shed light here on a new way to construct Gamma-Gamma Lancaster probabilities through geometric compounding exponential random vectors.

Numerical illustrations

In this section, we aim at recovering the joint distribution of (S 1 , S 2 ) defined in (1.1). As the joint probability density and survival functions are not available in a closed form, we assess the accuracy of the polynomial approximation on well known bivariate distributions. Then, we propose an example of practical applications by performing computations of quantities of interest for risk management purposes. Note that we consider a case where no closed formula of the BPDF is available.

Approximation of common bivariate exponential distributions

We define the approximation error as the relative difference between the exact function and its approximation

∆f (x, y) = f (x, y) -f P roxy (x, y) f (x, y) , x, y ∈ R + , (5.1)
and we express it in percentages. We propose a two-dimensional graphical visualization. The color of the surface indicates the level of error at each point located by its coordinates (x, y). We compare the results of the polynomial approximation to the results of the moments-recovered method when approximating the joint PDF and the joint survival function of the considered distributions. The approximation of the joint PDF through the moment recovered method follows from

f α,α X,Y (x, y) = e -x-y Γ(α + 2)Γ(α + 2) Γ( αe -x + 1)Γ( α e -y + 1) × α-αe -x k=0 α -α e -y l=0 (-1) k+l f X,Y ( αe -x + k, α e -y + l) k!(α -αe -x -k)l!(α -α e -y -l) . (5.2)
This approximation tends towards the desired function when α and α tend to infinity. The approximation of the joint survival function is given by

F α,α X,Y (x, y) = αe -x k=0 α e -y l=0 α i=k α j=l α i i k α j j l (-1) i+j-k-l f X,Y (i, j).
(5.

3)

The approximation formulas (5.2) and ( 5.3) are derived in [START_REF] Mnatsakanov | Moment-recovered approximations of multivariate distributions: The Laplace transform inversion[END_REF]. It is recommended when using this method to evaluate (5.2) and ( 5.3) on a grid defined as

(x i , y j ) = ln α α -i + 1 , ln α α -j + 1 , i = 1, . . . , α, j = 1, . . . , α .
The final approximation follows from a simple interpolation procedure.

Expansion of the Downton Bivariate Exponential distribution

The Downton Bivariate Exponential distribution -DBV E(µ 1 , µ 2 , ρ), introduced in [START_REF] Dowton | Bivariate exponential distributions in reliabilty theory[END_REF], has been used to capture the joint lifetimes of two components for reliability purposes. These two components are assumed to collapse after a random numbers of shocks with exponential inter-arrival times. The joint PDF of a DBV E(µ 1 , µ 2 , ρ)-distributed random vector (X, Y ) is given by

f X,Y (x, y) = µ 1 µ 2 1 -ρ exp - µ 1 x + µ 2 y 1 -ρ I 0 2 √ ρµ 1 µ 2 xy 1 -ρ , x, y ∈ R + 2 (5.4)
where I 0 denotes the modified Bessel function of the first kind, µ 1 , µ 2 ≥ 0 and 0 ≤ ρ ≤ 1.

The bivariate Laplace transform of (X, Y ) is given by

f X,Y (s 1 , s 2 ) = µ 1 µ 2 (µ 1 -s 1 )(µ 2 -s 2 ) -ρs 1 s 2 . ( 5.5)
In this particular case, the study of the generating function of the coefficients is of interest to choose the parameters of the Gamma measures. We inject the bivariate Laplace transform expression (5.5) in the generating function defined in (3.11) to get

B(z 1 , z 2 ) = (1 -z 1 ) 1-r 1 (1 -z 2 ) 1-r 2 µ 1 µ 2 [z 1 (µ 1 -1/m 1 ) -µ 1 ] [z 2 (µ 2 -1/m 2 ) -µ 2 ] -z 1 z 2 ρ m 1 m 2 .
(

We set

r 1 = r 2 = 1, m 1 = 1 µ 1 and m 2 = 1 µ 2 , so B(z 1 , z 2 ) = 1 1-ρz 1 z 2 which implies that a k,l = ρ k δ kl , (5.7)
along with a good accuracy. For the numerical illustration, we set {µ 1 = 1/2, µ 2 = 2, ρ = 1/4}. Figure 1 shows the approximation error associated to the polynomial approximation.

An acceptable level of error is reached at almost every points of the bivariate probability density function with an order of truncation equal to 5. Figure 2 shows the approximation error of the moment-recovered based method. Increasing values of α and α are considered.

The approximation of the moment-recovered based method works quite well too. However, the polynomial expansion seems to be better suited to the problem due to the geometrical decay of the coefficients of the expansion. In order to provide a fair comparison, we perform an approximation of the Marshall-Olkin Bivariate Exponential distribution in the next subsection. The Marshall-Olkin Bivariate Exponential distribution -M OBV E(λ 1 , λ 2 , λ 12 ) has been introduced in [START_REF] Marshall | A multivariate exponential distribution[END_REF] with reliability interpretations. It aims at modelling the time-to-failure of a two components system that receives shocks impacting one or both of the components. Each shock can cause the failure of none, one or both of the components. The definition implies the possibility that the two components fail at the same time, which means that the distribution admits a continuous and a singular part. This type of distribution arises naturally within a two-dimensional framework. The Marshall-Olkin Bivariate Exponential distribution M OBV E(λ 1 , λ 2 , λ 12 ) has been originally characterized through its joint survival function

F X,Y (x, y) = exp [-λ 1 x -λ 2 y -λ 12 max(x, y)] . (5.8)
The BPDF associated with this distribution is divided into the sum of an absolutely continuous part with respect to the bivariate Lebesgue measure and another part being singular,

f X,Y (x, y) = λ 2 (λ 1 + λ 12 )F X,Y (x, y)1 x>y (x, y) (5.9)
+ λ 1 (λ 2 + λ 12 )F X,Y (x, y)1 x<y (x, y) (5.10)
+ λ 12 e -λ max(x,y) 1 x=y (x, y), (5.11)
where λ = λ 1 + λ 2 + λ 12 . The Laplace transform associated to the MOBVE distribution is

f X,Y (s 1 , s 2 ) = (λ -s 1 -s 2 )(λ 1 + λ 12 )(λ 2 + λ 12 ) + s 1 s 2 λ 12 (λ -s 1 -s 2 )(λ 1 + λ 12 -s 1 )(λ 2 + λ 12 -s 2 ) .
(5.12)

The distribution is not absolutely continuous with respect to the bivariate Lebesgue measure, the hypotheses of Proposition 1 are not satisfied. However, we manage to isolate the singular part. So we can expand the defective probability density function associated with the continuous part of the distribution and add the singular part to the polynomial expansion in order to recover the entire distribution. The continuous part has a joint defective PDF f C X,Y , which is the sum of (5.9) and (5.10), and admits a Laplace transform of the form

f C X,Y (s 1 , s 2 ) = s 2 (λ 1 + s 1 ) (λ 1 + λ 2 + s 1 )(λ 1 + λ 12 + λ 2 + s 1 + s 2 )
(5.13)

+ s 1 (λ 1 + s 2 ) (λ 1 + λ 12 + s 2 )(λ 1 + λ 12 + λ 2 + s 1 + s 2 )
.

(5.14)

The generating function of the coefficients of the expansion is tedious, making it difficult to choose relevantly the parameters of the expansion. Unlike in the case of the DBVE distribution, we cannot ensure that the chosen parametrization is the best. Approximations with different parametrization are compared in terms of accuracy.

• Parametrization 1:

{m 1 = 1 λ 1 +λ 12 , m 2 = 1 λ 2 +λ 12 , ν 1 = 1, ν 2 = 1} • Parametrization 2: {m 1 = 1 λ 1 , m 2 = 1 λ 2 , ν 1 = 1, ν 2 = 1} • Parametrization 3: {m 1 = 1 λ , m 2 = 1 λ , ν 1 = 1, ν 2 = 1}
We set {λ 1 = 1/2, λ 2 = 2, λ 12 = 1}. In view of a future probabilistic application, the error on the joint survival function is very interesting. Let us see how the polynomial method performs when it comes to approximating the joint survival function. We also apply the moment recovered method with the inversion formula dedicated to survival function approximation. On Figure 3, the relative difference between the exact value of the joint survival function and its polynomial approximation is plotted for different parametrizations and an order of truncation equal to 10. The first parametrization seems to be the most suited to the problem as the convergence towards the exact value is quicker. Figure 4 displays the relative error between the exact value of the joint survival function and its polynomial approximation with parametrizations 1 and different order of truncation.

On Figure 5, the difference between the exact value of the joint survival function and its moments approximation is plotted for increasing values of α and α . The results are satisfying for both methods, even if it looks like the moment recovered method does not perform well for x, y ∈ [0, 0.5]. One can note that the polynomial expansion with the first parametrization performs better than the moment-recovered based method. However, the moment-recovered based method has been enhanced in recent papers with the introduction of the scaled Laplace transform, see [START_REF] Mnatsakanov | A note on recovering the distribution from exponential moments[END_REF][START_REF] Mnatsakanov | Approximation of the ruin probability using the scaled Laplace transform inversion[END_REF]. These improvements have been made in the univariate case but might be extended soon to the bivariate case and comparison would be interesting in the future.

Approximation of the distribution of a bivariate claim amount distribution

The main concern of this paper is to deal with the distribution of .15) as parameters of the reference distribution, the generating function of the coefficients of the expansion defined in (3.11) takes the form

X 1 X 2 = N j=1 U 1j U 2j . ( 5
B(z 1 , z 2 ) = 1 1 + z 1 z 2 [p 2 -ρ(1 -p) 2 -p] .
(5.17)

The coefficients of the expansion are therefore a kl = p 2 -ρ(1 -p) 2 -p k δ kl . The geometric decay of the coefficients implies a fast convergence of the approximation towards the desired function. It is easily checked that

f U 1 ,U 2 1 2µ 1 (1 -p) , 1 2µ 2 (1 -p) < p,
making the chosen parametrization valid. The exact value of the survival function is unknown. A benchmark value is derived through Monte Carlo simulations using algorithms described in [START_REF] He | Efficient simulation of complete and censored samples from common bivariate distributions[END_REF]. In order to get a good proxy via Monte Carlo techniques, we produce 10 5 realizations of the couple (Y 1 , Y 2 ). Figure 6 shows the difference between the values of the survival function derived by Monte Carlo and by polynomial expansion.

One can see that we get close to the benchmark with an order of truncation equal to 5.

Figure 7 provides a 3D visualization of the joint probability density and survival function of the distribution.

Reinsurer's risk profile in presence of two correlated insurers

We now consider a reinsurer with two customers. We assume that the reinsurer accepted a stop-loss reinsurance treaty with limit b i in excess of priority c i from insurer i, where i = 1, 2. This means that if S i is the aggregated claims amount for insurer i, then the reinsurer must pay

Z = min [S 1 -c 1 ] + , b 1 + min [S 2 -c 2 ] + , b 2 .
For risk management or pricing purposes, the reinsurer is interested in computing P (Z > z) for different risk limits. To achieve that, one needs the joint distribution of (S 1 , S 2 ) and it is not enough to know the distribution of (S 1 + S 2 ), which would often be easier to obtain.

As an illustration, we assume that each insurer has to pay some specific claims amount, which correspond to the two independent components. We also take into account common shocks arising from events which cause claims for both insurers, like hail or storm episodes: for example, in the case of a hail episode, it might happen that 1 windshield has to be replaced for insurer 1 and that 2 windshields have to be replaced for insurer 2. Of course, for one particular event (the k th for instance) the two claims U 1,k (incurred for insurer 1) and U 2,k (incurred for insurer 2) are positively correlated, because they are influenced by the severity of the event causing the common shock. We therefore consider that the vector (S 1 , S 2 ) is described as

S 1 S 2 = N j=1 U 1j U 2j + N 1 i=1 V i N 2 i=1 W i , (5.18)
and has been defined in the introduction. The number of claims that incurs in the two portfolios is a counting random variable denoted by N , and governed by a negative binomial distribution N B(a, p). The severities of those claims are modeled through a sequence of i.i.d random vectors {(U 1,j , U 2,j )} j∈N , DBV E(µ 1 , µ 2 , ρ)-distributed, and independent

= a 2 = 1, p 1 = p 2 = 3/4, α 1 = α 2 = 1, β 1 = β 2 = 1}.
These sequences are mutually independent and independent from N 1 and N 2 . To recover the distribution of (S 1 , S 2 ) using the polynomial technology, we use a bivariate polynomial approximation to recover the distribution of (X 1 , X 2 ), and a univariate polynomial approximation to both of the components of the second part of the right-hand side of (5.18). The polynomial expansion in a univariate context has been presented in a previous work, we refer the reader to [START_REF] Goffard | A polynomial expansion to approximate the ultimate ruin probability in the compound Poisson ruin model[END_REF]. The joint probability density function of the second part is obtained by multiplying the marginal PDF due to the independence of the two components. The BPDF of (S 1 , S 2 ) is derived by convoluting the polynomial approximations of the different parts. One must pay attention to the singularities that arise from the non-zero probability of having {N = 0}, {N 1 = 0} and {N 2 = 0}. To expand the common part, we set m 1 = 1 µ 1 (1-p) , m 2 = 1 µ 2 (1-p) , r 1 = 1, r 2 = 1 as parameters of the polynomial expansion in view of the good results obtained in Subsection 5.2. The PDF of a geometric compound distribution with exponential claim sizes is available in a closed form, so we do not need to use univariate polynomial expansions. In order to appreciate the accuracy of our method, we need benchmark values. A sample of 10 5 Monte Carlo simulations of (S 1 , S 2 ) is produced. On Figure 8, the relative difference between the survival function obtained via Monte Carlo and polynomial expansion is plotted.

The level of error is still acceptable. Figure 9 displays the joint probability density and survival functions.

As there exists a symmetry between the two portfolios, we assume that the same reinsurance treaty is proposed to each insurer with the same priority and limit. We set {b Figure 10 displays two graphics. On the left-hand side, we plot the survival function of Z delivered by the two methods. We can see that the two curves overlap. On the right-hand side, we plot the relative difference between the two approximations to better appreciate their proximity. Some values of P (Z > z) are given in Table 1.

The results are quite promising. The polynomial expansion seems to be well-suited to carry out numerical calculations. The operational advantage of numerical methods over Monte Carlo techniques lies in the gain of computation time that they provide. It is difficult to quantify it as it varies given the software and the coding skills of the user. Note that, in the polynomial expansion case, the approximation of the PDF is almost immediate, although some computation time might be needed to derive survival function for instance.

Conclusion

We have considered polynomial expansions of bivariate densities with well defined Laplace transforms. We have focused our attention on densities on R + × R + and have considered Laguerre polynomials expansions. MOBVE and DBVE distributions have been studied and numerical illustrations lead to satisfactory results. An application to a bivariate aggregate claims amount distribution with dependence has been proposed. This work is a numerical step that ensures the efficiency of the bivariate polynomial approach. Interesting applications follow from the fact that the approximants admit a tractable form. The approximation of the BPDF can turn into a nonparametric estimation of the BPDF. Indeed, the coefficients of the expansion can be replaced with their empirical counterparts when data are available. The statistical extension of this numerical method will be at the center of a forthcoming paper.

ee

  sy f (y)dy. (3.21) ≥ 1 s [f (x)e sx -f (a)e sa ] . (3.22) Thus, we deduce that ∀x ≥ a, f (x) ≤ A(s)e -sx , where A(s) = s f (s) + f (a)e sa .As f (s, t) < +∞, we havef X,Y (0, t) = ty f X,Y (x, y) < +∞.(3.23)As y → e ty f X,Y (x, y) is a non-negative function, applying Fubini-Tonelli's theorem yields +∞ 0

Figure 1 :

 1 Figure 1: Error map for the joint PDF of a DBV E(1/2, 2, 1/4) distribution approximated by polynomial expansion

Figure 2 :

 2 Figure 2: Error map for the joint PDF of a DBV E(1/2, 2, 1/4) distribution approximated by the moment-recovered method

Figure 3 :

 3 Figure 3: Error map for the joint survival function of a M OBV E(1/2, 2, 1) distribution approximated by polynomial expansion with different sets of parameters.

Figure 4 :

 4 Figure 4: Error map for the joint survival function of a M OBV E(1/2, 2, 1) distribution approximated by polynomial expansion with different sets of parameters.

Figure 6 :Figure 7 :

 67 Figure 6: Discrete error map for the joint survival function of compound distribution [N B(1, 3/4), DBV E(1, 1, 1/4)]

  1 = b 2 = 4, c 1 = c 2 = 1} . Recall that the reinsurance cost is modeled by a random variable defined as Z = min [S 1 -c 1 ] + , b 1 + min [S 2 -c 2 ] + , b 2 . (5.19)

Figure 8 :Figure 9 :

 89 Figure 8: Error map for the joint survival function of (S 1 , S 2 ) obtained by polynomial expansion

Figure 10 :

 10 Figure 10: Survival function of the reinsurance cost obtained by polynomial expansion (blue) and Monte Carlo simulations (red) ; Relative difference between the survival function obtained with Monte Carlo and the one obtained by polynomial expansion

  The second part of the right-hand side of 5.18 represents the claim amounts specific to each insurer. The numbers of claims are modeled by two independent counting random variables N 1 and N 2 that are governed by two negative binomial distribution B(a i , p i ), where i = 1, 2. The sequences {V i } i∈N and {W i } i∈N are formed by i.i.d nonnegative random variables that represent the claim sizes. We assume that they are exponentially distributed Γ(1, β i ), where i = 1, 2. We set {a 1

Table 1 :

 1 Survival function of the reinsurance cost obtained by polynomial expansion and Monte Carlo simulations

	P(Z>z)

Acknowledgements

The authors would like to thank Pierre Miehe and Gilbert Macquart for their useful comments. The authors are also really greatful to an anonymous referee for his comments and remarks which improve greatly the quality of this paper. This work is partially funded by the AXA research fund and the research chair Actuariat Durable, sponsored by Milliman.

We assume here that the sequence of random vectors (U 1i , U 2i ) is i.i.d. and DBV E(µ 1 , µ 2 , ρ)distributed. The number of claims is governed by a negative binomial distribution N B(a, p).

In this case, the corresponding bivariate Laplace transform

)

In view of Corollary 1, the parameters must be chosen under the constraints

where the couple

We set {µ 1 = 1, µ 2 = 1, ρ = 1/4} as parameters of the DBVE distribution and {a = 1, p = 3/4} as parameters of the negative binomial distribution. In this particular case, if we set

, r 1 = 1, r 2 = 1