Thibault Désert

Thierry Authié

Sébastien Trilles

Modelling of the ionosphere by neural network for equatorial SBAS

The estimation of the ionosphere delay and associated confidence interval constitutes the major issue to reach APV1 availability performance level for single frequency SBAS above the equatorial area.

The ionosphere is a complex physical system which dynamics is particularly disturbed at the Geomagnetic Equator while mid-latitude regions are quieter. Classical methods to compute ionosphere delays, such as those implemented in EGNOS and the WAAS, are specific to a smooth ionosphere behavior and are not really adapted to follow high spatial and temporal gradients, such as those observed in the equatorial area. Thus innovative methods, having flexible and reactivity qualities, shall be defined and adapted to propose efficient equatorial SBAS.

Classically in SBAS concept, the knowledge of the ionosphere delay is obtained by a set of lines-of-sight between the network of ground stations and the navigation satellite constellation. Each line of sight intersects the ionosphere layer, assumed infinitely thin, and the dual-frequency combination allows to compute, at first order, the ionosphere delay that affects the GNSS measurements.

From this set of heterogeneous information, locally sampled irregularly on the sphere and changing over time, we propose to build an interpolating method to calculate the ionosphere delay on a point of interest by an adaptive mesh, unlike fixed grids usually used.

In this paper, we introduce an interpolation method based on the definition of a flexible network that can adapt to the spatial location of the data. It is therefore proposed to create self-organizing maps -as the Kohonen mapdefined by a mesh that fits the data. This network ends up sticking to the data by "learning" in real time: the mesh becomes denser and denser in the presence of many measurements and relaxes otherwise. This technique increases the granularity of the ionosphere delay information to compute, in particular to be able to describe the local plasma bubbles or depletions, if they are observable.

The adaptive array technology "learning" has been widely studied in the field of modeling neural networks. Their main advantage is to be able to reach an optimal state based on the information they process. The experimentation based on this technique shows a very good behavior in the case of strongly disturbed ionosphere conditions and the preliminary results are promising to bring the expected robustness to deploy SBAS in equatorial area.

INTRODUCTION

The upper layers of Earth's atmosphere, which constitute the ionosphere, perturb the propagation of GNSS signals [START_REF] Goodman | Space Weather & Telecommunications[END_REF], [START_REF] Arbesser-Rastburg | Effects on satellite navigation[END_REF]. These disturbances degrade the accuracy of the distance measurements performed by GNSS receivers, and thus the user location. Besides the diffraction effects, the ionosphere changes the propagation of GNSS signals. Indeed the carrier wave and its modulation (codes) travel at different effective velocities. The code propagation effective velocity decreases and becomes slightly lower than the speed of light, resulting in a delay in distance measurement. The ionosphere delay depends essentially on the signal frequency and density of free electrons in the ionosphere. To address this degradation, SBAS broadcasts to the aviation user ionospheric correction parameters and integrity data improving the position computation.

The SBAS Minimum Operational Performance Standard document (MOPS, [START_REF]Minimum operational performance standards for Global positioning system / wide area augmentation system airborne equipment[END_REF]) specifies the way the augmentation data should be transmitted to the aviation users. This document specifies that the ionosphere shall be considered as a thin layer, located at 350 km altitude in the WGS reference frame. Consequently all ionosphere effects are assumed to be concentrated and applied on a single point, the Ionospheric Pierce Point (IPP) defined as the intersection between the thin layer and the line of sight that connects the aviation user and the GNSS satellite. The MOPS also specifies a discretization of the ionosphere layer by a regular square grid (5 deg per 5 deg except around the pole areas) whose vertices are called Ionosphere Grid Point (IGP). The Message Type 26 of [START_REF]Minimum operational performance standards for Global positioning system / wide area augmentation system airborne equipment[END_REF] contains the Grid Ionospheric Vertical Delay (GIVD) and the Grid Ionospheric Vertical Error (GIVE, a 1.10 -7 confidence interval) on each IGP defined and identified by the Message Type 18 of [START_REF]Minimum operational performance standards for Global positioning system / wide area augmentation system airborne equipment[END_REF]. Finally the MOPS defines a mapping function, that converts the vertical ionosphere delay into a slant one, and an interpolation method that allows aviation users to compute the ionosphere delay along the line of sight at the IPP from the surrounding IGPs.

In the SBAS context, the ionosphere correction message (MT26) is the one that has the most impact on the system performance, mainly user availability, for single frequency service levels. Consequently the development of an equatorial SBAS faces several difficulties that come from the ionosphere modeling capability. The equatorial area faces strong ionosphere dynamics, such as large spatial and temporal gradients, scintillation effects [START_REF] Kintner | PAULA GPS and ionospheric scintillations[END_REF], local plasma bubble [START_REF] Kelley | Convective ionospheric storms: A review[END_REF], [START_REF] Huang | Occurrence probability and amplitude of equatorial ionospheric irregularities associated with plasma bubbles during low and moderate solar activities (2008-2012)[END_REF]. So the challenge is to define a suitable ionosphere model able to track such gradient events so as to fit the TEC field observed through the measurements and to minimize the GIVD error.

For equatorial SBAS, other issues can arise by applying the MOPS computation at user level. First the MOPS bilinear interpolation function can introduce an error in the computation of the User Ionospheric Vertical Delay (UIVD) at each IPP. This error increases in case of high non-linear ionosphere dynamics (common above the equatorial area) and is amplified due to the large scale of the IGP distribution. However the SBAS has to protect the user against this error such that the User Ionospheric Vertical Error (UIVE), resulting from the interpolation of four GIVE associated to the surrounding IGP, shall contain the error between the real vertical delay and the one rebuilt using bilinear interpolation. Second the mapping function can produce errors in the estimation of the User Ionospheric Slant Delay (UISD) along the line of sight mainly at low elevation. Nevertheless these questions are not discussed here, the purpose is focused on the minimization of the GIVD error by using adaptive map.

PROBLEM CHARACTERIZATION

The characteristic scale depends on the ionosphere conditions (geomagnetic storms or quiet periods). At low latitudes the ionosphere presents a global dynamics with, for the area of interest (i.e. the equatorial area), spatial gradients of the order of 30-50 TECU over 1000 km of distance [START_REF] Hocke | A review of atmospheric gravity waves and travelling ionospheric disturbances: 1982-1995[END_REF] (1 TECU can be seen as a 16 cm delay for the L1 signal). Thus, changes in gradients of less than 1000 km distances are common features of the ionosphere at these latitudes.

