Erwan Salaun

Julien Apeloig

Frédéric Grisch

Charles-Edouard Yvonnet

Baptiste Nicolas

Frédéric Dionnet

Optical Investigation of Ignition Timing and Equivalence Ratio in Dual-Fuel CNG/Diesel Combustion

Dual-fuel engines are recognized as a short-medium term solution to reduce fuel consumption and pollutant emissions of CI engines, while maintaining high energy efficiency. Methane (CH 4 ) was chosen as it offers the best compromise between its heating value and H/C ratio. The high auto-ignition temperature of CH 4 requires auto-igniting a small quantity of liquid diesel before it initiates the combustion of the mixture. Therefore, new engine operations need to be specifically developed. This investigation explores the impact of time sequences of injection of the liquid fuel on the ignition of homogenous methane/air mixture. Experiments were performed on a Rapid Compression Expansion Machine (RCEM), to reproduce the operating and dynamic conditions encountered in a diesel engine cycle, allowing visualizations of fuel injection and combustion processes through a transparent piston. For the purpose of this work, the RCEM was modified to operate under dual-fuel conditions, while controlling the amount of diesel and methane-gas fueled. Experiments were performed for a wide range of equivalence ratios of the premixed charge. The study of the liquid fuel penetration and its consequence on igniting the homogenous charge was achieved using high-speed optical diagnostics. High-speed Schlieren technique (~22 kHz) was applied to characterize the diesel spray penetration, as well as the in-cylinder liquid fuel distribution. High-speed shadowgraphy and OH*-chemiluminescence techniques were used to determine ignitions delays. Moreover, the latest diagnostic was used to analyze the flame structure propagation and the heat release evolution.

Introduction

Over the last decades, the development of modern design of compression-ignition engines has largely evolved to enhance higher output efficiencies, greater specific-power outputs, increased reliability and reduction of pollutant emissions. There is an increased need, whether for environmental, economical, or resource preservation reasons, to operate natural or renewable energy sources. Compared to conventional diesel engines, which is characterized by a high efficiency but a high level of soot particles production, and compared to premixed charge gasoline engines, limited by low efficiencies induced by knocking constraints and pumping losses, gas-fueled diesel engines present higher efficiencies thanks to the reduction of pumping losses and heat transfer [START_REF] Karim | Dual fuel Diesel engines[END_REF][START_REF] Ando | Lean burn combustion engines[END_REF][START_REF] Soylu | Examination of combustion characteristics and phasing strategies of a natural gas HCCI engine[END_REF][START_REF] Takashima | Evaluation of Engine Performance and Combustion in Natural Gas Engine with Pre-Chamber Plug under Lean Burn Conditions[END_REF][START_REF] Weaver | Dual fuel natural gas/diesel engines: technology, performance, and emissions[END_REF][START_REF] Gibson | Comparison of propane and methane performance and emissions in a turbocharged direct injection dual fuel engine[END_REF][START_REF] Poonia | Experimental investigation of the factors affecting the performance of a LPGdiesel dual fuel engine[END_REF][START_REF] Papagiannakis | Experimental investigation concerning the effect of natural gas percentage on performance and emissions of a DI dual fuel diesel engine[END_REF][START_REF] Stewart | An experimental study of the dual-fuel performance of a small compression ignition diesel engine operating with three gaseous fuels[END_REF][START_REF] Alla | Effect of pilot fuel quantity on the performance of a dual fuel engine[END_REF]. Additionally, the use of low-cetane gaseous fuels such as Compressed Natural Gas (CNG) or methane offers attractive properties in regard to their higher octane number giving higher resistance to knock and to their short carbon chains allowing a systematic reduction in the particulate matter and the NOx emissions [START_REF] Karim | Dual fuel Diesel engines[END_REF][START_REF] Weaver | Dual fuel natural gas/diesel engines: technology, performance, and emissions[END_REF][START_REF] Gibson | Comparison of propane and methane performance and emissions in a turbocharged direct injection dual fuel engine[END_REF]. On the contrary, their use in diesel engines requires a high-energy source such as spark plug or pilot-fuel injection to achieve ignition. Previous studies reported in literature reveal that the second solution, defined as dual-fuel combustion operation mode with pilot injection of liquid fuel as igniter is particularly promising. Additionally, the technological conversion of a diesel engine operating in a dual-fuel combustion mode requires few modifications to existing hardware and offers the advantage to switch back to conventional diesel operating mode in case of shortfall in gas supply.

In a general manner, dual-fuel combustion can be considered as made of three main overlapping phenomena: injection [START_REF] Lata | Analysis of ignition delay period of a dual fuel engine with hydrogen and LPG as secondary Fuels[END_REF][START_REF] Carlucci | Combustion and emissions control in diesel-methane dual fuel engines: the effects of methane supply method combined with variable in-cylinder charge bulk motion[END_REF][START_REF] Gatts | An experimental investigation of incomplete combustion of gaseous fuels of a heavy-duty diesel engine supplemented with hydrogen and natural gas[END_REF][START_REF] Goldsworthy | Combustion behaviour of a heavy duty common rail marine diesel engine fumigated with propane[END_REF] and ignition of the pilot diesel-fuel [START_REF] Karim | Dual fuel Diesel engines[END_REF], ignition of the fuel-to-air mixture in the vicinity of the pilot spray and combustion of the remainder of the primary fuel-to-air mixture by flame propagation [START_REF] Kusaka | Combustion and Exhaust Gas Emission Characteristics of a Diesel Engine Dual-Fueled with Natural Gas[END_REF][START_REF] Badr | An Examination of the Flame Spread Limits in a Dual Fuel Engine[END_REF]. Efforts to interpret these complex phenomena need to be supported by a profound understanding of the relevant physical and chemical parameters in order to draw conclusions about the underlying coupled processes governing dual-fuel combustion. In this sense, various studies have previously been conducted to characterize combustion and exhaust emission of dual-fuel compression-ignition engines with diesel-fuel and CNG [START_REF] Dronniou | Optical Investigation of Dual-fuel CNG/Diesel Combustion Strategies to Reduce CO2 Emissions[END_REF][START_REF] Sahoo | Effect of Engine Parameters and Type of Gaseous Fuel on the Performance of Dual-Fuel Gas diesel Engines -A Critical Review[END_REF][START_REF] Serrano | Exploring the Potential of Dual Fuel diesel-CNG Combustion for Passenger Car Engine[END_REF][START_REF] Ishiyama | Improvement of Performance and Reduction of Exhaust Emissions by Pilot-Fuel-Injection Control in a Lean-Burning Natural-Gas Dual fuel Engine[END_REF][START_REF] Königsson | Characterization and Potential of Dual Fuel Combustion in a Modern Diesel Engine[END_REF]. A notable problematic feature of dualfuel combustion is its performances trade-off, especially when using small amounts of pilot fuel. It is manifested by a reduction of the NOx-particulates with no change in the total efficiency but it is revealed by an increase of specific energy consumption, cyclic variations, and levels of unburned hydrocarbons and CO relative to the diesel engine operation. These limitations arise from the fact that the flame cannot propagate fast and far enough within the time available to consume the gaseous-fuel/air mixture [START_REF] Karim | Dual fuel Diesel engines[END_REF][START_REF] Elliot | Dual fuel combustion in diesel engines[END_REF][START_REF] Ding | Emissions and fuel economy of a prechamber diesel engine with natural gas dual fueling[END_REF]. Such effects can be attenuated by increasing the gaseous-fuel concentration to allow flame propagation into the combustion chamber, or by taking other actions such as: optimizing the stratification of gaseous fuel admission, heating of the intake charge, optimizing the pilot injection timing, changing the pilot fuel, etc. [START_REF] Karim | An examination of some measures for improving the performance of gas fueled diesel engines at light load[END_REF]. To be successful, such implementation requires a detailed knowledge of the sensitivity to these parameters for which dual-fuel operation begins to be favorable. Most of the mentioned experimental investigations were generally performed on single cylinder, direct injection, naturally aspirated research engines, which provide better flexibility of control and economy of testing. On the other hand, fundamental investigations concerning enhanced ignition with a pilot fuel system on gas engines are sparse in the open literature. For instance, few experiments were performed on RCEM to investigate the dual-fuel combustion characteristics of lean natural gas/air mixtures using diesel fuel as ignition source [START_REF] Egúsquiza | Experimental Investigation of the Natural Gas / Diesel Dual-Fuel Combustion Using a Rapid Compression Machine[END_REF][START_REF] Kikusato | A Study on the Characteristics of Natural Gas Combustion at a High Compression Ratio by Using a Rapid Compression and Expansion Machine[END_REF]. Moreover, these studies were generally performed in low-load operating conditions, avoiding then to draw conclusions on the benefits of dual-fuel combustion when applied in real diesel combustion engines.

7/20/2015

Based on these observations, the objective of the present work aims at exploring the impact of the time sequence of the liquid-fuel injection on auto-ignition of the homogenous charge of methane/air and its impact on the flame propagation in operating conditions similar to diesel engines at partial-and full-load operating conditions. Experiments were performed on an "optical" RCEM, providing a valuable reproduction of various heterogeneous combustion modes and especially operating conditions encountered in diesel engine cycles at partial-load, between -40 and +20 CAD around Top Dead Center (TDC). The RCEM was modified to operate under dual-fuel conditions, while controlling the amount of either diesel or methanegas fuels. The primary fuel is methane, which is ignited by the microinjection of liquid-diesel fuel. The cylinder pressure condition selected in the current study ranges between 90 and 120 bar, corresponding to diesel engine-like conditions at partial-loads. The measurements consisted in studying the effects of the methane/air equivalence ratio, pressure and temperature operating conditions on autoignition, inflammation and flame propagation inside the combustion chamber. To this purpose, various optical measurement techniques were systematically employed as measuring tools. In particular, the study of the pilot liquid-fuel penetration and its consequence on ignition of the premixed charge was accomplished using high-speed optical diagnostics. High-speed Schlieren diagnostic was applied to characterize the diesel spray tip penetration as well as the in-cylinder liquid fuel distribution on dual-fuel ignition. Additionally, high-speed shadowgraphy and OH*chemiluminescence techniques were used to determine ignition delays. Complementary high-speed OH*-chemiluminescence experiments were finally carried out to investigate the effect of fuel injection on the timing sequence of the combustion propagation inside the combustion chamber.