Another physical phenomenon that affects the ionosphere behavior above the equatorial areas are the plasma bubbles. A plasma bubble is a strong electron density contrast area, starting from the lower part of the ionosphere and rising in altitude (high level of ionosphere), and follows the magnetic field lines in the North West direction [START_REF] Kelley | Convective ionospheric storms: A review[END_REF]. The magnetic field strongly controls the generation and behavior of bubbles. In Africa, the magnetic field is smoother than on the South American sector (South Atlantic Anomaly), but the plasma bubbles are more common there, a phenomenon that is not yet fully explained [START_REF] Huang | Occurrence probability and amplitude of equatorial ionospheric irregularities associated with plasma bubbles during low and moderate solar activities (2008-2012)[END_REF].

Finally scintillation effects cause amplitude fading on code phase measurements (10dB observed above Kourou's station [START_REF] Secretan | SAGAIE: A Network of GNSS station to investigate the ionosphere behaviour in sub-Saharan regions[END_REF]) and phase jitter on the carrier phase measurements. The effect of scintillation on the carrier phase can generate cycle slips and potentially a loss of lock. Scintillations are small-scale plasma irregularities that compose the ionosphere, which translate into a rapid change in the refractive index for the radio signals [START_REF] Kintner | PAULA GPS and ionospheric scintillations[END_REF]. Plasma bubbles are at the origin of the scintillation to lower latitudes. Around equatorial area, scintillation is a nocturnal phenomenon that occurs in the first part of the night.

STATE OF THE ART ON ESTIMATION TECHNIQUES

Let us consider first the methods used in EGNOS and the WAAS for a short description.

The EGNOS system uses an approximation of the ionosphere thin layer, by a polyhedron which facets are triangles (Triangular Interpolation, TRIN). The technique has been developed by Mannucci and presented in [START_REF] Mannucci | A new method for monitoring the Earth's Ionospheric Total Electron Content Using the GPS Global Network[END_REF]. Each vertex of the TRIN mesh has a constant solar local time and computes the vertical ionosphere delay with a Kalman filter. After computing the vertical TEC (Total Electron Content) on the nodes of the TRIN, the GIVD are deduced by a second triangular interpolation (see [START_REF] Alleau | Comparative Ionosphere Electron Content Estimation Method in SBAS Performances[END_REF] for a detailed description). The advantages are first to keep a ionosphere history and second to avoid the Kalman filter following high dynamic night-and-day transition. The limitation of the TRIN is that the method assumes that the ionosphere evolves linearly over a facet. In [START_REF] Trilles | Adaptative ionospheric electron content estimation method[END_REF] some new ideas have been developed to propose an adaptive TRIN mesh. In this context the algorithm creates additional TRIN nodes around a detected local gradient and removes it when normal conditions are recovered.

The WAAS has implemented a Kriging estimation, that belongs to the class of geostatistical interpolators [START_REF] Blanch | An Ionosphere Estimation Algorithm for WAAS Based on Kriging[END_REF]. This interpolation provides the best unbiased estimatorin the sense of minimal variance -of the variable to interpolate using linear combinations of neighboring values. In this method, the variable is assumed to be composed of a deterministic part and a random part that is not a Gaussian white noise process, but contains a spatial dependency structure. The advantages are first to solve the GIVD in a single iteration and minimize the number of interpolation operations, and second to provide a method that reduces the GIVE without impacting integrity. One limitation of Kriging is that the method assumes a spatial isotropy of the data. To cope with this difficulty the algorithm may be more complex by defining two variograms, one according to the latitude and the other according to the longitude.

Other techniques exist and are widely used for modeling the ionosphere layer. The spherical harmonics decomposition is a well-known technique to model a thin ionosphere layer. The accuracy of the model depends on the order of the decomposition, and the resolution with high order level depends on the quantity of observables. This approach is efficient for middle latitude areas where the ionosphere behavior is quite smooth. The GIVD computation is given directly by the analytic value of the model at each IGP location.

A Taylor approximation consists in modeling locally the thin ionosphere layer by a polynomial (of degree two classically). The polynomial coefficients are fit using the vertical TEC at IPPs. This technique can be used to compute directly the GIVD on each IGP by adjusting a parabola above each IGP. In this context, the GIVD is identified as the zeroth order coefficient of the polynomial.

Finally, another technique can be developed coming from the image processing and signal reconstruction field [START_REF] Alleau | Comparative Ionosphere Electron Content Estimation Method in SBAS Performances[END_REF].

The Adaptive Normalized Convolution method can be viewed as a Local Weighted Mean Least Square; the weighting kernel is given by a local gradient estimation. First the method determines the gradient of samples inside a local ball of a Gaussian filter around the desired point to interpolate (IGP). Second the method deforms the base of the filter in the direction of the lowest gradient in order to interpolate among values evolving slowly. The method is designed to restitute the contour of a rapid change in the data. It is also a direct IGP interpolation from data measurements taking into account the anisotropy of the data.

In middle latitude areas (in particular in the ECAC area, European Civil Aviation Conference) all methods provide good results and no interpolation method is found to be significantly better than the other (see [START_REF] Alleau | Comparative Ionosphere Electron Content Estimation Method in SBAS Performances[END_REF]). In particular there is no specific and local geographic pattern where one method provides better results in term of accuracy.

The context of the equatorial area is very different due to rapid changes in the ionosphere dynamics. Some physical events, as plasma bubbles, may appear in regions not observed by measurements, leading wide-averaging techniques to produce GIVD errors that may not be covered by the GIVE. These considerations led TAS to study adaptive mesh, which vertices remain close to the data each time.

THEORETICAL DESCRIPTION OF SELF ORGANIZING MAPS

The idea of having an adaptive mesh is to have its nodes constantly moving according to the inputs, in our case the Ionospheric Pierce Points. This way, the mesh is refined in regions where IPP density is high, and looser in the regions where density is lower. This allows having a much finer spatial definition where there is more information. The approach we used to create an adaptive mesh is to use Self Organizing Maps (SOM).

Self-Organizing Maps have been introduced by T.

Kohonen in [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF], and are a kind of artificial neural network. Artificial neural networks create nonlinear mapping between inputs and outputs. They are made of a set of units, or neurons, who are given a weight (a scalar or a vector). With a series of inputs and associated outputs, they "learn" how to map the two by readjusting their weights. Artificial neural networks can be divided in two kinds: with supervised learning, in which the network must first be "taught" with a training set of input/outputs before it is able to tackle real problems, and unsupervised learning, in which the neural network adapts itself continually without supervision. The Self Organizing Maps belong to the latter.

More precisely, the Self Organizing Maps are made of a set of neurons with connections between one another, thus forming a mesh (fig. 1).

Each neuron is assigned a weight, also called a synaptic vector, which dimension is the same as the inputs.

Figure 1: Representation of a Self-Organizing Map neurons forming a rectangular mesh

When an input data occurs, a competitive learning phase is performed, in which the neuron most similar to the input point is searched. The similarity criterion is usually defined as the Euclidean distance between the input point and the synaptic vector. The most similar neuron is called the Best Matching Unit (BMU). The neurons synaptic vectors are then updated to match more closely the input point. The update amplitude decreases with the distance between a neuron and the Best Matching Unit, therefore a neighborhood function is defined on the network, and is usually the Euclidean distance between synaptic vectors. The update formula for a synaptic vector is:

← + , -
with w i the synaptic vector of neuron i, I the input data vector, and theta the neighborhood function between neuron i and the best matching unit s. Several kinds of neighboring functions can be defined. In the original version of SOMs, every neuron is updated and the neighboring function decreases exponentially with the distance. Other versions of SOMs work more locally, by updating only the Best Matching Unit and its neighbors, with an update amplitude greater for the BMU than for its neighbors.

With this process iterated over every input data, the Self Organizing Map neurons are updated so that they are more resembling to the input data that were fed to the network. The network is thus deformed to match the input data.

Self-Organizing Maps can thus be seen as a way of interpolating data by mapping a continuous input space to a discrete space made by the network neurons. To each vector in the continuous input space, an interpolated value is defined by the synaptic vector of the Best Matching Unit. The learning phase of the network aims at having the best interpolation possible.

Of course, it is possible to have a more sophisticated interpolation formula with a combination of the Best Matching Unit and its neighbors. The synaptic vector update formula must be adapted to match the interpolation formula. This is the base principle of Continuous Self Organizing Maps [START_REF] Aupetit | Function Approximation with Continuous Self-Organizing Maps using Neighboring Influence Interpolation[END_REF]. In this version, for each vector in the input space, a set of neurons made of the best matching and some of its neighbors are "activated" and given a weight depending on their distance to the input vector, so that the interpolation is a weighted mean of the surrounding neurons, giving a continuous interpolating function.

We have applied the Self Organizing Maps to the ionosphere estimation problem, by means of an adaptation of the TRIN grid built as a polyhedron that discretizes the thin ionosphere layer.

The construction of this grid starts with one of the five Plato solid, the icosahedron (20 regular triangular facets, 12 nodes), that is discretized by considering a new point in the middle of each edge and forming new sub triangles (four new in each previous facet). This process leads to a more complex polyhedron, automatically not regular, and provides:

2 + 10 • 4 grid points (node), and 20 • 4 triangular facets.

where L is the level of discretisation. We started from the current TRIN grid, which level is 3.

The TRIN vertices will be the neurons of the Self Organizing Map, with the edges defining the neighboring relationship of the neurons.

Figure 2 : Original regular mesh over the globe

The synaptic vectors for each neuron is made of the position of the vertices on the sphere. According to the location of the Ionospheric Pierce Points, the synaptic vectors will be updated to look more like the IPP, which means the TRIN vertices will be moved to fit the IPP distribution. It is then expected that the TRIN nodes will be more concentrated where the IPP density is high, giving a better spatial resolution in the TEC estimation. The nodes will conversely be more spaced where the IPP are rare, but since there is much less information in these regions, a reliable finer spatial resolution cannot be obtained anyway.

As for the TEC value estimation, two options are possible:

1. The TEC value is part of the synaptic vector and is estimated using the same learning phase as for the node locations on the sphere. There are many different adaptation formulas for , , that make sense, we tested many of them and the best one according to GIVD accuracy for our type of data is the following:

, = | | ! < 2 , = . $ %.| | &'ℎ )
With * % being the distance between the neuron and the IPP , and +)) is the estimation of the error at {-.// , 0 .// , 1 .// } : the difference between the input 3+4 .// and the interpolation. The variability of the data is included with ! being the standard deviation to reduce the impact of an error in case of a great variability.

2. The TEC value is not part of the synaptic vector, and is estimated another way. In our case, it is estimated with a Kalman filter just as with a static mesh. The one used in EGNOS considers the state of the synaptic weights 5 678 = 5 6 and the covariance matrix 9 678 = 9 6 + : 6 with : 6 being the process noise matrix. In this case, since the mesh is timeevolving, we have to update the node synaptic values by interpolating them.

; 5 678 = <5 6 9 678 = <9 6 < + : 6 With A being the matrix that has in each line, the weights of the different nodes of the mesh for interpolating the TEC value at its new position. The synaptic vector only contains nodes positions.

Both methods have been tested. A high level analysis of each method is given below.

For method 1, the TEC is treated in the same way as the node locations, the Self Organizing Map learns the IPP distribution along with the IPP TEC values. All of it is done with the same process and does not require an additional method to estimate the TEC values. However, in this problem, we deal with uncertainties in the input data, which come from measurements flawed with errors. Also, an SBAS needs to provide users with a confidence interval (GIVE) on the TEC value estimated on the IGP. It is therefore necessary to have a measure of the estimation quality. This is the main weakness of the Self Organizing Maps, as there is no clear way to process the input variances in the learning phase, nor to have an estimation variance after the learning phase has updated the synaptic vectors. Several attempts have been made to deal with such problems, giving another unsupervised training neural networks, called the Generative Topographic Mapping (GTM) [START_REF] Kiviluoto | S-Map: A network with a simple self-organization algorithm for generative topographic mappings[END_REF]. However, even with GTM, dealing with input data variances is not straightforward, and still does not give an output estimation variance. This method shall be called "SOM" below.

For method 2, two processes need to be implemented and interlaced. The learning phase of the SOM moves the TRIN nodes on the sphere, and then an iteration of the Kalman filter is performed to estimate the TEC value on the TRIN nodes. Of course, the propagation step of the Kalman filter must take into account the node displacements. That being said, the Kalman filter allows taking the IPP TEC variances into account in a straightforward way, and also provides an output TEC variance. This allows the TEC estimation to be optimal in the sense of the Kalman filter optimality properties, thus minimizing the GIVD error. This method shall be called "Adaptive Kalman" below.

RESULTS AND DISCUSSION

As the concerns of this study is the equatorial zone, the experimentations have been made above the geographic zone delimited by [-10 deg, +30 deg] in latitude and [-25 deg, +30 deg] in longitude that corresponds to ASECNA service area (air navigation African agency). The results have been obtained using the data coming from both the SAGAIE network and IGS stations located inside the ASECNA area. The SAGAIE network is a set of GNSS stations that cover the West African Region [START_REF] Secretan | SAGAIE: A Network of GNSS station to investigate the ionosphere behaviour in sub-Saharan regions[END_REF]. The stations currently deployed are all installed on major airports of the ASECNA area: Dakar (1), Lomé (2), Ouagadougou (3), Douala (4), N'Djamena [START_REF] Goodman | Space Weather & Telecommunications[END_REF]. The network considered for the experimentation is presented in figure 3.

Figure 3: GNSS network (SAGAIE + IGS) selected inside ASECNA area

During the experimentation, the methods SOM and Adaptive Kalman have been tested to try to improve the ionosphere estimation already implemented. Every method implemented relied on the same adaptive mesh fitting the distribution of the IPP.

To illustrate the tailoring, the mesh over the equatorial ASECNA zone is presented in Figure 4. We injected 36 hours of sample data in our algorithms and the red points represent the IPP at the last period of time. Compared to the original mesh on Figure 2, there are more nodes on the equatorial part, implying that the triangles are wider in other areas.

The challenge regarding the adaptive mesh is to define a node displacement tuning able to ensure the stability of the mesh all the time. In other words, the move of all nodes shall be controlled to avoid degeneracy of the triangles along time. The tuning presented here preserves the global shape of the initial TRIN, the triangles are only stretched or tightened so that the global geometric properties (isosceles, equilateral) are kept.

In order to compare our maps to a reference, EGNOS RMS map -the Root Mean Square of the error between the estimate and the truth data-is presented on Figure 6. However, EGNOS IGP are not monitored during the whole of the day. The monitoring proportion goes from 3% to 100% for one IGP on the area we study. Therefore, the maps generated with our methods and presented below are only considering the error at the EGNOS monitoring time.

Using the adaptive mesh, the Self Organizing Maps method is analysed first. Figure 7 shows the RMS map over the equatorial IGP on ASECNA using method 1. The SOM adaptation provides a really good fit to the input data as long as there are enough measurements with a small variability. However, several IGP are poorly estimated at some points in the area with a great variability. Since the variability of input data cannot be taken into account in the SOM adaptation, the Kalman filter seemed a better way to estimate the TEC.

The Kalman filter offers a strength regarding the variability of the data, recognizing the significance of an IPP compared to the past ones for the TEC fitting. Moreover, the Kalman formulation has a selfmanagement of measurement noise and provides automatically a formal covariance for all elements of the state vector. It is thus possible to associate a confidence level for each measurement during the process of the TEC computation. The computation time consumed by the Kalman process is equivalent.

Figure 8 shows the RMS map on ASECNA using the adaptive Kalman filter method.

The results are better than EGNOS using the Adaptive Kalman filter. Indeed, 31 of the 38 IGP on the map are more precisely interpolated on average with the adaptive mesh method. The adaptive mesh method does not provide a GIVE estimation at the moment, so the results are compared to EGNOS GIVD values on EGNOS monitoring period.

In order to compare the GIVD error evolution during the whole period, we pick the 6 most monitored IGPs (between 94% and 100% of monitoring time) and display each method results. Figure 5 shows the GIVD error evolution of EGNOS (red), the SOM (blue) and the Adaptive Kalman filter (green).

CONCLUSIONS AND PERSPECTIVES

Two TEC tailoring methods have been implemented and tested. Both based on adaptive mapping, they provide better GIVD error RMS than the method implemented in EGNOS. First, the Self Organizing Map using an exponential function practically always provides a very good fit to the input data. Nevertheless, GIVD error spikes can occasionally occur on a few IGP when input measurements are very noisy.

In order to take into account the error inherent to a measurement, an adaptive Kalman filter is best suited for the task. This second method also provides better GIVD errors than EGNOS and avoid GIVD error spikes thanks to its capacity to deal with the covariance contrarily to the SOM method.

Therefore a new method to improve GIVD accuracy is proposed. The accuracy improvement allows to reduce the GIVE while maintaining the same integrity margins hence improving availability. It is shown that this method is effective on an equatorial region. These results support the idea that equatorial SBAS is possible.

1 3 4 Figure 4 :

 144 Figure 4 : Mesh (in blue) adapted to 36h of input data, and IPP (in red) at a given date

Figure 5 :

 5 Figure 5 : GIVD error evolution on the 6 most monitored IGP Below is a table with the GIVD error RMS for these IGP, with the three methods: IGP(10,10) IGP(10,5) IGP(10,0) EGNOS 1.35 m 1.57 m 1.80 m SOM 0.94 m 0.92 m 1.02 m Kalman 0.75 m 0.86 m 0.91 m IGP(5,10) IGP(5,5) IGP(5,0) EGNOS 1.67 m 1.66 m 1.63 m SOM 1.38 m 1.39 m 1.30 m Kalman 1.33 m 1.32 m 1.23 m

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6 : EGNOS RMS map

ACKNOWLEDGMENTS

The authors would like to thank Michel Monnerat at TAS, Hugues Secretan at CNES and Roland Kameni at ASECNA for the availability of SAGAIE data used for the experimentations.

is currently a Thales Alenia Space trainee. Soon to be Thales Alenia Space on the EGNOS project. He is a specialist in Navigation Performances, working at Thales Alenia Space on the EGNOS Project as Algorithm Navigation Design Authority.


Dominik Kortschak

Esterina Masiello

Pierre Ribereau

Ruin problem in a changing environment and application to the cost of climate change for an insurance company

Keywords: Non homogeneous Poisson process, ruin probability, asymptotics, cost of climate change

In this paper we obtain asymptotics for ruin probability in a risk model where claim size distribution as well as claim frequency change over time. This is a way to take into account observed and/or projected changes, due to climate change, in some specific weather-related events like tropical storms for instance. Some examples will be presented in order to illustrate the theory and start a discussion on the possible cost of climate change for an insurance company who wants to remain financially solvent.

Introduction

Climate change is of great concern for insurers because of increasing in frequency and intensity of extreme weather events. Tropical cyclones, for example, represent a significant threat to coastal population and infrastructure as well as marine interests (shipping and offshore activities for example). More frequent and/or intensive events would affect insurance systems because of increasing losses also related to economic growth, a greater concentration of people and wealth in periled areas and rising insurance penetration. Detection of long term past trends in measures of tropical cyclone activity is constrained by data availability, quantity and quality: for example, historical tropical cyclone records are known to be heterogeneous owing to changes in measurement practices over time. Nevertheless, the IPCC (Intergovernmental Panel on Climate Change) Fifth Assessment Report (AR5, WG1)( [START_REF][END_REF]) highlighted that there is very strong evidence that frequency and intensity of the strongest tropical cyclones have increased in the North Atlantic basin since the 1970s. As to the future, a broad range of modeling studies, resumed in the IPCC Special Report on Extreme Events ( [START_REF] Ipcc | Managing the risks of extreme events and disasters to advance climate change adaptation[END_REF]), project decreases or no change in the overall global tropical cyclone frequency 1 but also project a substantial increase in the frequency of the most intense storms in some ocean basins by the end of 21st century. This is clearly worrying as the strongest storms are generally responsible for the majority of damage. IPCC Special Report on Extreme Events ( [START_REF] Ipcc | Managing the risks of extreme events and disasters to advance climate change adaptation[END_REF]) also announces, because of climate change, a likely increase in mean maximum wind speed in tropical cyclones as well as in heavy rainfalls associated with tropical cyclones which is expected to intensify cyclone's impacts.

According to the IPCC Fifth Assessment Report ( [START_REF][END_REF]), direct and insured losses from weather-related disasters have substantially increased in recent decades. With respect to the future, the projected increase in the frequency of most intensive tropical cyclones will result in higher direct economic losses and loss variability. This will challenge insurance systems to offer coverage for premiums which are still affordable while at the same time requiring more risk-based capital (AR5, WG2). It is clear that if insurance coverage is to be maintained, insurers would need more risk-based capital to indemnify catastrophic losses and remain financially solvent.

Apart from climate change related trends, tropical storms and hurricanes also exhibit very specific seasonal patterns. The peak season for these weather events in the Southern Hemisphere is January to March while in the Northern Hemisphere, most tropical storms and hurricanes develop in June to November. Figure 1 describes the yearly estimated probability of hurricane occurrence in the Florida gulf for the period 1851 to 2013. Clearly there is a pattern in such a kind of climatic event: tropical storms and hurricanes rarely form during the months December to April in the North Atlantic and North Pacific basins. Indeed, Atlantic basin tropical cyclone activity either before 1 June or after 30 November is nearly negligible. Both storms and hurricanes show a strong maximum in mid-September with most cyclones occurring between 1 August and 31 October, when hurricane season reaches its peak. The same kind of seasonal pattern may also be observed for losses: in general, most expensive hurricanes will occur when the probability of occurrence of the hurricane is higher. For example, Katrina (125 Milliards dollars) occurred at the end of August, beginning of September (2005), Sandy (68M) and Wilma (29.3M) occurred at the end of October (2012 and 2005), Ike (37.5) occurred at the beginning of September (2008) and so on. Some exceptions are represented by Beryl (148 000 dollars, end of May 2012), Arthur (78 millions, end of May 2008), Olga (45 millions, mid December 2007) and Odette (8 millions, beginning of December 2003) but they were classified as tropical storms and not hurricanes.

In the light of the above comments, we will assume in our paper that not only frequencies but also observed losses are changing over time. For example, losses may exhibit seasonal variations corresponding to the meteorological characteristics of the considered weather events as well as an increasing trend due to climate change. Ignoring changing hazard conditions results in biased estimates of expected loss, loss variability and risk capital requirements and this may affect the solvency of an insurance company. Ruin problems in a periodic environment have been studied in Asmussen and Rolski (1994) [START_REF] Asmussen | Risk theory in a periodic environment: the Cramér-Lundberg approximation and Lundberg's inequality[END_REF] and in [START_REF] Rolski | Stochastic processes for insurance and finance[END_REF] [START_REF] Rolski | Stochastic processes for insurance and finance[END_REF]. A periodic risk process is studied in [START_REF] Chukova | Compound counting process in a periodic random environment[END_REF] [6] from a reliability theory point of view. A risk process with periodic claim intensity is considered by [START_REF] Morales | On a surplus process under a periodic environment: a simulation approach[END_REF] [START_REF] Morales | On a surplus process under a periodic environment: a simulation approach[END_REF] who provides a practical simulation methodology to evaluate ruin probabilities and Lu and Garrido (2005) [START_REF] Lu | Doubly periodic non-homogeneous poisson models for hurricane data[END_REF] who study estimation method.

Lin (2014) [START_REF] Lin | Second order tail behaviour for heavy-tailed sums and their maxima with applications to ruin theory[END_REF] studies second order asymptotic results for the sum of not necessarily identically distributed heavy tailed random variables with applications to ruin theory. A time inhomogeneous risk model has also been studied in Kortschak et al.

(2015) [START_REF] Kortschak | Ruin problems with worsening risks or with infinite mean claims[END_REF] where claim size distribution changes over time in an unfavorable direction but there is no change in claim frequency. In this paper we study the asymptotics of the ruin probability when both the frequency and the severity depend on time. Our aim is to start a discussion on the cost of climate change for an insurance company in a sense that will be specified later in the paper. Finding a model for climate change is not an easy task. Therefore we analyze a simplified insurance portfolio on a long term basis which in some sense resembles climate change. We will assume that the insurance company works on a yearly basis. This means that at the beginning of each year the company sets a premium that is adequate compared to the expected claims for the next period. Further the insurance company has to provide enough solvency capital in order to cover the signed risks and remain solvent.

In the long run we will expect that the paid claims are approximately the same as the expected claims and this is reflected by an increased premium. Since it is natural to think that customers are willing to accept premiums that are related to their expected claim amounts, we will assume that the increase in claims is basically paid by the customers and hence is no real cost for the insurance company. On the other hand, because of the increased risk we will assume that the needed risk capital will increase over time. Now, this increase in capital is basically what the climate change will cost the insurance company.

The paper is organized as follows. In Section 2, we present a risk model that allows for frequency and severity to depend on time. Asymptotics for ruin probability under the proposed model are obtained in Section 3. In Section 4 we present some examples and focuse on the cost of adapting to the impacts of climate change for an insurance company. Conclusions are provided in Section 5.

The risk model

As mentioned above, we want to define a risk model which allows to take into account a changing environment due to climate change. Once the model defined and some results on ruin probability obtained, we will try to define a way to calculate the cost of climate change for the insurance company in a hypothetical situation. Our mathematical framework is the following. Assume that X 1 , X 2 , . . . are independent and identically distributed (iid) regularly varying random variables (rv's) of index α > 1 and distribution function (df) F . The risk process will be defined as

S t = Nt i=1 µ(T i )X i -ct
where N t is a nonhomogeneous Poisson process with intensity λ(t), µ is a bounded function on [0, T ], T i is the time of the i-th claim and c is the premium intensity. Further denote with

Y i = µ(T i )X i the i-th claim size, T i = T i -T i-1 the inter arrival times, µ m = min 0≤t≤T µ(t), µ M = max 0≤t≤T µ(t) and Λ(t) = t 0 λ(x)
dx the hazard function of the process. The function µ(T i ), which is representing a multiplicative factor on the claim size, will be chosen in such a way to ensure a severity changing over time. We are interested in the finite time ruin probability defined by

ψ(u, T ) = P sup 0≤t≤T S t > u
where u is the insurer's initial capital. To evaluate ψ(u, T ) we will use an asymptotic approximation. Note that for fixed N t , Nt i=1 µ(T i )X i has the same distribution as the same sum where the T i 's are replaced by iid random variables (T i ) with density λ(x)/Λ(t). Hence we can deduce a first order asymptotic approximation by standard methods like in [2, section X.4]

ψ(u, T ) ∼ T 0 P(X 1 > u/µ(x))λ(x)dx. (1)
Since this asymptotic does not depend on the premium rate c, Section 3 will be devoted to the calculation of the error rate for this approximation. As we will see (Theorem 3.1), the error rate will depend on the premium rate. Indeed, this is important because we would like to integrate changes in premiums (due to changes in risk) in our calculations of the solvency capital for example (see Section 4). Since the error rate obtained in Theorem 3.1 is not explicitly given, bounds will be derived in Corollary 3.2 and the upper bound, given below, will be used as an approximation for the ruin probability in Section 4, for the applications:

ψ(u, T ) ≈ E T 0 λ(t)F (u/µ(t))dt + f (u) T 0 λ(t) µ(t)E [X 1 ] t 0 µ(y) α λ(y)dy + µ(t) α E [X 1 ] t 0 µ(x)λ(x)dx -ct dt.

Asymptotic Results

In this section we provide the approximation for the error term in the asymptotic approximation [START_REF] Alfieri | Global warming increases the frequency of river floods in europe[END_REF]. Note that similar results can also be found in [START_REF] Borovkov | On probabilities of large deviations for random walks. i. regularly varying distribution tails[END_REF] [START_REF] Borovkov | On probabilities of large deviations for random walks. i. regularly varying distribution tails[END_REF] and Lin (2014) [START_REF] Lin | Second order tail behaviour for heavy-tailed sums and their maxima with applications to ruin theory[END_REF].

Theorem 3.1. Let X 1 , X 2 , . . . be iid random variables with distribution function F and survival function F . F is regularly varying of index α > 1 and has regularly varying density f . Further let 0 < µ m ≤ µ(t) ≤ µ M be a bounded function (on [0, T ]). Then

ψ(u, T ) -E [N T ] T 0 P(X 1 > u/µ(x)) λ(x) Λ(T ) dx ∼ f (u) T 0 λ(t)µ(t) α E [X 1 ] t 0 λ(x)µ(x)dx -ct dt + T 0 λ(t)µ(t)E Nt k=1 µ(T k ) α ∞ sup T k ≤s≤t c(t-s)-(S N t ,-k -S Ns,-k ) µ(t)
F (y) dy dt .

Remark 3.1. Lin (2014) [START_REF] Lin | Second order tail behaviour for heavy-tailed sums and their maxima with applications to ruin theory[END_REF] provides error rates for the probability that the maximum of the cumulative sums of independent but not identically distributed random variables exceeds a threshold u. This problem is quite close to the problem studied in this paper. Hence, we may write the error rate in Theorem 3.1 in a similar way to the rate provided in Lin (2014) [START_REF] Lin | Second order tail behaviour for heavy-tailed sums and their maxima with applications to ruin theory[END_REF]:

f (u) E [X 1 ] E   N T k=1 µ(T k ) α N T i =k µ(T i )   -cE T N T N T k=1 µ(T k ) α + E N T -1 k=1 k l=1 µ(T l ) α 0 ∨ min k<m≤N T c(T m -T k ) - m i=k+1 µ(T i )X i .
Proof. Without loss of generality (w.l.o.g.) we will assume that µ M > 1. Conditioning on the event that ruin occurs at time t we get

ψ(u, T ) = T 0 λ(t)E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0} dt.
To evaluate the above integral note that

E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0} = E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0, Ns i=1 µ(T i )X i ≤u/2} (2)
+ E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0, Ns i=1 µ(T i )X i >u/2} . (3)
At first we analyze the term (2). Note that

1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0, Ns i=1 µ(T i )X i ≤u/2} = 1 { Ns i=1 µ(T i )X i ≤u/2} .
Further a Taylor expansion leads to (0

≤ |ξ| ≤ |ct -Nt i=1 µ(T i )X i |) E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 N t i=1 µ(T i )X i ≤u/2 = F (u/µ(t))P Nt i=1 µ(T i )X i ≤ u/2 + E Nt i=1 µ(T i )X i -ct µ(t) f u -ξ µ(t) 1 N t i=1 µ(T i )X i ≤u/2
.

By Potter bounds ([4]

) we get that there exists an > 0 and a C 1 > 1 so that

F (u/µ(t))P Nt i=1 µ(T i )X i > u/2 ≤ F (u/µ(t))E N T P X 1 > u 2N T µ M ≤ C 1 F (u) 2 µ 2α+ E N 1+α+ T ,
where µ = sup 0≤s≤t µ(t). Since f is regularly varying we get for some C 2 > 0 the bound

Nt i=1 µ(T i )X i -ct µ(t) f u -ξ µ(t) 1 N t i=1 µ(T i )X i ≤u/2 ≤ C 2 f (u)(2µ M ) 1+α+ Nt i=1 µ(T i )X i -ct µ(t)
.

By dominated convergence we get that