Experimental Apparatus

The following presentation of the experimental apparatus focuses on two aspects: the RCEM, its operating principle and modes, and the optical diagnostics, their setups and characteristics.

Rapid Compression Expansion Machine

RCEM Operating Principle

The operating principle of the RCEM (Figure 1.a) can be separated in two parts presented in Figure 1.b; the experiment area (on the left) and the driving area (on the right). The experiment area consists of the liquid-fuelled injector and the combustion chamber enclosed by the cylinder wall and the transparent piston head. This design permits to successfully reproduce the operating conditions of a single cylinder diesel engine [START_REF] Egúsquiza | Experimental Investigation of the Natural Gas / Diesel Dual-Fuel Combustion Using a Rapid Compression Machine[END_REF][START_REF] Kikusato | A Study on the Characteristics of Natural Gas Combustion at a High Compression Ratio by Using a Rapid Compression and Expansion Machine[END_REF][START_REF] Schlatter | Experimental Study of Ignition and Combustion Characteristics of a Diesel Pilot Spray in a Lean Premixed Methane/Air Charge using a Rapid Compression Expansion Machine[END_REF]]. The RCEM sizing specifications are summarized in Table 1. The machine is instrumented to control and adjust the liquid-fuel injection timing and volume, the cylinder pressure, and the piston stroke. The latest is controlled thanks to the driving area. It uses hydraulic and pneumatic chambers coupled with solenoid valves to provide the energy to move the piston. This mode of operation means that the centre of gravity is kept at every position of the piston, guaranteeing a strong acceleration of the piston without any exterior forces and therefore minimizing vibrations of the entire machine. Thus, vibration free operation of the machine favours the use of optical measurement diagnostics with a high reproducibility ratio. The RCEM also offers the opportunity to reproduce most different engine states and thermodynamic conditions encountered in conventional diesel engines. Due to its design specially devoted to perform single-shot runs, the RCEM is efficient to suppress the cycle-to-cycle variations of real engines [START_REF] Egúsquiza | Experimental Investigation of the Natural Gas / Diesel Dual-Fuel Combustion Using a Rapid Compression Machine[END_REF][START_REF] Kikusato | A Study on the Characteristics of Natural Gas Combustion at a High Compression Ratio by Using a Rapid Compression and Expansion Machine[END_REF][START_REF] Schlatter | Experimental Study of Ignition and Combustion Characteristics of a Diesel Pilot Spray in a Lean Premixed Methane/Air Charge using a Rapid Compression Expansion Machine[END_REF]. The single-shot mode enables the simulation of full-load operating conditions. Therefore, it can dissociate the fuel evaporation process from the aerodynamics inside the combustion chamber. This is of prime importance to provide a detailed understanding of the physical and chemical processes occurring before and during dual-fuel combustion. The optical access is achieved through an optical window in the bottom of the piston with a diameter, D visu , of 45 mm and 16.6 mm in depth coupled with a 45° tilted mirror (Figure 2). An additional transparent ring can be placed as an optical cylinder liner to allow the penetration of a laser sheet inside the combustion chamber. The materials of the windows are synthetic quartz glass. Thanks to the single-shot operating mode, the RCEM is more efficient than a "classical" visualization high-pressure test cell to observe and analyze the basic phenomena involved. This is true since the RCEM is able to reproduce the dynamic motion of the piston with a 95% confidence, as shown in Figure 3. The valid region (-40 -+20 CAD relative to the TDC) is obtained by comparing the piston stroke of a conventional diesel engine and the one of the RCEM. These results are shown for two operating conditions, 1300 and 2800 rpm, which produce compression ratios, , of 10 and 25, respectively. Unlike an optical internal combustion engine, the RCEM allows a maximum pressure around 250 bar, permitting high compression ratios, as well as simulating full-load engine conditions.

7/20/2015 Finally, both piston and cylinder walls can be heated up to 493K (120°C) to be representative of real engine conditions. Being originally designed to only test conventional-like injection diesel engines, the RCEM was modified for the current study.

RCEM Modifications for Dual-Fuel Combustion

The main adaptation the RCEM required to perform dual-fuel combustion took place on the intake air/swirl pipe. It was adapted for to inject the gaseous CH 4 inside the combustion chamber. To guarantee perfect mixing between air and CH 4 , the gaseous fuel is injected straight into the intake pipe that feeds the compression chamber. This modification made necessary to synchronize the dualfuel injection with the operation of the RCEM. Methane is injected to the combustion chamber by compressed air that supplies the charge gas pressure of the machine. This process assures complete mixing and transportation to the RCEM combustion chamber before the test cycle. The gas is supplied by a high-pressure bottle equipped with a pressure regulator. A plenum is used to feed the CH 4 injector to the desired pressure. This system is regulated by a differential pressure regulator system, which combines a Bronkhorst EL-PRESS ® pressure meter/controller and a Bronkhorst EL-FLOW ® mass flow meter/controller. It supplies the injector with regulated constant pressure above the intake pressure (variable during the intake), as represented in Figure 4.

Intake and exhaust are controlled with respect to the piston stroke, outside of the valid region period. The valves are opened and closed using a specifically developed routine, which ensures that the valves are closed during the experiment. Concerning the exhaust, the valve is open when the piston is at the BDC. Then, air is injected as a purging procedure in order to clean the combustion chamber and reset it for another test.

Optical Diagnostics

As described earlier, the RCEM allows the use of optical measurements techniques [START_REF] Zhao | Laser Diagnostics and Optical Measurement Techniques in Internal Combustion Engines[END_REF] to enhance our understanding of the effect of the Start Of Injection (SOI), and of the equivalence ratio of the charge, φ, on dual-fuel combustion. It is found that for the diagnosis of the injection and combustion phenomena, a timeresolved sequence of images is of great interest, especially if strong time variation is observed. For this reason, three methods which deliver image sequences at high-repetition rates are considered.

First of all, the Schlieren imaging technique was used to study the time evolution of the pilot liquid fuel injection inside the combustion chamber. Second, the shadowgraph imaging diagnostic was used to analyze the ignition sequence. Third and last, the OH*chemiluminescence imaging technique was implemented to study the ignition sequence and the flame behavior after ignition. This technique offers a comparative and complementary analysis of the phenomena with respect to the widely-used shadowgraph imaging technique.

Schlieren Imaging for Diesel Spray Characterization

The diagnostic used to visualize the time evolution of the liquid and gaseous diesel fuel penetration inside the combustion chamber was a high-frequency Schlieren system allowing the acquisition of images at a frame rate up to 100 kHz. The setup shown in Figure 5 uses an Nd-Yag high-speed laser source emitting at a wavelength of 532 nm. The laser beam is routed to the combustion chamber by an optical system specially adapted for this RCEM. The light is reflected on the cylinder head bottom, mirrorpolished for this purpose. This reflected light is passed through a narrow bandpass filter centered at 532 nm (±20 nm), and collected with a digital high-speed Phantom v7 video system with an acquisition frequency of 22.1 kHz and a resolution of 288x288 pix².

The bandpass-filter is used to limit the noise from the flame emission and therefore to protect the camera sensor. The laser, camera and injection system are triggered and synchronized with the RCEM microcontroller and delay pulse generators.

7/20/2015 Figure 6 illustrates a typical sequence of Schlieren images highlighting the fuel penetration inside the combustion chamber. The white circle at the center of the images is the injector, from which the liquid diesel is exiting. Diesel in its liquid phase is represented in black while the fuel vapor is shown in white.

Shadowgraph Imaging for Ignition Delay Identification

The shadowgraph imaging technique is an optical diagnostic similar to the Schlieren imaging diagnostic except that it cannot discriminate between liquid and vapor fuel. This diagnostic was used to visualize the time sequence of the autoignition, and therefore to obtain the ignition delay since the start of the pilot diesel-fuel micro-injection.

For the experiment, high-speed shadowgraph imaging technique was used. A Spectra-Physics Hg-Xe high power light source (500 W) was combined with the same digital high-speed Phantom v7 video system and was operated in the same conditions as for the Schlieren imaging technique (Figure 7). For instance, typical shadowgraph images recorded during a run are presented in Figure 8. On the latest, the diesel fuel, in liquid and vapor phases, is shown in black, while the white parts are representative of the combustion zone containing: radicals, combustion products, and soot. Hence, the ignition delay analysis is based on the time difference between the start of injection and the first image in which a white spot, characteristics of ignition, is observed. This diagnostic was only restricted to visualize the autoignition and not the flame propagation because the detection system collects a broadband-light emission induced not only by flame front radicals (like OH*) but by all the chemical species and particulates transported in the exhaust gases. It is important to note that the two aforementioned imaging techniques collect the full light spectrum. Another important observation is the dark spot on the right hand side displayed on the images in Figures 6 and8, which is only attributed to the presence of the pressure sensor and not to the phenomena observed.

OH* Chemiluminescence Imaging for Flame Structure and Propagation Analysis

In hydrocarbon combustion, the OH*, an excited radical, is produced in large quantities in the flame front. This naturally excited radical is largely used in the literature to study the position and the structure of flame in various configurations.

Results and Discussion

The following section is dedicated to the presentation of the methodology used to evaluate the quantity of pilot fuel which will be injected into the combustion chamber of the RCEM during dual-fuel combustion experiments. Thereafter, results on the phenomena observed during dual-fuel experiments such as fuel injection, spray visualization, ignition delay, flame structure, and flame propagation are presented and discussed.

Preliminary study

To prepare the dual-fuel experiments, investigations were performed to evaluate the minimal quantity of diesel fuel required to solely generate autoignition sites and reduce as much as possible its contribution the global heat released during combustion process.

Definition of Injection

The injection device selected for dual-fuel experiments is chosen among the devices available on the market and present worldwide on diesel combustion engines. The RCEM is equipped with a 6-hole Bosch common-rail injector placed at the center of the cylinder head.

The fuel pump pressure can be varied up to 1600 bar.

This study is centered on apprehending the overall phenomena rather than focusing on characterizing the injection, or the best performances possible. Therefore, resulting data analysis gives values averaged over the 6 fuel sprays.

To perform experiments in the RCEM in engine-like conditions, the pilot diesel-fuel is injected inside the combustion chamber with a reduced injection pressure and injection duration compared to the ones used in conventional diesel engine operations. To inject the smallest amount of pilot fuel, it is necessary to adapt the drive current pulse. In the present study, the high-pressure fuel injector is actuated with a specific current curve in order to comply with the requirement for defined and reproducible fuel-injection processes.

The conventional Engine Control Unit (ECU) is replaced by a power driving module (IPoD Coil Power Driver EFS 8427) on the RCEM. It uses magneto-electric technology to deliver the booster voltage required to bring the current up to a high value as quickly as possible in the booster phase. The opening current is applied for 350 µs, followed by a 50-µs holding current. In parallel, the injection pressure is limited to 400 bar. After calibration experiments, the diesel injected quantity is measured around 1.4 mm 3 . This volume is in accordance with the manufacturer data and published data [START_REF] Dronniou | Optical Investigation of Dual-fuel CNG/Diesel Combustion Strategies to Reduce CO2 Emissions[END_REF].

To evaluate the reproducibility of this diesel-fuel micro-injection, high-speed Schlieren visualizations were performed. The microinjection sequence was performed five times with the cylinder only charged with air. For this preliminary study, the start of injection was fixed at -30 CAD relative to the TDC. The penetration lengths presented in Figure 10 were measured from these tests. These results ensured a 95% confidence on the reproducibility of liquid fuel injection. The same confidence was observed by analyzing the ignition delays from the same data acquisitions. Performances of the RCEM are compared between a compression cycle of the charge only with no diesel injection, and a cycle including only the diesel micro-injection. The goal of this comparison is to ensure that the diesel micro-injection is only producing the energy required to ignite the charge, and nothing more, in order to maximize the diesel substitution rate. In the same time, the chargeonly case validates the need of diesel injection for its ignition. In Figure 11 are reported two parameters of interest, the measured cylinder pressure (lines) and the heat release (dashed lines), which is calculated from the cylinder pressure, temperature, and piston displacement data. This

Selection of test cases

Multiple preliminary tests were performed to verify and ensure the entire apparatus operability range and safety of the RCEM, plus to 7/20/2015 select the test cases to be performed during the dual-fuel experiments. Table 2 summarizes the parameters and variables selected for the RCEM, air, pilot fuel and methane fuel. As shown in Table 2, the equivalence ratio of the methane/air mixture, which is known to be a key parameter in terms of impact on combustion efficiency, was ranged between 0.8 and 1.2. The lower limit of 0.8 is explained by the appearance of knocking possibly linked to the high operating pressure chosen in this study. Significant knocking was observed for lower equivalence ratio, generating heavy stress on the optical RCEM. The upper limit of 1.2 was chosen to provide a global tendency of the equivalence ratio effect around the stoichiometric value. The end goal of the study is to provide an enhanced understanding of the phenomenological processes and not to optimize the engine performances. Table 3 displays the operating conditions explored in this study, and the associated diagnostics. Additionally, the oxygen and the methane molar fraction of the charge are presented for each equivalence ratio. Moreover, the cylinder pressure, P cyl , and the corresponding isentropic-based computation of the temperature, T isen , are given at the start of the diesel fuel injection (SOI).

Spray Visualization

The results presented in this section are deduced from the Schlieren images. The data processing of these images were performed with a home-made routine. The first outcome presented in Figures 12 and13 is the measurement of pilot liquid-fuel penetration length. From the analysis of the luminosity of the spray images, the length of the visible part of the spray was obtained. This length, defined as spray penetration, was assumed to be equal to the distance from the nozzle to the furthest point on the spray axis where the luminosity was above a given threshold. It was assumed that the spray luminosity varies smoothly along the spray, and the effect of cluster shedding from the spray tip was ignored.

The effects of the SOI for an equivalence ratio of 0.8, and of the equivalence ratio for a SOI at -30 CAD, are displayed, respectively, in Figures 12 and13. The lengths of the liquid-fuel penetration into the cylinder are the average values over the six fuel-sprays and over five runs performed for each test case. These averaged values are also fitted using second-order polynomials to show the time evolution tendencies. The SOI characterizes the cylinder conditions (pressure and temperature) at the start of the pilot fuel injection. An observation of the results of Figure 12 shows a decrease of the penetration length as the SOI is closer to the TDC. Figure 13 presents a similar result as the equivalence ratio is decreasing. In this case, the main change in the cylinder is the charge density, which is increasing as φ is decreasing. Both results are corroborating this effect as the cylinder charge pressure, temperature, and density are increasing as the SOI is closer to the TDC. This behavior is consistent with previous studies on liquid-fuel penetration in diesel engine [START_REF] Espey | The Effect of TDC Temperature and Density on the Liquid-Phase Fuel Penetration in a D. I. Diesel Engine*[END_REF][START_REF] Siebers | Scaling Liquid-Phase Fuel Penetration in Diesel Sprays Based on Mixing-Limited Vaporization[END_REF][START_REF] Pastor | Fuel effect on the liquid-phase penetration of an evaporating spray under transient diesel-like conditions[END_REF].

Ignition Delay

The shadowgraph imaging and the OH*-chemiluminescence imaging techniques were used to measure the ignition delay. Ignition delay is defined as the time interval between the start of injection and the autoignition. In the experiments, the ignition delay was measured based on a method assessing the light intensity variation. Both measurements techniques were used in order to assess their 7/20/2015 comparable nature. The interest of this study focuses only on the ignition delay between the triggered pilot fuel injection and autoignition, for both techniques. Each result is based on a five-run dataset presenting no variation whatsoever in the number of frames between injection and ignition. The results for both techniques are presented in Figure 14. These results uncertainties are respectively ±45 µs and ±67 µs for the shadowgraph and the OH*chemiluminescence imaging techniques. These uncertainties come from the delay between two successive images. The representation in bar graphs illustrates both the effect of the fuel injection timing and the equivalence ratio on the variation of ignition delay, obtained from the shadowgraph (left) and the OH*-chemiluminescence (right).

First of all, the OH*-chemiluminescence experiments were not systematically performed for all the twelve test cases because the eight test cases presented in Figure 14 shows excellent agreement with measurements performed with the shadowgraph imaging technique. The small differences between ignition delays measured with both techniques arise from the different repetition rates of images acquisition; 22.1 kHz for shadowgraph imaging against 15 kHz for OH* chemiluminescence imaging. The linked bias is around 20 µs. The results comparison allowed to assess that the OH*chemiluminescence technique gives comparable results to the ones from the shadowgraph imaging. Therefore, future work could be simplified to implement only the Schlieren and the OH*chemiluminescence imaging techniques instead of all three. However, this is true if and only if there is no need for precise measurements of this ignition delay. Shadowgraph imaging technique will always be more precise as there is no intensifier required. When focusing on analyzing the results for a specific SOI, the emerging observation is that for increasing premixed equivalence ratio, the ignition delay increases as well. Part of the explanation comes from the variation of the oxygen molar concentration,  O2 (cf. Table 3). The increase in oxygen concentration improves the diesel dilution, and therefore decreases the local equivalence ratio toward a faster igniting point.

It is interesting to note that the less oxygen concentration the higher the pressure must be to generate an ignition of the dual-fuel mixture. Finally, focusing on the effect of the injection timing, and therefore looking at a single equivalence ratio, it is observed that the ignition delay decreases as the SOI is closer to the TDC. This result can be correlated to the temperature and pressure increase in the combustion chamber.

Flame Propagation

The goal of this section is to contribute to a better understanding of the dual-fuel combustion. The main idea is to follow the flame propagation of the burning methane/air mixture, after autoignition by the pilot fuel. Results from the OH*-chemiluminescence imaging are presented hereafter. In order to visualize the injection and ignition processes along the simulated engine cycle, Figure 15 regroups the cylinder pressure of the dual-fuel cases with the SOI = -30 CAD, and φ = 0, 0.8, 1.0, and 1.2. In addition, the start of injection and the beginning of the various ignitions are shown. The timeline is centered at the TDC.

The next figures [START_REF] Badr | An Examination of the Flame Spread Limits in a Dual Fuel Engine[END_REF][START_REF] Dronniou | Optical Investigation of Dual-fuel CNG/Diesel Combustion Strategies to Reduce CO2 Emissions[END_REF] show sequences of OH*-chemiluminescence images with a time interval of 200 µs. This time interval corresponds to the third of the acquisition frequency (15 kHz), meaning that only one out of three instantaneous OH* images is presented. These sequences display the first 2.4 ms following the ignition of the premixed charge, which is taken as the time reference. For example, the zero-times, presented Figure 17, correspond to the ignitions shown in Figure 15 for the various premixed equivalence ratio. As a reminder, the OH* radical is significantly spawned during combustion and is one of the most relevant chemical tracers of the flame front. Hence, the intensity of the OH* emission is considered proportional to the local heat release produced by combustion. The following results show data with the same color-scale, allowing a direct comparison of the heat release between the different test cases. Finally, it is important to notice the visual illusion produced by the displacement of the diopter, linked to the piston motion. This has the effect of showing a movement of the flame center, due to the position of the tilted 45°-mirror.

Figure 16 shows the comparison of the spatial OH* chemiluminescence evolution with respect to the SOI for the stoichiometric equivalence ratio, while Figure 17 depicts the effect of the methane/air equivalence ratio on the OH* distribution for SOI at -30 CAD.

As an observation of results displayed in Figure 16, it is possible to note that combustion is more efficient with SOI closer to the TDC.

7/20/2015

The areas of high combustion intensity around the pilot sprays are much wider when SOI is moved closer to the TDC. This effect can be attributed to the cumulative effects of the reduced ignition delay with SOI and larger values of pressure and temperature inside the cylinder. A shorter ignition delay reduces the accumulation of large amounts of evaporated diesel fuel, favoring the sequential ignition of the diesel fuel and methane that is entrained by the pilot spray in the close vicinity of the spray combustion. At the same time, the combined increase of pressure and temperature inside the cylinder leads to an increase in reaction rate, hence leading to a better combustion. Despite favorable conditions, no flame propagation in the premixed charge can be sustained, probably due to low levels of turbulence inside the cylinder which prevent a good entrainment of the hot product gases towards the unburned methane/air mixture. As shown in Figure 17, the reduction of equivalence ratio at a fixed condition of pressure and temperature also increases the combustion efficiency.