E T 0 λ(t)E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 N t i=1 µ(T i )X i ≤u+ct dt = E T 0 λ(t)F (u/µ(t))dt + f (u) T 0 λ(t)µ(t) α E Nt i=1 µ(T i )X i -ct dt + o f (u)) + O(F (u) 2 .
Next we have to provide the asymptotic expansion for (3). Note that for i = j

P X i > 1 4µ M N T , X j > 1 4µ M N T ≤ C 3 E N 2α+ T µ 2α+ M F (u) 2 .
Denote with S Ns,-k = Ns i=1,i =k µ(T i )X i and

A k = Nt i=1 µ(T i )X i > u/2, S Nt,-k ≤ u/4 .
Then, by dominated convergence, we get

E F u + ct -Nt i=1 µ(T i )X i µ(t) 1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0, Ns i=1 µ(T i )X i >u/2} = E Nt k=1 F u + ct -Nt i=1 µ(T i )X i µ(t) 1 {inf0≤s<t u+cs-Ns i=1 µ(T i )X i >0,A k } + O(F (u) 2 ) = E   Nt k=1 inf T k ≤s≤t u+cs-S Ns,-k µ(T k ) u/2-S N t ,-k µ(T k ) F u + ct -S Nt,-k -µ(T k )y µ(t) f (y)dy1 {SN t ,-k ≤u/4}   + O(F (u) 2 ) = E Nt k=1 1 µ(T k ) u/2+ct sup T k ≤s≤t c(t-s)-(S Ns,-k -S Ns,-k ) F y µ(t) f u + ct -S Nt,-k + x µ(T k ) dy1 {SN t ,-k ≤u/4} + O(F (u) 2 ) = f (u)E Nt k=1 µ(T k ) α µ(t) ∞ sup T k ≤s≤t c(t-s)-(S Ns,-k -S Ns,-k ) µ(t) F (y) dy + o(f (u)) + O(F (u) 2 ).
Corollary 3.2. Let X 1 , X 2 , . . . be iid random variables with distribution function F that is regularly varying of index α > 1 and has regularly varying density f . Further let 0 < µ m ≤ µ(t) ≤ µ M be a bounded function (on [0, T ]) then