As previously observed, the reduction of the ignition delay promotes the fast ignition of diesel fuel and methane surrounding the pilot sprays. Furthermore, the penetration of the fuel sprays in a premixed charge with lower density seems to enhance the inflammation phase as well as the subsequent premixed combustion, resulting in higher OH* signals over the volume of the combustion chamber. This effect can be attributed to the existence of higher turbulence levels which widely diffuse the hot products of combustion towards the unburned methane/air mixture. It seems also plausible that the flame propagation towards the unburned premixed charge will be dependent of the mixture stratification and in particular the local equivalence ratio distribution around the fuel sprays. To gain additional information from the evolution of the flame in the visible part of the combustion chamber, all pixels of each OH*chemiluminescence image were summed to obtain global sense of the combustion evolution. For instance, the influence of methane/air equivalence ratio is shown in Figure 18. The OH* signals exhibit reverse behavior; the less there is methane in the cylinder, the higher the OH* intensity is. This result suggests a strong influence of the premixed phase on the OH* emission. For high equivalence ratio (φ = 1.2), the global OH* signal shows a one-peak timing evolution spreading on several milliseconds. This evolution is mainly attributed to the heat released by combustion of the six fuel sprays. On the contrary, higher and longer OH* signal with a two-peak evolution is observed with reducing methane content, demonstrating a longer and better combustion. As shown in Figure 18, this phenomenon is exacerbated for the lowest equivalence ratio investigated, and a net two-peak evolution is clearly visible. The existence of a second peak arises from the development of combustion reactions controlled by the mixing of hot products and the remaining unburned methane/air mixture (visible on Figure 17, .φ = 0.8, at 1.6-1.8 ms).

Conclusions

An extended experimental investigation was conducted to obtain a deeper understanding of ignition, inflammation and combustion processes occurring in dual-fuel CNG/diesel combustion systems. The dual-fuel combustion concept consists of burning a homogeneous methane/air charge ignited by a pilot micro-injection of liquid diesel fuel. The operating pressure conditions selected for the current study are ranging between 90 and 120 bar, corresponding to cases of diesel engine at partial-loads. For this purpose, a Rapid Compression Expansion Machine (RCEM) with optical access through the moving piston was modified to operate in the defined conditions without changing its main configurations. Measurements were taken with various advanced high-speed optical diagnostics including Schlieren, shadowgraphy and OH*-chemiluminescence techniques. The measurements consisted of studying the effects of the methane/air equivalence ratio, pressure, and temperature operating conditions on the auto-ignition, inflammation, and flame propagation of the methane/air gaseous charge inside the combustion chamber.

According to the results of this research, the following conclusions can be outlined.

7/20/2015

The short duration adopted to inject the small amount of diesel fuel allows the methane/air charge to ignite and to propagate flames in the combustion chamber for all operating conditions tested. From the OH*chemiluminescence and shadowgraph imaging, increasing autoignition delay of the pilot diesel spray is observed with decreasing temperature and pressure, independently of the methane/air equivalence ratio. Furthermore, the auto-ignition delay is also significantly reduced in time when the equivalence ratio is reduced.

The work also revealed that the penetration length of the diesel spray tends to be higher when the equivalence ratio is increasing. The same effect is also observed when the SOI is moving away from the TDC. The smaller the gas density, the longer is the penetration length.

These effects (auto-ignition and penetration length) associated with the use of a 6-hole injector produce multi-ignition sites on a large surface area of the combustion chamber, which enhance flame propagation. However, an observation of the time evolution of the spatial OH* chemiluminescence distribution shows that flame propagation is a function of the gaseous charge. For high-equivalence ratios ( =1 and 1.2), the flame progresses only near the pilot diesel spray region, limiting the consumption of the gaseous methane/air charge. By contrast, the cumulative effect of reducing the ignition delay and decreasing the spray penetration in lower equivalence ratio cases enhances autoignition and flame propagation inside the cylinder. This, along with increasing turbulence generated by the spray penetrating the gaseous charge of lower density, improves mixing of the liquid diesel and the gaseous charge, hence, improving the flame propagation in the full section of the combustion chamber.

Increasing the level of turbulence such as the ones encountered in classical combustion engines could then provide major benefits for utilizing the CNG/diesel dual-fuel operating mode.

It is clear that additional work is required to develop the dual-fuel concept, particularly for full-load charge applications. The results presented in this study do not cover all the underlying physical processes governing the dual-fuel combustion. However, they provide a good understanding of the auto-ignition of methane/air mixture by a pilot micro-injection of diesel fuel. Specific experiments need to be conducted in engine-like operating conditions in which local turbulence levels will be higher than the ones used in this experiment. The goal would be to validate the idea that the turbulence will permit the complete burning of the CH 4 /air charge. Another important component of these experiments would be to measure pollutants (in situ and ex situ) in order to correlate the results of ignition and flame propagation on the pollutant emissions.

Figure 1 .

 1 Figure 1. Image (a) and schematic (b) of the RCEM.

Figure 2 .

 2 Figure 2. RCEM configuration for optical diagnostics.

Figure 3 .

 3 Figure 3. Diesel engine vs. RCEM piston strokes.

Figure 4 .

 4 Figure 4. RCEM intake modifications for dual-fuel combustion experiments.

Figure 5 .

 5 Figure 5. RCEM optical configuration for Schlieren imaging.

Figure 6 .

 6 Figure 6. 22.1 kHz diesel injection Schlieren imaging sequence (t ≈ 45 µs); SOI = -30 CAD & φ = 0.8; Dvisu = 45 mm.

Figure 7 .

 7 Figure 7. Shadowgraph imaging setup.

Figure 8 .

 8 Figure 8. Sample high-speed shadowgraph imaging sequence (t ≈ 45 µs); SOI = -20 CAD & φ = 0.8; Dvisu = 45 mm.

Figure 9 .

 9 Figure 9. High-speed OH*-chemiluminescence imaging setup. A high-speed CMOS camera (LaVision HSS8 / Photron Fastcam SA5) coupled with an external image intensifier (Lavision HS-IRO) is used to record the chemiluminescence of OH*. The broadband flame emission collected through the piston window is then filtered through a narrow bandpass-filter (Semrock BrightLine FF01-315/15-

Figure 10 .

 10 Figure 10. Time evolution of the penetration length of the diesel microinjection in air; SOI = -30 CAD & φ = 0.

  figure compares three different cycles: a charge only, in blue (no SOI & φ = 1.0), a diesel micro-injection only, in red (SOI = -30 CAD & φ = 0), and a dual-fuel, in black (SOI = -30 CAD & φ = 1.0). The result indicates a limited effect of the diesel micro-injection on the RCEM performances, but reveals its importance to ignite the charge.

Figure 11 .

 11 Figure 11. Comparison of measured cylinder pressure (lines) and computed heat release (dashed lines) for three cycles: charge only -blue (no SOI & φ = 1.0), diesel micro-injection only -red (SOI = -30 CAD & φ = 0), and dual-fuel -black (SOI = -30 CAD & φ = 1.0).

Figure 12 .

 12 Figure 12. Comparison of the effect of the start of injection on the liquid diesel penetration length analyzed from Schlieren imaging; φ = 0.8.

Figure 13 .

 13 Figure 13. Comparison of the CH4/air equivalence ratio on the liquid diesel penetration length analyzed from Schlieren imaging; SOI = -30 CAD.

Figure 14 .

 14 Figure 14. Comparison of ignition delays as a function of the SOI, for various φ , obtained from shadowgraph imaging (left) and OH* chemiluminescence (right).

Figure 15 .

 15 Figure 15. Pressure cycle comparison with respect to the equivalence ratio; SOI = -30 CAD & φ = 0, 0.8, 1.0, and 1.2.

Figure 16 .Figure 17 .

 1617 Figure 16. Comparison of SOI effect on OH*-chemiluminescence evolution; SOI = -40, -30, and -20 CAD & φ = 1.0; Dvisu = 45 mm.

Figure 18 .

 18 Figure 18. Post-ignition evolution of global OH*-chemiluminescence for the cases presented in Figure 17; SOI = -30 CAD & φ = 0.8, 1.0, and 1.2.

Table 1 .

 1 RCEM specifications.

	Cylinder Volume at BDC	[cm 3 ]	886
	Cylinder Volume at TDC	[cm 3 ]	27
	Piston Diameter	[mm]	84
	Piston Stroke	[mm]	160
	Combustion Chamber Volume	[cm 3 ]	913

Table 2 .

 2 RCEM, air, and fuel parameters and variables used during the experiments realized for this study.

	RCEM Parameters
	Compression Ratio		16:1
	Wall Heating Temperature	[K]	363
	Air Parameters
	Admission Pressure	[bar]	2
	Admission Intake Temperature	[K]	293
	Fuel Parameters
	Diesel Injection Pressure	[bar]	400
	Diesel Injection Time	[µs]	400
	Diesel Injection Volume	[mm 3 ]	1.4
	CH4 Injection Pressure	[bar]	2
	Fuel Variables
	Start Of Injection	[CAD]

-40, -30, -20 CH4/Air Equivalence Ratio 0, 0.8, 1.0, 1.2

Table 3 .

 3

			CH4/Air Equivalence Ratio	Pcyl [bar]	Tisen [K]
			0	0.8	1.0	1.2	@ SOI	@ SOI
		-40	* X	*	* X	*	20	565
	SOI [CAD]	-30 -20	* X * X	* X *	* X * X	* X *	32.5 55	650 755
	CH4 [Mol %]	0	7.7	9.5	11.2
	O2 [Mol %]	20.9	19.3	19.0	18.6

Study matrix and the diagnostics applied, associated with CH4 and O2 concentration in the charge. * = High-speed Schlieren and shadowgraph imaging. X = High-speed OH* chemiluminescence imaging.

Acknowledgments

This work was supported by Institut CARNOT-ESP (Energie et Systèmes de Propulsion) which is gratefully acknowledge.

Contact Information

Pr. Frédéric


domain: qfin stat

text:
Yang Lu
email: luyang000278@gmail.com

Flexible (Panel) Regression Model for Bivariate Count-Continuous Data with Insurance Application

Keywords: mixed data, polynomial expansion, random effect, sequential forecasting/pricing

Flexible (panel) regression models for bivariate count-continuous data with an insurance application

Yang Lu

Introduction

This paper introduces a flexible regression model for mixed count-continuous panel data, with flexible correlated bivariate random effects (or unobserved heterogeneity). Such a model is useful for a wide range of applications, such as:

• number/cost of visits to physicians in health economics [see e.g. [START_REF] Duan | A comparison of alternative models for the demand for medical care[END_REF]],

• number/amount of shopping in marketing [see e.g. [START_REF] Paull | A generalized compound poisson model for consumer purchase panel data analysis[END_REF]],

• number/cost of trips in tourism forecasting [see e.g. [START_REF] Englin | Estimating social welfare using count data models: an application to long-run recreation demand under conditions of endogenous stratification and truncation[END_REF]].

Our paper will be focused on an area which necessitates sequential forecasting, that is the analysis of longitudinal insurance claims. Indeed, in many other applications, only cross-sectional data are used and the scientific interest is often limited to explaining the count and continuous variables by observable individual characteristics (or covariates). In Insurance, regulations require the yearly insurance premium to be updated regularly, in order to take into account each policyholder's own past history. 1 In this context, the bivariate count-continuous panel data arises since data are typically aggregated annually and are available in terms of annual claim count and the corresponding total claim cost. While the traditional pricing method has focused on the counts only [see e.g. [START_REF] Dionne | A generalization of automobile insurance rating models: the negative binomial distribution with a regression component[END_REF]], there are several benefits to consider both total cost and count, as well as their interdependence. First, it increases the forecast precision of the total cost, which is the ultimate variable of interest for an insurance company. Second, when only the count variable is taken into account, the pricing system penalizes unfairly customers with small previous claims compared to those who file large claims. This induces an efficiency loss for the insurance market [see [START_REF] Einav | Optimal mandates and the welfare cost of asymmetric information: Evidence from the uk annuity market[END_REF]], and gives individuals the incentive to (strategically) not hide minor accidents, if the expected reimbursement is less than the future premium increase2 . Therefore, of vital importance is a model that i) takes into account both observable covariates, as well as random effects; ii) captures, in a flexible way, the dependence between the count and continuous (size) variable; iii) allows for tractable Bayesian forecasting formulas. Our paper proposes a model unifying these properties.

More precisely, we introduce a bivariate regression model with correlated random effects.

They induce not only contemporaneous dependence between the two components, but also serial correlation, which serves as the basis for the forecasting and pricing of future claims. Our joint distribution of the heterogeneities is semi-parametric, on the contrary to the standard gamma heterogeneity assumption. This leads to flexible (marginal and joint) distributions for the observable bivariate count-continuous variable.

The proposed joint density function of the random effects takes the form of polynomial expansion with respect to a reference density and can approximate any distribution function arbitrarily well. It has a similar spirit as the model of [START_REF] Gurmu | Semiparametric estimation of count regression models[END_REF]; [START_REF] Gurmu | Flexible bivariate count data regression models[END_REF]] for cross-sectional count data. Our contribution to this latter literature is twofold. First, we propose a simple, and formal justification of this method. This also leads us to a comparison with another stream of literature, which employs a different polynomial expansion method to approximate the conditional transition density of the asset price dynamics. In particular, it is shown that our approach requires a much weaker assumption, and ensures positivity of the density function.