ψ(u, T ) -E [N T ] T 0 P(X 1 > u/µ(x)) λ(x) Λ(T ) dx f (u) T 0 λ(t) µ(t)E [X 1 ] t 0 µ(y) α λ(y)dy + µ(t) α E [X 1 ] t 0 µ(x)λ(x)dx -ct dt and ψ(u, T ) -E [N T ] T 0 P(X 1 > u/µ(x)) λ(x) Λ(T ) dx f (u) T 0 λ(t)µ(t) ct µ(t) F ct µ(t) + E X 1 1 X 1 > ct µ(t) t 0 µ(y) α λ(y)dydt + f (u) T 0 λ(t)µ(t) α E [X 1 ] t 0 µ(x)λ(x)dx -ct dt.
Proof. Just note that 0 ≤ sup

T k ≤s≤t c(t -s) -(S Ns,-k -S Ns,-k ) µ(t) ≤ ct µ(t) and ∞ a F (x)dx = aF (a) + E X 1 1 {X 1 >a} .
As already mentioned, the upper bound obtained in Corollary 3.2 will be used to approximate the ruin probability in the following section where two examples will be presented in the purpose of estimating the cost of climate change for an insurance company in a fictitious portfolio.

Some examples

Climate change is changing the risk of the insurance portfolio and, as a consequence, insurance companies have to adjust their solvency capital over time in order to fulfill solvency requirements. In other words, climate change is representing a cost for insurance companies that we would like to evaluate on some examples thanks to the theoretical results on ruin probability obtained in Section 3.