Second, we propose an equivalent parameterization of the polynomial expansion-based density function. This allows us to derive significantly simpler expressions for the likelihood function and, more importantly, the forecasting formulas. This should further increase the appeal of this semi-parametric class of models in empirical studies.

The rest of the paper is organized as follows. Section 2 presents the model. Section 3 provides an empirical illustration using a database of vehicle insurance claims. Section 4 concludes. Proofs are gathered in Appendix.

The model

The general setting

For each individual i, we denote by N i,t , t = 1, ..., T a sequence of count variables and Y i,t taking values in [0, ∞[, as well as a set of covariates X i . For expository purpose we assume the latter to be time-invariant and the extension to the time-varying case is straightforward. In insurance applications, the count variable corresponds to the number of claims during period t, and the continuous variable corresponds to the total cost of these N i,t claims, respectively. As a convention, Y i,t is 0 if N i,t equals zero. We assume that conditional on (static) unobserved heterogeneities U i,1 , U i,2 and observable covariates X i , the count variable N i,t is Poisson distributed with parameter λ i U 1,i , and Y i,t is Gamma distributed, with shape parameter δN i,t , and rate parameter

β i U i,2
, where:

λ i = exp(d 1 X i ), β i = exp(d 2 X i ). (1)
From now on, for expository purpose, we will assume that (X i , (N i,t ) t , (Y i,t ) t ) are i.i.d. across individuals and will thus omit the index i.

Thus we have:

E[Y t |U 1 , U 2 , N t ] = δNt βU2
, which is proportional to N t . Moreover, variable Y t has the following (conditional) compound representation:

Y t = Nt j=1 Z j,t , (2)
where Z j,t , j = 1, .. Moreover, the dependence between U 1 and 1/U 2 can be either positive or negative. In the first case, individuals with a higher claim frequency tend to have more severe claims. This can be explained by the fact that both the count and total cost reflect, to some extent, the skill of the driver. Alternatively, this association can also be negative. One potential explanation is the so-called "hunger for bonus" phenomenon. Indeed, some policyholders strategically choose not to report claims whose cost is smaller than the cumulated future premium increase [see e.g. [START_REF] Lemaire | Bonus-malus systems in automobile insurance[END_REF]]. This would lead to a higher (resp. lower) claim frequency, but a lower (resp. higher) cost per claim for honest individuals (resp. arbitrageurs who do not report small claims). The evidence of such negative dependence has recently been found by [START_REF] Garrido | Generalized linear models for dependent frequency and severity of insurance claims[END_REF].

As a comparison, the existing literature on mixed discrete-continuous regression usually considers only either of these two types of dependence. For instance, [START_REF] Fitzmaurice | Regression models for a bivariate discrete and continuous outcome with clustering[END_REF]; [START_REF] Yang | Regression models for mixed poisson and continuous longitudinal data[END_REF]; de [START_REF] De Leon | Copula-based regression models for a bivariate mixed discrete and continuous outcome[END_REF] specify the joint density between a discrete and a continuous variable via the marginal distribution of the discrete variable and the conditional distribution of the continuous variable given the discrete one. This corresponds to the first type of dependence mentioned above. Due to the lack of random effects, these models are not suitable when extended to a panel context, as serial correlation cannot be introduced between observations of different periods.

On the other hand, [START_REF] Catalano | Bivariate latent variable models for clustered discrete and continuous outcomes[END_REF]; [START_REF] Sammel | Latent variable models for mixed discrete and continuous outcomes[END_REF]; [START_REF] Gueorguieva | A correlated probit model for joint modeling of clustered binary and continuous responses[END_REF] introduce dependence via random effects only. However, these papers either concern binary discrete variable and not count variable, or involve rather non tractable, simulation-based algorithm for estimation and forecasting [see e.g. [START_REF] Sammel | Latent variable models for mixed discrete and continuous outcomes[END_REF]; [START_REF] Dunson | Bayesian latent variable models for clustered mixed outcomes[END_REF]].

Let us now propose a specification for the joint density of (U 1 , U 2 ), before studying its implications for the estimation and forecasting.

The model for the random effects

We assume that the joint density function of (U 1 , U 2 ) is:

g(u 1 , u 2 ) = 1 M e -c1u1-c2u2 u α1-1 1 u α2-1 2 J j=0 J k=0 b j,k u j 1 u k 2 2 , ( 3
)
where:

• coefficients (b j,k ) are real numbers, and b 0,0 = 1 for identification purpose.

• the normalization constant is equal to:

M = J j1,k1,j2,k2=0 b j1,k1 b j2,k2 Γ(j 1 + j 2 + α 1 )Γ(k 1 + k 2 + α 2 ) c j1+j2+α1 1 c k1+k2+α2 2 = 2J j=0 2J k=0 Γ(j + α 1 )Γ(k + α 2 ) c j+α1 1 c k+α2 2 Π j,k = W 1 (0) ΠW 2 (0), (4)
where matrix (Π j,k ) 0≤j,k≤2J is defined by:

Π j,k = j1+j2=j 0≤j1,j2≤J k1+k2=k 0≤k1,k2≤J b j1,k1 b j2,k2 , (5)
and vector functions

W 1 = (W 1,j ) 0≤j≤2J , W 2 = (W 2,k
) 0≤k≤2J are given by:

W 1,j (s 1 ) = Γ(j + α 1 ) (c 1 + s 1 ) j+α1 , W 2,k (s 2 ) = Γ(k + α 2 ) (c 2 + s 2 ) k+α2 , ∀0 ≤ j, k ≤ 2J, s 1 , s 2 ≥ 0. (6)
Thus g(u 1 , u 2 ) is a linear combination4 of gamma product densities:

f J (u 1 , u 2 ) = 1 M 2J j=0 2J k=0 e -c1u1 u α1+j-1 1 e -c2u2 u α2+k-1 2 Π j,k = e -c1u1 u α1-1 1 e -c2u2 u α2-1 2 X(u 1 ) ΠX(u 2 ) W 1 (0) ΠW 2 (0) ,
where vector functions X(u) := (1, u, u 2 , ..., u 2J ) for each u.

The resulting model is semi-parametric in the sense that the distribution of the unobserved heterogeneity is flexible and can approximate any distributions, whereas the regressors λ and β involve a finite dimensional parameter d. Similar polynomial expansions have been previously

proposed by [START_REF] Gallant | Semi-Nonparametric Maximum Likelihood Estimation[END_REF]; [START_REF] Gurmu | Semiparametric estimation of count regression models[END_REF]; [START_REF] Bierens | Semi-Nonparametric Interval-Censored Mixed Proportional Hazard Models: Identification and Consistency Results[END_REF]. Its background is the orthonormal projection of √ g 0 in an appropriate L 2 space, using a polynomial basis, whereas the square ensures the positivity of the density [see e.g. [START_REF] Gallant | Semi-Nonparametric Maximum Likelihood Estimation[END_REF]]. However, unlike the conventional approach which explicitly introduces these orthogonal polynomials, in this paper we use, without loss of generality, the canonical polynomials. This is due to the fact that both the orthogonal polynomials and the canonical polynomials are basis of the set of all polynomials, but the latter is much easier to deal with when used to the square (since the product of two canonical polynomials is still a canonical polynomials). This choice has the advantage of greatly simplify the computation of the likelihood function in the full model.

From the representation of combination of gamma densities, we derive the marginal distribu-tions of U 1 and U 2 :

g 1 (u 1 ) = e -c1u1 u α1-1 1 X(u 1 ) ΠW 2 (0) W 1 (0) ΠW 2 (0) , ( 7
)
g 2 (u 2 ) = e -c2u2 u α2-1 2 W 1 (0) ΠX(u 2 ) W 1 (0) ΠW 2 (0) , ( 8
)
which are also linear combinations of gamma densities. Thus we have the following property:

Proposition 1. The two components U 1 and U 2 are independent if and only if all the coefficients b jk are separable, that is, if for all j and k, we have:

b jk = b j,0 b 0,k .
In other words, our model can accommodate for the special case where the two random effects are independent.

Proof. See Appendix.

Similarly, the joint Laplace transform has also a closed form expression:

L(s 1 , s 2 ) = ∞ 0 e -u1s1-u2s2 g(u 1 , u 2 )du 1 du 2 = W 1 (s 1 ) ΠW 2 (s 2 ) W 1 (0) ΠW 2 (0) , ∀s 1 , s 2 ≥ 0.
The density of (N t , Y t ). Conditional on U 1 , U 2 , the joint distribution of (N t , Y t ) has two components:

• The first is the point mass at (0, 0), with probability:

P[N t = Y t = 0 | U 1 , U 2 ] = e -λU1 .
• One continuous component with respect to the measure µ 1 ⊗ µ 2 , where µ 1 is the Lebesgue measure on positive integers, and µ 2 is the Lebesgue measure on positive real numbers.

This component has a density:

f (n, y | U 1 , U 2 ) = λ n U n 1 e -λU1 n! β δn U δn 2 y δn-1 e -βU2y Γ(δn) , ∀n ≥ 1, y > 0. ( 9
)
By integrating out the distribution of (U 1 , U 2 ), we get the unconditional distribution of (N 1 , Y 1 ).

The elementary probability of the degenerate component is:

P[N t = Y t = 0] = E[e -λU1 ] = 1 M 2J j=0 2J k=0 Γ(j + α 1 )Γ(k + α 2 ) (c 1 + λ t ) j+α1 c k+α2 2 Π j,k = V 1 (0, λ)ΠV 2 (0, 0) W 1 (0) ΠW 2 (0) , ( 10
)
where vector functions V 1 , V 2 are defined by:

V 1,j (n, s 1 ) = Γ(j + α 1 + n) (c 1 + s 1 ) j+α1+n , V 2,k (n, s 2 ) = Γ(k + α 2 + δn) (c 2 + s 2 ) k+α2+δn , ∀n, s 1 , s 2 , j, k.
Similarly, the density of the continuous component is:

f (n, y) = λ n β δn y δn-1 n!Γ(δn) E[e -λU1-βU2y U n 1 U δn 2 ] = λ n β δn y δn-1 n!Γ(δn) 1 M 2J j=0 2J k=0 Γ(j + α 1 + n) (c 1 + λ 1 ) j+α1+n Γ(k + α 2 + δn) (c 2 + β 1 y) k+α2+δn Π j,k , = λ n β δn y δn-1 n!Γ(δn) V 1 (n, λ) ΠV 2 (n, βy) W 1 (0) ΠW 2 (0) ∀y > 0, n ∈ N {0}. ( 11
)
Accounting for the heavy-tail of the loss distribution. One of the desirable properties of loss distributions in Insurance is their heavy-tail. Let us first prove that our model allows for heavy-tailed loss distribution:

Proposition 2. The distribution of the total cost Y t has a heavy tail. In particular, it has a finite mean if and only if:

α 2 > 1. (12)
Proof. See Appendix.

From now on we assume this condition to be satisfied, since otherwise the risk is not insurable.

Estimation and forecasting algorithm

The log-likelihood function. Let us assume that for each individual i, we observe (N i,t ) 1≤t≤T , (Y i,t ) 1≤t≤T and covariates X i . These pairs (N i,t , Y i,t ) are independent conditional on (U 1 , U 2 ).

Then the log-likelihood function of the model is:

(θ) = I i=1 log P N i,t = n i,1 , Y i,t = y i,t , t = 1...T | X i ∝ I i=1 log E e -Λ 1iT U1 U N iT 1 U δN iT 2 e -Λ 2iT U2 + t,Ni,t =0 δN i,t log y i,t β i,t -log Γ(δN i,t ) ∝ I i=1 log V 1 (N iT , Λ 1iT ) ΠV 2 (N iT , Λ 2iT ) W 1 (0) ΠW 2 (0) + t,Ni,t =0 log λ Ni,t i,t β δNi,t y δNi,t Γ(δN i,t ) ,
where N iT is the sum of counts, Λ 1T := T t=1 λ = T λ is the cumulative intensity, and Λ 2T := T t=1 βy t is the normalized cumulative cost, and the set of parameters is equal to

θ = (α 1 , α 2 , (b j,k ) j,k , δ, c 1 , c 2 , d 1 , d 2 ).
Forecasting formulas. Conditional on past claim history (N t ), (Y t ), t = 1, ..., T , the distribution of (U 1 , U 2 ) is:

f T (u 1 , u 2 |(N t ), (Y t ), t = 1, ..., T ) ∝ f (u 1 , u 2 ) T t=1 u Nt 1 e -λu1 T t=1 u δNt 2 e -βu2yt
Thus the expected total cost at time T + 1, given claim count and cost history is:

E[Y T +1 | (N t , Y t ), t = 1, ..., T ] = δ β t λ t E[U 1 /U 2 | (N t , Y t ), t = 1, ..., T ] = δλ β E[U N T +1 1 U δN T -2 2 e -Λ 1T U1-Λ 2T U2 ] E[U N T 1 U δN T -1 2 e -Λ 1T U1-Λ 2T U2 ] = δλ β 2J j,k=0 Π j,k Γ(j+α1+N T +1) (c1+Λ 1T ) j+α 1 +N T +1 Γ(k+α2+δN T -1) (c2+Λ 2T ) k+α 2 +δN T -1 2J j,k=0 Γ(j+α1+N 1T ) (c1+Λ 1T ) j+α 1 +N 1T Γ(k+α2+δN T ) (c2+Λ 2T ) k+α 2 +δN T Π j,k (13)
= δλ β V 1 (N T + 1, Λ 1T ) ΠV 2 (N T -1/δ, Λ 2T ) V 1 (N T , Λ 1T ) ΠV 2 (N T , Λ 2T ) .
We can also rewrite this expectation as a weighted average, indeed from (13) we have:

E[Y T +1 | (N t , Y t ), t = 1, ..., T ] = δλ β 2J j,k=0 (j + α 1 + N T )(c 2 + Λ 2T ) (c 1 + Λ 1T )(k + α 2 + δN T -1) ω j,k (14)
with:

ω j,k = Π j,k Γ(j+α1+N T ) (c1+Λ 1T ) j+α 1 +N T Γ(k+α2+δN T ) (c2+Λ 2T ) k+α 2 +δN T 2J j ,k =0 Γ(j +α1+N T ) (c1+Λ 1T ) j +α 1 +N T Γ(k +α2+δN T ) (c2+Λ 2T ) k +α 2 +δN T Π j ,k
.

In terms of insurance pricing, the term 14) is the product of two terms:

(j+α1+N T )(c2+Λ 2T ) (c1+Λ 1T )(k+α2+δN T -1) in (
• The first one j+α1+N T c1+Λ 1T is the posterior expectation of U 1 conditionally on previous counts only, if U 1 was gamma distributed;

• The numerator of the second term

c2+Λ 2T
k+α2+δN T -1 is, up to an additive constant, linear combination of previous claim costs whereas the denominator is, up to an additive constant, the cumulative claim counts. Thus this second term measures, to some extent, the average cost per claim during the first T periods.

Note also that due to the dependence between U 1 and U 2 , the posterior expected cost is generically not the product of expected future claim count and expected future average cost per claim 5 , that is:

E[N T +1 | (N t , Y t ), t = 1, ..., T ]E[Y T +1 /N T +1 | (N t , Y t ), t = 1, ..., T ] = E[Y T +1 | (N t , Y t ), t = 1, ..., T ].

Numerical application

In this section we fit our model to an insurance database, and analyze the potential bias in terms of pricing, when an independence assumption between U 1 and U 2 is imposed.

The data

Our data concerns an Australian vehicle insurance portfolio. It contains only one period of observation, hence in the following, the time index t will be omitted. It has initially been analysed by De Jong and Heller (2008) using generalized linear models and is available free of charge at the following website: www.businessandeconomics.mq.edu.au/our_departments/Applied_Finance_and_Actuarial_ Studies/research/books/GLMsforInsuranceData/data_sets

The database consists of 65324 (vehicle damage) policies, of which 4423 (6.8 percent) had at least one claim. For each policyholder, we observe the total claim amount (continuous valued, in 10000 $) and the total claim count (integer valued), as well as the following covariates 6 :

• type of the vehicle (categorical, with 6 different values 7 )

• vehicle's age (categorical, with 4 different values)

• gender (male or female)

• driver's area of residence (categorical, with 6 different values)

• driver's age (categorical, with 4 different values)

The following table summarizes the distribution of the count variable.

5 Except when U 1 , U 2 are independent. In this case, coefficients b j,k = b j b k are separable and as a consequence coefficients Π j,k , ω j,k are also separable.

6 Although the value of the vehicle is also available, it has not been taken into account. This is in line with the existing literature [see e.g. [START_REF] Pinquet | Allowance for the age of claims in bonus-malus systems[END_REF]; [START_REF] Czado | A mixed copula model for insurance claims and claim sizes[END_REF][START_REF] De | Generalized linear models for insurance data[END_REF]] and is motivated by several reasons: i), the value of the vehicle is partially correlated with other covariates, in particular the type of the vehicle; ii) for some policies, the recorded vehicle value is recorded as zero, which suggests a lack of credibility of this variable. One simple alternative would be to divide the amount of the claim by the value of the vehicle to get the "normalized claim cost". However this alternative has not been chosen due to certain policies with a zero vehicle value.