Example 1. We consider a model where the X i 's are Pareto distributed rv's whose survival function can be written as

F (x) = (1 + x) -α .
We take α = 2.5. We assume that the intensity of the non homogeneous Poisson process λ(t) is equal to the multiplicative factor µ(t) and exhibits a cyclic behavior (corresponding to seasonal patterns) as well as a long term evolutionary trend. This trend should represent an increasing risk because of climate change. In practice the trend should be chosen on expected change in risk because of climate change. We take

µ(t) = λ(t) = 1 4 + 1 10 sin 2π t T + 0.075 120 t T
2 so that the claim size distribution is changing over time, not only the frequency. More precisely, we are representing cases where the higher the frequency is, the higher the losses are. As already explained in the introduction, this kind of behavior is quite common for hurricanes for example. We will compare this model with a similar model without a trend where the intensity function includes a trigonometric component but not a polynomial one:

µ * (t) = λ * (t) = 1 4 + 1 10 sin 2π t T .
This may correspond to the case where we admit there is no climate change. See Figure 2 for a plot of µ(t) (respectively λ(t)) and µ * (t) (respectively λ * (t)). In this paper we assume that the insurance company works on a yearly basis. Let us denote by T 0 the starting time of the year. We take T = 12, corresponding to the number of months in a year. Note that the expected claim in one year is given by

E   N T +T 0 i=N T 0 +1 µ(T 0 + T i )X i   = E [X 1 ] T +T 0 T 0 µ(T 0 + t)λ(T 0 + t)dt.
See Figure 3 for a plot of the expected claim amount as a function of time both with and without a trend. Obviously, because of the assumptions in our risk model, we observe an increase in the expected claim amount over the time when we consider a trend.

We assume that the insurance company is collecting risk-commensurate premiums at a constant rate c per unit time where and θ is the relative safety loading. We take θ = 10% in our example. Clearly, premiums are adjusted at the beginning of every year and, because of climate change, are increasing over time. We assume that customers are ready to accept to pay increasing premiums as premiums are commensurate to risk. But increasing claim amounts also demand to adjust risk capital every beginning of the year.

c = (1 + θ) 1 T E   N T +T 0 i=N T 0 +1 µ(T 0 + T i )X i   0
The needed solvency capital u for a year is estimated so that the finite time ruin probability is small enough, i.e. ψ(u, T 0 , T 0 + T ) = 1/200. As an approximation for the ruin probability we will use the asymptotic upper bound given in Corollary 3.2. To evaluate the needed capital for later years we just start the process at the starting time of the considered year and recalculate the premium. This provides us with a series of needed solvency capital for the future. The resulting solvency capitals for the different years are represented in Figure 4. Moreover, in Figure 5 we plot the relative solvency capital, i.e the needed solvency capital over the expected claim amount. We can observe that, in case of trend, the relative solvency capital is decreasing over time. A possible explanation is that with more claims there is a bigger diversification effect that reduces the relative solvency capital.