7 There are initially 13 different types of vehicles. Seven of which (bus, convertible, motorized caravan, minibus, coupe, panel van, and roadster) account for, in total, only 2200 policies, or 3.5 % of the sample. These policies are dropped for illustration purpose and we are left with five major types of vehicles (utility cars, sedan, trucks, station wagons, hardtop, as well as hatchbacks). Percentage 90% 95 % 97% 98 % 99% 100% Quantile 0.49 0.80 1.06 1.28 1.70 5.59

Number of claims

Table 2: Some quantiles of the cost variable, in 10000 dollars.

We can see, for instance, that the 99% quantile of the cost variable is approximately two (resp.

• Model M1, with J = 1, that is

f (u 1 , u 2 ) = 1 M e -c1u1-c2u2 u α1-1 1 u α2-1 2 (1 + b 01 u 1 + b 10 u 2 + b 11 u 1 u 2 ) 2 ,
• Model M2, with J = 2, that is:

f (u 1 , u 2 ) = 1 M e -c1u1-c2u2 u α1-1 1 u α2-1 2 (1+b 01 u 1 +b 02 u 2 +b 10 u 2 +b 11 u 1 u 2 +b 12 u 1 u 2 2 +b 20 u 2 1 +b 21 u 2 1 u 2 +b 22 u 2 1 u 2 2 ) 2 .
• In order to compare with the unconstrained Model M2, and quantify the induced potential actuarial pricing error9 , we also estimate Model M2bis, that is the constrained version of M2, under the independence condition b j,k = b j,0 b 0,k [see Proposition 1]. More precisely we have:

f (u 1 , u 2 ) = 1 M e -c1u1-c2u2 u α1-1 1 u α2-1 2 (1 + b 01 u 1 + b 20 u 2 1 ) 2 (1 + b 10 u 2 + b 02 u 2 2 ) 2 .
Due to the large number of parameters, the numerical optimization of the log-likelihood function is highly complicated. Thus it is essential to find an initial value for the optimization that is close to the real value of the parameter. Such an initial value can be found via two seemingly unrelated10 non-linear least square estimation, also called pseudo maximum likelihood

[see e.g. [START_REF] Gourieroux | Pseudo maximum likelihood methods: Applications to poisson models[END_REF]]. Indeed, we can remark that:

E[N i | X i ] = exp(d 1 X i )E[U 1 ] (15) E[Y i | X i ] = E E[Y i | X i , N i ] | X i = exp (d 1 -d 2 )X i E[U 1 /U 2 ]. ( 16
)
Thus d 1 and d 1 -d 2 , as well as

E[U 1 ], E[U 1 /U 2 ]
can be estimated by minimizing:

I i=1 (N i -exp(d 1 X i )E[U 1 ]) 2 ,
and

I i=1 Y i -exp (d 1 -d 2 )X i E[U 1 /U 2 ] 2 .
Once the initial values of the regression parameters are obtained, we estimate the different models by maximum likelihood. As expected, we can check that the estimators of the regression coefficients are quite similar across different models. Somehow surprisingly, the sign of the regression coefficients before the same covariate might be different for the count and the cost variable. In other words, the covariates might have opposite effects on the two response variables; this confirms also that it is important to accommodate for flexible dependence between the random effects U 1 , U 2 .

The following table reports the estimates of the remaining parameters, that are δ, as well as parameters characterizing the joint distribution of the unobserved heterogeneities. Table 5: Estimates of the remaining parameters for the three models. The symbol "-" indicates that a parameter is not involved in a particular model. Table 6: Information criteria of the four models.

We can observe that Model 2 is the best model in terms of BIC, although the difference between M2 and M1 is rather small. On the other hand, we can see that the benchmark model M0, as well as the constrained model M2bis have a BIC that is much larger than their more flexible competitors. To further illustrate this, let us analyze the heterogeneity distribution for the estimated model M2. The following figures plot the p.d.f. of the marginal distributions of U 1 and that of U 2 . We can see that Model M2 indicates a significant difference of the marginal distributions of both U 1 and U 2 from the gamma distribution assumption. Let us now plot the iso-densities of the joint distribution of (U 1 , U 2 ). Figure 4: Iso-densities of the joint distribution of (U 1 , U 2 ) obtained from Model 2. The lighter a region is, the larger the joint p.d.f. is in that region.

We can observe that the joint p.d.f. of (U 1 , U 2 ) have two distinct local maxima. This is to be compared with the joint p.d.f. in Model M0, in which U 1 , U 2 are independent gamma (in this case, there is only one global maximum). This implies that the joint p.d.f. of (U 1 , U 1 /U 2 ) has also two modes, one mode corresponds to individuals with both larger values of U 1 and 1/U 2 , or in other words, the riskier group; the other one corresponds to individuals with smaller values of U 1 and 1/U 2 , that are the safer group. This positive association between U 1 and 1/U 2 is also reflected by the correlation coefficient between these two terms, which is equal to 0.33.

Impact of wrong dependence assumptions on actuarial pricing

Let us now study the impact on the actuarial premium of the incorrect assumption of independence between U 1 and U 2 . Given the claim counts N 1 , N 2 , ..., N T , as well as total costs Y 1 , Y 2 , ..., Y T , we compute next year's actuarial fair premium, that is the conditional expected future cost Y T +1 of the claim, under Model M2, as well as the nested model11 M2b, which assumes independence between U 1 and U 2 .

For illustration purpose, we consider a policyholder for which the values of the covariates are all equal to 0, that is, X 1 = X 2 = .... = X 17 = 0. In this case, we have λ t = β t = 1 for all t, and Λ 1T = T is equal to the number of periods, whereas Λ 2T = Y T is equal to the cumulative claim cost. As a consequence, in this special case, the premium is a function of cumulative counts

N T = N 1 + • • • + N T , and cumulative costs Y T = Y 1 + • • • + Y T ,
N T = 0, Y T = 0, T = 0 147 $ 147 $ N T = 0, Y T = 0, T = 1 144 $ 142 $ N T = 0, Y T = 0, T = 3 139 $ 134 $ N T = 1, Y T = 300 $, T = 1 159 $ 163 $ N T = 1, Y T = 300 $, T = 3 153 $ 154 $ N T = 1, Y T = 2000 $, T = 1 263 $ 419 $ N T = 1, Y T = 2000 $, T = 3 253 $ 396 $ N T = 1, Y T = 5000 $, T = 1 454 $ 904 $ N T = 1, Y T = 5000 $, T = 3 436 $ 853 $
Table 7: Examples of posterior premium at period T + 1, given claim counts and costs of the first T periods.

We observe that the two models provide the same a priori premium (first row of the previous table), that is when no previous claim history is available. This is expected, since the a priori premium is equal to E[U 1 /U 2 ], which can be estimated consistently by the same quasi likelihood method. Once past claim history begins to accumulate, the premium given by the two models start to differ quite significantly. From the second and third row, we can see that Model 2 provides more premium reduction than Model 2, when 0 past claim has been reported during the past observation period. This is expected, since in Model 2, the unobserved heterogeneities U 1 and 1/U 2 are positively correlated. Thus the absence of claims indicates that the individual has, in average, a smaller value of U 1 , and, as a consequence, a smaller value of 1/U 2 . From the fourth and fifth row, when one accident of value 300 $ is reported (which is significantly smaller than the average claim cost, for those who report claims), Model 2 charges a higher premium increase than Model M2bis. This can be explained by the following. The presence of a claim indicates that the individual has a larger value of U 1 , but the small size of the claim indicates a smaller value of 1/U 2 . In Model M2bis, where dependence is not taken into account, these two effects are well compensated. This is less the case in Model 2, where due to the positive dependence between U 1 and 1/U 2 , the two effects are only partially compensated, hence the higher premium in Model 2. For the same reason, in the sixth and seventh row, where a claim of average value (2000 $) is reported, or in the eighth and ninth row, where a claim of large value (5000 $) is reported, Model M2 charges a higher premium than Model M2bis. Moreover, the larger the claim size, the more important the difference between the premium given by the two models. To summarize, compared to Model M2, Model M2bis overcharges low-risk customers, but undercharges high-risk ones. This constitutes a huge disadvantage on a competitive market.

Indeed, an insurer using Model M2bis, say A, is likely to lose low-risk customers (who prefer its competitor, say B, operating under Model M2 and hence provide lower premium for these customers). On the other hand, Insurer A is more likely to take in high-risk customers (who leave insurer B to seek lower premium), who are undercharged by insurer A compared to their real risk. As a consequence, insurer A will face solvency issues in the long-run, since the premium income is likely to be insufficient to cover the realized claims. Therefore, Model M2bis provides not only a less satisfactory fit12 , but also significant pricing bias.

Conclusion

We have proposed a general, flexible random effect panel data model for bivariate count-continuous variables, that is especially adapted for insurance forecast and pricing. The model allows for convenient, matrix-based expressions for the likelihood function, as well as the forecast formulas.

Finally, we have demonstrated the benefit such a model offers for the forecasting of future insurance claims, compared to existing parametric specifications.

The simulated data and summary statistics of the estimates

The Data Generating Process (DGP) we propose is the following: in the joint density formula of (U 1 , U 2 ), we take:

α 1 = α 2 = 2, c 1 = c 2 = 1, J = 1, b 00 = 1, b 01 = b 10 = 0, b 11 = 0.5,
where b 00 is set to one by convention. Note also that for illustration purpose, we have assigned nonnegative values to all the parameters b j,k . The advantage of this set of parameters is that the resulting joint density is a mixture distribution and is therefore easy to simulate. 14

Then we set the two regressors in (1) as:

λ = exp(β 1 X 1 ), β = exp(β 2 X 2 ),
where β 1 = 0.5, β 2 = 0.6, and X 1 , X 2 are i.i.d. N (0, 1) distributed. Finally, parameter δ is set to 1. of the estimates along with their standard deviation.

14 Nevertheless, even when some of the coefficients are negative, it is still possible to use a acceptance-rejection algorithm to simulate samples of (U 1 , U 2 ). See [START_REF] Gallant | A nonparametric approach to nonlinear time series analysis: estimation and simulation[END_REF] for details.

All of the parameter estimates are significant, except those of b 01 and b 10 . This is expected since their true values are exactly zero.

Impact on pricing

Next, let us examine how the sampling error might affect the forecasting performance. To this end, for each of the 50 simulated databases, we use the estimated model to compute the conditional expectation E[Y 2 |Y 1 , N 1 , X] using (13). We then compute the sample mean and variance of E[Y 2 |Y 1 , N 1 , X]. The following figure plots the histogram of these estimated forecasts, for X 1 = X 2 = 0, N = 1 and Y = 1.6. Of these 50 samples of E[Y 2 |Y 1 , N 1 , X], the largest error is around 0.3, which is about 10 % of its theoretical value.

It is useful to compare this simulation exercise with the real data application of Section 3.

There are several differences between the previous simulation and the insurance data. Firstly, the insurance data contains more than 60000 observations, that is more than 6 times larger than the simulated datasets, which should contribute to lower the sampling error of the maximum likelihood estimator. Secondly, as is often the case in practice, the real data we analyzed contains much more covariates than the simulation experiment. 15 One can reasonably argue that this spells a larger number of parameters and thus can adversely affect the precision of the maximum likelihood estimator. Nevertheless, we argue that in practice this impact can be mitigated by the

Figure 1 :

 1 Figure 1: Histogram of total claims.

  59 = b 01 b 10 = -0

Finally

  , let us provide a comparison of the four models in terms of information criteria.

Figure 2 :Figure 3 :

 23 Figure 2: Marginal p.d.f of U 1 obtained from Model M2.

Figure 5 :

 5 Figure 5: Histogram of the estimated forecasts.

Table 1 :

 1 Number of policies with the given number of claims Number of policies according to the number of claimsThe average number of claims per individual is approximately 0.073, whereas the variance of the number of claims is equal to 0.077. Thus the count variable is slightly over-dispersed.Figure1plots the histogram of Y |N > 0, that is the distribution of total claim cost, for those who reported claims. We can remark that this distribution has both a heavy left tail, corre-

	0	60901
	1	4149
	2	255
	3	17
	4	2

sponding to a large proportion of minor-sized claims, as well as a heavy right tail, corresponding to extreme value claims.

Table 4 :

 4 The following table reports first the estimates of the regression coefficients. Estimates of the regression coefficients.

			Regression coefficient			Regression coefficient
	Variable		for the count variable			for the cost variable
		M0	M1	M2	M2bis	M0	M1	M2	M2bis
	X 1	-0.0201 -0.0235 -0.0226 -0.0217 0.166	0.178	0.152	0.176
	X 2	-0.158	-0.146	-0.142	-0.153	0.436	0.468	0.421	0.456
	X 3	-0.114	-0.132	-0.125 -0.0120 0.426	0.453	0.409	0.431
	X 4	-0.131	-0.146	-0.126	-0.128	0.340	0.369	0.324	0.352
	X 5	-0.287	-0.264	-0.259	-0.271	0.416	0.436	0.462	0.427
	X 6	-0.185	-0.169	-0.199	-0.175	0.260	0.245	0.236	0.254
	X 7	0.418	0.403	0.426	0.423	-0.361 -0.335 -0.374 -0.365
	X 8	0.254	0.229	0.260	0.249	-0.120 -0.135 -0.114 -0.126
	X 9	0.228	0.196	0.264	0.221	-0.052 -0.046 -0.038 -0.044
	X 10	0.198	0.185	0.201	0.193	-0.031 -0.029 -0.036 -0.033
	X 11	0.0156	0.0192	0.0144	0.0164	0.107	0.095	0.114	0.112
	X 12	0.0872	0.0893	0.0869	0.086	0.107	0.113	0.093	0.112
	X 13	0.225	0.218	0.197	0.230	0.0896 0.0878 0.0852 0.905
	X 14	0.0838	0.0866	0.0842	0.0893	0.606	0.572	0.632	0.664
	X 15	0.146	0.154	0.142	0.153	0.106	0.098	0.120	0.114
	X 16	0.172	0.169	0.182	0.184	0.254	0.268	0.224	0.263
	X 17	0.358	0.326	0.335	0.349	-0.026 -0.018 -0.014 -0.029
	X 18	0.214	0.228	0.241	0.223	-0.127 -0.132 -0.142 -0.128
	X 19	0.243	0.226	0.238	0.239	0.122	0.136	0.114	0.126

  and T only. The following table provides some examples of the values of the premium for different values of N T , Y T , and T .

	Model M2b Model 2

Table 8 :

 8 The following table reports the theoretical mean and variances of the two variables N, Y . Summary statistics of the couple (N, Y ).We can see that for this DGP, the count variable is quite heavily overdispersed. As pointed out by a referee, this situation is commonly encountered in practice for low count data.The following table reports the sample mean and standard deviation of the estimates, obtained respectively from the 50 simulated databases.

	Variable Mean Variance Variance/Mean
	N	7.6	30.6	4.04
	Y	3.2	29.3	9.2
	Variable Theoretical value Empirical Mean Standard Deviation
	β 1	0.500	0.501	0.12
	β 2	0.600	0.598	0.11
	α 1	2.00	2.02	0.13
	α 2	2.00	2.10	0.14
	c 1	1	1.06	0.07
	c 2	1	0.97	0.08
	b 01	0	0.01	0.06
	b 10	0	0.01	0.05
	b 11	0.5	0.052	0.2
	δ	1	0.97	0.2

Table 9 :

 9 Table

Such pricing method is called bonus-malus.

This is called hunger for bonus phenomenon.

In particular, if U 2 follows gamma marginal distribution, then the distribution of each individual cost is called generalized beta of the second kind (GB2) [see e.g.[START_REF] Mcdonald | Some generalized functions for the size distribution of income[END_REF];[START_REF] Cummins | Applications of the gb2 family of distributions in modeling insurance loss processes[END_REF];[START_REF] Gouriéroux | The econometrics of risk classification in insurance[END_REF];[START_REF] Frangos | Design of optimal bonus-malus systems with a frequency and a severity component on an individual basis in automobile insurance[END_REF]].

But with possibly negative weights, since the entries of matrix B need not to be necessarily nonnegative.

For instance, if one categorical variable can take 4 different values, then it is transformed into three linearly independent binary variables.

The current actuarial literature usually assumes independence between U 1 and U 2 , see e.g.[START_REF] Tzougas | Optimal bonus-malus systems using finite mixture models[END_REF].

That is to say, by neglecting the dependence between N and Y . See e.g.[START_REF] Zellner | An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias[END_REF] for a discussion.

The choice of the alternative model M2bis is motivated by the fact that M2bis is flexible when it comes to modeling the marginal distributions of U 1 and U 2 . Thus Models M2 and M2bis differ only in their ability of capturing the dependence. One could also compare, say, Models M2 and M0, but in this case the pricing errors would be more difficult to analyze, since both the marginal distributions and the dependence structure contribute differ between the two models.

In terms of BIC, for instance.

This relatively small number of replications is explained by the fact that simulating each database, as well as computing the maximum likelihood estimator, is quite computationally intensive.

In our simulation experiment, we have deliberately chosen only two covariates, in order not to mix the error induced by adding too many regression coefficients, as well as those induced by the novel specification of the random effects (U 1 , U 2 ). Such a choice is rather standard in this literature and has been previously made in[START_REF] Gurmu | Semiparametric estimation of count regression models[END_REF];[START_REF] Gurmu | Flexible bivariate count data regression models[END_REF] for the analysis of count data.

four) times larger than the 95 % (resp. 90 %) quantile, whereas the biggest claim is another three times larger. The heavy tail of the distribution of cost implies that standard models based on gamma-distributed claim costs [see e.g. [START_REF] Czado | A mixed copula model for insurance claims and claim sizes[END_REF]] may not be appropriate for this data.