To sum up, an insurance company has to provide solvency capital every year and, in a model with trend, the solvency capital is increasing every year. Clearly, insurance companies have to follow. The increase of the solvency capital can only partly be financed by the earnings through safety loadings, hence extra money has to be inserted. This needed extra money can be seen as the true cost of the trend for an insurance company and will be denoted as cost of trend in the rest of the The accumulated cost of trend is represented in Figure 6 (black line) together with the solvency capital without trend (red line) and the cost of trend (green line). We can observe that the cost of trend is increasing over time but stays smaller than the solvency capital without a trend. Finally we want to relate the yearly cost of trend to the expected claim amount. In Figure 7 we provide the yearly cost of trend over the expected claim amount together with the collected premiums over the expected claim amount. We can observe that the relative cost increases in the beginning but decreases after approximately 15 years. The cost of trend can reach up to 40% of the expected claim amount which is much higher than the chosen safety loading (red line). So we can conclude that, in this specific example, the cost of trend is a significant cost for insurance companies.

Example 2 (about flood risk). In this example we will focus on flood hazard in France. Among flood hazards, river floods are considered one of the most important natural disasters in Europe whose impact can be devastating. We will use results by [START_REF] Alfieri | Global warming increases the frequency of river floods in europe[END_REF] [START_REF] Alfieri | Global warming increases the frequency of river floods in europe[END_REF] to obtain projections of changes in the frequency of flood hazard in France through the current century. Since we did not find usable informations on projected changes in claim size related to flood hazard in France, we decided to assume that the distribution of claim sizes will stay constant over time. On the other hand, the claim frequency will be supposed to increase in a way that we are going to precise. [START_REF] Alfieri | Global warming increases the frequency of river floods in europe[END_REF] [START_REF] Alfieri | Global warming increases the frequency of river floods in europe[END_REF], investigate changes in the frequency of extreme events with return periods equal or larger than 100 years, in three time We will use a linear approximation of the frequencies in Table 1 in order to obtain an intensity function, λ(t) = -0.463 + 0.000238 * (1990 + t), corresponding to the projected increase in frequencies. Again, we assume that the X i 's are Pareto distributed rv's with α = 2.5. Since we are only considering events with a return period of at least 100 years we will use the conditional distribution given that we are bigger than the 100 year event i.e we will use x > x 0 with x 0 = 0.01 -1/2.5 -1 ≈ 5.3. The survival function of X can then be written as:

F (x) = 100(1 + x) -2.5 , x > x 0 .
Moreover we will assume that in the model without trend, both the intensity function λ(t) and the multiplicative function µ(t) are constant. More precisely, λ * (t) = 0.01 and µ * (t) = 1. Remark that we also assume µ(t) = µ * (t) so that claim size is not time dependent in this example. We will use T = 1. Other assumptions required for the calculation of quantities of interest, are the same as in Example 1.

Because of the assumptions in our example, the solvency capital needed at the beginning of every year is increasing over time in case of trend (see Figure 8). On the other hand, the relative solvency capital (solvency capital needed at the beginning of every year over the expected claim amount) in case of trend is decreasing because of possible diversification effects (see Figure 9). Again, the accumulated cost of trend is increasing but it is, in this example, smaller than the needed solvency capital (see Figure 10). Further, the cost of trend is slightly decreasing. Finally, in Figure 11, we can observe that the relative cost is increasing at the beginning but then immediately decreasing. Anyhow, it is higher than the safety loading, so representing a significant cost for the insurance company.

Conclusions

In this paper a generalization of the classical risk model has been proposed to better represent more realistic cases where weather-related claims not only present a seasonal pattern but are also affected by climate change. The model we propose is essentially a non-homogeneous Poisson process with severity depending on time. This allows for environments where the claim size is higher or lower depending on the season. For example, we noted that the most expensive hurricanes occur in the peak season. Our model allows to take into account this kind of trend in weather events as well as other kinds of events. Under the proposed risk model, we obtain asymptotics for finite time ruin probability. Theoretical results are then used to illustrate on some examples how much climate change would cost to an insurance company who will experience events, like tropical storms or floods, which are more

Figure 1 :

 1 Figure 1: Estimated probability of hurricane occurrence in the gulf of Florida in one year conditionally to the occurrence of the hurricane (1851-2013).

Figure 2 :

 2 Figure 2: λ(t) (respectively µ(t)) with (in black) and without (in red) a trend.

Figure 3 :

 3 Figure 3: Average claim amount for the different years with (in black) and without (in red) a trend.

Figure 4 :

 4 Figure 4: Solvency capital needed at the beginning of every year with (in black) and without (in red) a trend.

Figure 5 :

 5 Figure 5: Solvency capital needed at the beginning of every year over the expected claim amount, with (in black) and without (in red) a trend.

Figure 6 :

 6 Figure 6: Accumulated cost of trend (black line), solvency capital without trend (red line), and the cost of trend (green line).

Figure 7 :

 7 Figure 7: Yearly cost of trend over the expected claim amount (black line) and collected premiums over expected claim amount (red line).

Figure 8 :

 8 Figure 8: Solvency capital needed at the beginning of every year with (in black) and without (in red) a trend.

Figure 9 :

 9 Figure 9: Solvency capital needed at the beginning of every year over the expected claim amount, with (in black) and without (in red) a trend.

Figure 10 :

 10 Figure 10: Accumulated cost of trend (black line), solvency capital without trend (red line), and the cost of trend (green line).

Figure 11 :

 11 Figure 11: Yearly cost of trend over the expected claim amount (black line) and collected premiums over expected claim amount (red line).

Table 2

 2

	in Alfieri et al. (2015)

Table 1 :

 1 Mean annual exceedance frequency of the 100-year return period peak flow for France.

and more frequent and whose consequences are heavier and heavier. In our simplified portfolios, the cost of climate change seems to be significant for insurance companies. This field seems worth exploring on more complex models.