The following table summarizes the available covariates. All covariates are discrete and transformed into binary variables 8 . For identification issue, the constant is not included as a covariate, and the means

1 if the area of residence is "C" X 5

1 if the area of residence is "D" X 6

1 if the area of residence is "E" X 7

1 if the driver's age is "1" (youngest) X 8

1 if the driver's age is "2" X 9

1 if the driver's age is "3" X 10 1 if the driver's age is "4" X 11

1 if the driver's age is "5"

1 if the vehicle is a station wagon

Estimation

In this subsection we estimate the model for three different values of J:

, which corresponds the density of independent gamma variables;

Appendix 1: proofs

Proof of Proposition 1.

We have independence if and only if g(x, y) = g 1 (x)g 2 (y), which is equivalent to:

where P 1 = E1(u1) ΠW2( 0) W1( 0) ΠW2(0) , P 2 = W1(0) ΠE( u2) W1( 0) ΠW2(0) are polynomials. By taking u 2 = 0, we have,

Thus we have finally:

which yields b jk = b j,0 b 0,k for integers j, k between 0 and J.

Proof of Proposition 2.

The marginal density of y has two components, one mass at zero, and one continuous component on R >0 . The density of this component is:

By equation ( 11), the dominant term of this function is proportional to y -1-α2 , when y goes to infinity. Thus E[Y t ] exists if and only if α 2 > 1.

Appendix 2: A simulation study

In order to get a sense of the finite sample behavior of the estimation method as well as its impact on pricing, in this section we propose a Monte-Carlo simulation exercise. For a fixed vector of parameters θ, we simulate 50 databases 13 , which one consisting of 10000 i.i.d. samples of (N 1 , Y 1 , X 1 ). Then we conduct maximum likelihood estimation on these 100 simulated databases.

fact that these regression parameters can be quite easily estimated separately by nonlinear least square [see equations ( 15) and ( 16)]. Also, Lasso-like methods can also be used to diminish the number of regression parameters, although it is clearly out of the scope of the present paper.