# Modeling the Tropical Cyclone risk in the North-West Pacific Basin

Octave Mesnard, Mathis Joffrain

<figure><img src="image_1.png" /></figure>

## To cite this version:

Octave Mesnard, Mathis Joffrain. Modeling the Tropical Cyclone risk in the North-West Pacific Basin. 2023. hal-04239442

HAL Id: hal-04239442

https://hal.science/hal-04239442v1

Submitted on 26 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Actuarial thesis :
Modeling the Tropical Cyclone risk in the
North-West Pacific Basin

Author :
Octave MESNARD

Supervisor :
Mathis JOFFRAIN

ESSEC professor :
Marie KRATZ
<figure><img src="image_3.png" /></figure>
# Résumé

Les catastrophes naturelles ont un impact majeur sur la société et l'économie. Parmi ces catastrophes, les cyclones tropicaux constituent l'un des risques principaux auxquels sont confrontées les communautés côtières.

Dans le bassin Pacifique Nord-Ouest, les cyclones tropicaux ont causé 36,5 milliards de dollars de pertes économiques pour la seule année 2020. Pour comprendre et quantifier ces risques, les assureurs et réassureurs non-vie s'efforcent d'utiliser les meilleurs outils qui sont, à ce jour, des logiciels commerciaux. Dans certains cas, des modèles internes peuvent être développés à l'aide des dernières avancées scientifiques et des données les plus récentes.

Cette étude développe un modèle interne de catastrophe qui évalue les dommages aux bâtiments causés par les cyclones tropicaux : Le chapitre 1 présente les caractéristiques des cyclones tropicaux ; les chapitres 2 et 3 examinent les trajectoires des cyclones et la modélisation du champ de vent à l'aide de techniques d'apprentissage automatique et le chapitre 4 étudie les impacts financiers et la prise en compte du changement climatique dans la modélisation.

**Mots-clés** : Cyclone tropical, Catastrophe naturelle, Apprentissage automatique, Vulnérabilité, Changement climatique
Abstract

Natural disasters have a big impact on society and the economy. Among these disasters, tropical cyclones are one of the major risks coastal communities face.

In the Western North Pacific (WNP) basin, Tropical Cyclones caused USD 36.5 Bn in economic losses in 2020 alone [14]. To understand and quantify these risks, P&C insurers and reinsurers have been striving to use the best tools and to date, these have been commercial. In certain cases, internal models can be developed with the latest state-of-the-art science and data.

This study presents an internal catastrophe model that assesses the building damage cyclonic winds cause : Chapter 1 characterizes the wind hazard ; Chapters 2 and 3 examine cyclone trajectories and windfield modeling with machine learning techniques ; and Chapter 4 studies the financial impacts and climate change modeling per se.

**Keywords** : Tropical cyclone, Natural catastrophe, Machine learning, Vulnerability, Climate change
Acknowledgements

First of all, I would like to thank the teams from Group Risk Management entity of AXA GIE for giving
me the opportunity to write this study. They form a skillfull and fulfilling work environment to carry out these
projects, to help to work with efficient computing equipment, and more generally to discover the modeling of
all types of perils.

I would particularly like to thank Mathis Joffrain, who was my tutor at AXA and who, among other things, drives the modelling of Tropical cyclone forward within the Group Risk Management. First of all, I thank him for his support and management during my experience at AXA GRM. Second, for its accompaniment throughout the writing of the thesis by giving valuable advice on the tone and form to adopt as well as to guide me on the perspective of this work.

Then, I want to acknowledge the work of all modellers who previously contributed on the modelling of
Tropical Cyclone at AXA GRM and helped to start from a solid basis for this thesis. They include, Edern,
Madeleine, Mathis, Hugo, Anh and all the other ones I did not have a chance to meet or hear about.

I would particularly like to thank Marie Kratz, who, in addition to being my academic tutor for the writing of this dissertation, is the director of the actuarial track at ESSEC. I want to underline the importance of this course as a student and the incredible luck of having been able to be trained through dense and instructive courses at both ESSEC and ISUP.

Finally, I would like to extend my thanks to Marianne Le Guen and Timothée Mesnard who carefully
reviewed the work and gave me a lot of insights on the understanding of the study for outsiders.
List of Acronyms

**TC** Tropical Cyclone

**YLT** Year Loss Table

**MPI** Maximal Potential Intensity

**OEP** Occurrence Exceedance Probability

**QRF** Quantile Regression Forest

**PCA** Principal Component Analysis

**FEMA** Federal Emergency Management Agency

**SST** Sea Surface Temperature

**NWP** North West Pacific

**DR** Damage Ratio

**RP** Return Period
Table des matières

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Introduction to Tropical Cyclones</td><td>8</td></tr><tr><td>1.1</td><td>Natural Catastrophe modeling at AXA</td><td>8</td></tr><tr><td>1.1.1</td><td>Organization and missions</td><td>8</td></tr><tr><td>1.1.2</td><td>External catastrophe models</td><td>8</td></tr><tr><td>1.1.3</td><td>Internal catastrophe models</td><td>10</td></tr><tr><td>1.2</td><td>Description of Cyclonic events</td><td>11</td></tr><tr><td>1.2.1</td><td>Introduction to cyclonic events</td><td>11</td></tr><tr><td>1.2.2</td><td>Cyclone's life cycle</td><td>11</td></tr><tr><td>1.2.3</td><td>Cyclone's secondary perils</td><td>12</td></tr><tr><td>1.3</td><td>Data used</td><td>12</td></tr><tr><td>1.4</td><td>Scope</td><td>12</td></tr><tr><td>1.4.1</td><td>Spatial</td><td>13</td></tr><tr><td>1.4.2</td><td>Temporal</td><td>13</td></tr><tr><td>1.5</td><td>Modelled parameters</td><td>13</td></tr><tr><td>1.5.1</td><td>Location parameters</td><td>13</td></tr><tr><td>1.5.2</td><td>Intensity parameters</td><td>15</td></tr><tr><td>1.5.3</td><td>Correlation matrix between variables</td><td>17</td></tr><tr><td>1.6</td><td>Limits</td><td>18</td></tr><tr><td>1.6.1</td><td>Limitation in the data</td><td>18</td></tr><tr><td>1.6.2</td><td>Limitations in the approach</td><td>20</td></tr><tr><td>1.7</td><td>Key takeaways</td><td>20</td></tr><tr><td>2</td><td>Track Generation</td><td>22</td></tr><tr><td>2.1</td><td>TC events frequency</td><td>22</td></tr><tr><td>2.1.1</td><td>Block bootstrap methodology</td><td>22</td></tr><tr><td>2.1.2</td><td>Time dependence structure and size of the bootstrap blocks</td><td>23</td></tr><tr><td>2.1.3</td><td>Frequency validation</td><td>24</td></tr><tr><td>2.2</td><td>TC Track generation : generation of location parameters</td><td>25</td></tr><tr><td>2.2.1</td><td>Initiation</td><td>25</td></tr><tr><td>2.2.2</td><td>Track propagation</td><td>28</td></tr><tr><td>2.2.3</td><td>Termination</td><td>31</td></tr><tr><td>2.2.4</td><td>Validation</td><td>31</td></tr><tr><td>2.3</td><td>TC Track generation : generation of intensity parameters</td><td>34</td></tr><tr><td>2.3.1</td><td>Maximum sustained wind speed : MWS</td><td>34</td></tr><tr><td>2.3.2</td><td>Central pressure : CP</td><td>39</td></tr><tr><td>2.3.3</td><td>Radius of maximum winds RMAX</td><td>39</td></tr><tr><td>2.3.4</td><td>Angle of maximum winds : AMAX</td><td>44</td></tr><tr><td>2.4</td><td>Tracks simulation procedure</td><td>44</td></tr><tr><td>2.4.1</td><td>Procedure</td><td>44</td></tr><tr><td>2.4.2</td><td>Remarks on the procedure</td><td>45</td></tr><tr><td>2.5</td><td>Limits</td><td>45</td></tr><tr><td>2.5.1</td><td>Historical based approach for the track generation</td><td>45</td></tr><tr><td>2.5.2</td><td>Underestimation of the MWS</td><td>45</td></tr><tr><td>2.6</td><td>Key takeaways</td><td>46</td></tr></tbody></table>
3 Footprint Generation 47

3.1 Prediction of windspeed fields : framework and challenges 47

3.1.1 Problem framework 47

3.1.2 Challenges 48

3.1.3 Data used 48

3.2 Dimension reduction 48

3.3 Prediction of reconstructed windspeed fields 50

3.3.1 Reduced problem framework 50

3.3.2 Windfields reconstruction 51

3.3.3 From windfield to wind footprint generation 52

3.4 Adjustment of wind footprints 53

3.4.1 Roughness effect 53

3.4.2 Gust effect 54

3.5 Footprints validation 55

3.6 Limits 57

3.6.1 Lack of observation data for wind footprint validation 57

3.6.2 Methodology for gust and roughness effect evaluation 57

3.7 Key takeaways 57

4 Insured Losses Generation 58

4.1 Exposure data 58

4.2 Damage state and damage ratio 58

4.2.1 Damage states 59

4.2.2 Damage ratios 59

4.3 Fragility curves 59

4.4 Adapation of fragility curves 60

4.5 Simulation procedure for Year Loss Table (YLT) and Occurrence Exceedance Probability (OEP) derivation 62

4.6 Towards a framework for climate change impact assessment 63

4.7 Limits 65

4.8 Key takeaways 65

5 Conclusion and Perspectives 67
Introduction

The North West Pacific (NWP) basin has the world highest annual frequency of tropical cyclones. [3] These
are a major challenge for insurance companies covering property damage in the region. Cyclonic events can cause
high economic and insured losses when hitting the land with a certain intensity. The 2013 Haiyan typhoon struck
the Philippines, Vietnam, Laos and the south of China, resulting in 13 billion dollars of economic losses.[1]

For P&C insurers highly exposed in Asia, such as AXA, understanding and modeling tropical cyclone risk in this region is crucial. This risk can be modelled through commercial catastrophe models that offer a complete well-established range of solutions since the 90s. AXA uses these model results for yearly capital requirement evaluations and reinsurance calibrations.

However, since the 2010s, the wave of open data and open science have made developing « in house » catastrophe models easier. The scientific community per se has notably extensively shared databases and open-source libraries to provide a reliable validated basis for modeling tropical cyclone risk.

This study (i) describes a model building methodology that capitalizes on these new open data and library
sources and (ii) evaluates it in application to damage risk from cyclonic wind in the Philippines.

The present work provides to the actuarial community a transparent and applicable modeling method for
cyclonic wind damage evaluation, and a major update of the actuarial thesis [9]. It relies on existing developments
but also introduces new methodologies that strongly upgrade AXA catastrophe modelling standards.

**Chapter 1** provides the study's operational context (1.1), a description of tropical cyclone events (1.2), and focuses on the data used (1.3), by introducing the scope (1.4), the main parameters modelled to represent these events (1.5) and the data sources limitations (1.6).

**Chapter 2** describes a methodology to simulate a catalogue of stochastic tracks. It studies (2.1) the frequency of these tracks, the generation of (2.2) location parameters and (2.3) intensity parameters, and (2.4) the complete simulation procedure of simulation of the stochastic tracks catalogue.

**Chapter 3** applies the tracks from Chapter 2 to generate cyclones windspeed footprints. It introduces (3.1) the windspeed prediction problem and its associated challenges, presents (3.2) a solution based on a dimension reduction approach and covers (3.3) the generation of spatial grids of windspeed values called « footprints ». Chapter 3 then examines (3.4) the adjustment of the wind footprints to include terrain effect, (3.5) the validation of the generated footprints and (3.6) the limitations in the approach.

**Chapter 4** uses the windspeed footprints from Chapter 3 to evaluate economic losses associated to an insured portfolio in the Philippines. It describes (4.1) the exposure portfolio and two key concepts for economic losses generation : (4.2) the damage states with their associated damage ratios and (4.3) the fragility curves. Then chapter 4 details (4.5) the insured losses generation procedure and (4.6) the adaptation of this procedure to include climate change impacts.
# 1. Introduction to Tropical Cyclones

Chapter 1 describes the main parameters and quantities to represent Tropical Cyclone risk impact on insurers' portfolios.

## 1.1 Natural Catastrophe modeling at AXA

### 1.1.1 Organization and missions

This study has been conducted at AXA GIE, the headquarters of AXA Group, within P&C Group Risk Management department. The P&C Group Risk Management department has three main functions :

1. Catastrophe and Reinsurance modelling

2. Solvency II capital requirements modelling

3. Underwriting and reserves review

The Catastrophe and Reinsurance modelling team has three main missions :

- Collecting and monitoring AXA Group exposure to natural catastrophe risk

- Calculating group STEC (Short Term Economic Capital) pertaining to Natural Catastrophe risk, gross and net of reinsurance

- Providing modelling results for both internal and external reinsurance arrangements.

### 1.1.2 External catastrophe models

To fulfill these missions, AXA Group has used analyses from the following commercial softwares :

- AIR (www.air-worldwide.com)

- EQE (www.corelogic.com)

- RMS (www.rms.com)

These three companies operate since the 90s and are the leading players in the catastrophe modelling field. The models licensed by the companies have the same underlying structure with the four following modules :

#### Modules

- **Exposure module** is a list of geolocalized insured risks, together with their sum insured values. These sum insured value can be broken down in three coverages : Building, Content and Business Interruption. In this study, only the coverage building will be studied. Further characteristics of the insured risks can be provided, such as structure or business type. Eventually, exposure module contains contractual conditions associated to the listed risks.

- **Hazard module** is a catalogue of simulated natural peril events. These events are spatial grids of physical parameters. These physical parameters represent the intensity of each event, for instance flood depth for flood events, or windspeed for Tropical Cyclone (TC) events. Since these parameters are associated with spatial coordinates, they can intersect with the geolocalized risks from the exposure module.

- **Vulnerability module** is a set of functions that associate physical parameter values to destruction rates values. destruction rates are the ratio of the loss due to the natural peril event on a given risk over the sum insured of that risk
The interaction between these three modules can be summarized as follows :

<figure><img src="image_4.png" /><figcaption>FIGURE 1.1 – Components of a natural catastrophe model</figcaption></figure>

For any risk in the exposure module, intensities values can be derived by associated each events from the hazard module to this risk. Then these intensities can be translated into losses using the vulnerability functions from the vulnerability module.

The application of these modules leads to a standardized output called Year Loss Table (YLT). This output is key and is the main quantity to represent Tropical Cyclone risk impact on insurers' portfolios.

From the YLT, catastrophe modellers can extract many quantities that are crucial for their day to day operations and communications with other teams.

## YLT : Year Loss Table

An example of YLT is provided below :

<table><thead><tr><td>Year</td><td>Event ID</td><td>Loss</td></tr></thead><tbody><tr><td>1</td><td>1</td><td>100</td></tr><tr><td>2</td><td>2</td><td>0</td></tr><tr><td>3</td><td>3</td><td>500</td></tr><tr><td>3</td><td>4</td><td>120</td></tr><tr><td>4</td><td>5</td><td>200</td></tr></tbody></table>

- The column Year refers to a simulation number. It is called year because the simulated distribution represents one year of catastrophic event

- The column Event Id is the list of all simulated catastrophic events from the hazard module. Event ID 3 and 4 are associated to simulation number 3, which means both events occurred during the simulated year number 3.

- The column Loss is the economic impact by each events in euros.

The YLT can be reduced by taking the maximum loss over each year of simulation leading to the following table :
<table><thead><tr><td>Year</td><td>Max Loss</td></tr></thead><tbody><tr><td>1</td><td>100</td></tr><tr><td>2</td><td>0</td></tr><tr><td>3</td><td>500</td></tr><tr><td>4</td><td>200</td></tr></tbody></table>

### OEP: Occurrence exceedance probability

After sorting the previous table according to the loss column, the following table is derived:

<table><thead><tr><td>Exceedance Probability (EP)</td><td>Return Period (RP)</td><td>Loss</td></tr></thead><tbody><tr><td>100%</td><td>1</td><td>0</td></tr><tr><td>75%</td><td>1.33</td><td>100</td></tr><tr><td>50%</td><td>2</td><td>200</td></tr><tr><td>25%</td><td>4</td><td>500</td></tr></tbody></table>

where $EP(x) = \mathbb{P}(Loss > x)$ and $RP = \frac{1}{EP}$.

The simulated distribution in the table above is called the Occurrence Exceedance Probability (OEP), and is denoted by OEP. It can be used for:

- Defining extreme loss scenarios by extracting high quantiles from this distribution

- Defining contractual conditions in reinsurance contracts

The 99.5th quantile of the loss distribution is defined by:

$$VaR_{0.995}(Loss) = inf\{y \in \mathbb{R} : \mathbb{P}[Loss > y] \le 1 - 0.995\}$$

The 99.5% destruction rate of an insured portfolio is then defined by:

$$DRT(0.995) = \frac{VaR_{0.995}(Loss)}{A}$$

where A is the total sum insured of the portfolio.

The destruction rate is usually provided by peril, and reflects the impact by a given natural peril on the portfolio of the insurer.

### 1.1.3 Internal catastrophe models

The YLT mentioned earlier can also be calculated using models developed internally. Since the 2010s, the wave of open data and new science tools have made easier the development of « in house » catastrophe models.

This wave has three main components:

- Open source libraries in R and Python are freely and easily available. They considerably enriched the available tools, be it for machine learning (scikit learn), geodata processing (raster,sf) and statistical analyses.

- The scientific community has extensively shared databases, such as IBTRACS [20], or physical model results such as ERA5. [27]

- Start ups produced databases and research, such as INCYC [30].
In this context, AXA Group has developed a range of internal catastrophe models for a variety of perils and geographies, including Tropical Cyclone in the North West Pacific basin. The goal of this study is to describe a methodology that builds a YLT for Tropical Cyclone risk in the North West Pacific basin.

## 1.2 Description of Cyclonic events

### 1.2.1 Introduction to cyclonic events

Cyclones have been studied from the middle of the 19th century. The sailor Henry Piddington published in 1844 his book "The Horn-book for the Law of Storms for the Indian and China Seas". He is referred to as the founding father of the study of this meteorological phenomenon, to which it gave the name "cyclone" in reference to a snake winding itself or "kyklos" in Greek.

When they form, these events are called « tropical depression ». Under appropriate conditions, they might evolve in tropical storms and eventually being considered as cyclone, their most extreme form. Other terms are used to designate cyclones such as hurricanes (mostly in the US) or typhoon (specifically in Asia). Here, the term TC (Tropical Cyclone) will refer to these events.

### 1.2.2 Cyclone's life cycle

The life of a cyclone takes place in three phases : formation, intensification and dissipation for a lifetime that can last up to 15 days.

#### 1. Formation

Tropical cyclones usually form with :

* Water temperature exceeding 26 deg [16].

* Latitudes over 5 degrees [16].

* Stable wind, in both horizontal and vertical directions [46].

#### 2. Intensification

Conditions for the intensification of a Tropical Cyclone are similar to the formation conditions, and are described in [45]. It is called a rapid intensification when the windspeed increased by over 55 km/h in less then 24h. Hagibis (2019) is an example of this phenomenon.

#### 3. Dissipation

Tropical cyclones usually dissipate when reaching cold water or land [16]. If they survive long enough to move to upper latitudes (over 30 degrees for instance), their structure changes. This transformation is called « extra tropical transition ». For instance, Mawar (2012) had an extra tropical transition phase.

This study will focus on wind subperil. However, coastal flooding and pluvial flooding are consequent of cyclonic events too.
### 1.2.3 Cyclone's secondary perils

Although cyclones are mainly destructive due to strong winds, they can be combined with two other main sources of hazards: pluvial and coastal flooding. These sub-hazards are an integral part of TC risk modelling, even if this study will only focus on the wind component.

**Pluvial flooding** TC events can release massive amounts of rainfall in a very short period of time with two main damage consequences to property: First, the water causes flooding with high velocity flows. Second, the water can infiltrate buildings. For instance, Jebi caused such hazard in Japan, see [6]

**Coastal flooding** A storm surge is an abnormal rise in the sea generated by a TC. It can flood low lying coastal areas. For instance, Hato caused such hazard in Hong-Kong [21].

## 1.3 Data used

**IBTRACS** [20], is a global database for historical TC events, created by the National Centers for environmental information (NOAA). It provides complete and detailed information on all historical TC events across the world, by cross-referencing data recorded in all major meteorological centers. In particular, it compiles location information, such as past TC center s' geographical positions (latitudes and longitudes), and various intensity parameters at 6 hour intervals. IBTRACS heavily relies on observations, and will be mainly used in Chapter 2, for TC track simulation.

**INCYC** [30], is a global database of tropical cyclone Events, created by REASK, a commercial company, and relies on model simulation. In particular, Weather Research and Forecasting (WRF) numerical weather model [7] was used to simulate representative TC events across the world. WRF parameters for the simulation have been reviewed by various academics, including Céline Planche from UCA [48]. This database is formed of spatially regular gridded atmospheric data, which includes location and intensity variables for each events, such as wind or pressure. INCYC will be mainly used in Chapter 3, for TC footprint generations. spatial resolution is at 1km and temporal resolution is at 1h.

## 1.4 Scope
1.4.1 Spatial

The study focuses on the North West Pacific Basin (NWP), which belongs to a frame between 100 to 180 degrees of longitude and 0 to 70 degrees of latitude.

<figure><img src="image_5.png" /><figcaption>FIGURE 1.2 – IBTRACS - Historical events in North West Pacific Basin from 2008</figcaption></figure>

1.4.2 Temporal

The study relies on data from the range 1980-2018. This choice is justified based on IBTRACS [20] and INCYC data [30], and more explanations are given in the Appendix H.

Some stochastic event models such as [4] rely on a larger timeframe, starting in the 1950s, and do not follow
our choice.

1.5 Modelled parameters

The main parameters modelled in the study to characterize TC events are defined below. Some descriptive
statistics based on IBTRACS data are also provided.

1.5.1 Location parameters

*   **Bearing and Bearing changearta** (*resp. BEAR, dBEAR*) in degrees: Bearing is the direction of the TC event at a given position, here defined as the angle between the segment joining its last two positions and a vector pointing to the east. Angles are counted positively counterclockwise, and range in [−180, 180]. This definition is consistent with [8]. Bearing change is the variation of Bearing over six hours.
<figure><img src="image_6.png" /><figcaption>FIGURE 1.3 – Bearing angle illustration</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 1.4 – Bearing and bearing change - degree</figcaption></figure>

The bimodal partitioning of the bearing distribution (Fig. 1g) [20] shows that the two most common directions are either north/northeastern (BEAR ≃ 45 degrees) or south (BEAR ≃ -80 degrees). The bearing change dBEAR has a symmetric distribution centered in 0.

* **Angle of maximum winds (AMAX)**, in *radians* : is the angle between the segment joining the location of the MWS and the location of the CP, and the Bearing. AMAX is measured clockwise from the Bearing.

* **Distance to land (dist2coast)**, in *km* : This feature reflects the distance between the TC center and the nearest coast (negative when inland, positive when offshore).

<figure><img src="image_8.png" /><figcaption>FIGURE 1.5 – Distance to land - km</figcaption></figure>
A large number of points recorded in the IBTRACS database are located for cyclones quite close to the coast. There is also a very low number of negative values : cyclones disappear quite quickly after a landfall.

* **Forward speed and acceleration (resp. $F$, $dF$)**, in $m/s$ and $ms^{-2}$: Forward speed is the displacement speed of a TC. This parameter is derived from the distance covered by a cyclone between two 6h consecutive time-steps. Acceleration is the variation of Forward speed over time.

<figure><img src="image_9.png" /><figcaption>FIGURE 1.6 – Forward speed (in m/s) and acceleration $ms^{-2}$</figcaption></figure>

The distribution of the forward speed ($F$) is compact with a majority of the captured value belonging to the range [10-40] km/h. The acceleration ($dF$) distribution is centered and close to 0 with rare extreme values exceptions, showing that the forward speed evolution is smooth.

## 1.5.2 Intensity parameters

* **Maximum sustained wind speed (MWS)**, in $km/h$, is the maximum of 1 minute - moving averages of sustained wind speed over an observation period, and over the event area. This also can be referred as « rotational wind speed » and is a key intensity parameter to estimate losses from wind hazard.

<figure><img src="image_10.png" /><figcaption>FIGURE 1.7 – Maximum sustained wind speed - $km/h$</figcaption></figure>

The distribution of MWS of cyclones is asymmetric with a majority of records below 100 km/h and rare instances over 200 km/h. The maximal recorded windspeed is 315 km/h and was reached by TC Haiyan in 2013. Only three cyclones exceeded windspeed of 300 km/h since 1980 : TC Haiyan in 2013, TC Meranti in 2016 and TC Halong in 2019.

* **Saffir Simpson Scale** : is an intensity scale based on MWS
<table><caption>TABLE 1.1 – Saffir Simpson Categories</caption><thead><tr><td>Category</td><td>MWS</td></tr></thead><tbody><tr><td>Tropical Depression</td><td>≤ 63 km/h</td></tr><tr><td>Tropical Storm</td><td>63 - 118 km/h</td></tr><tr><td>CAT 1</td><td>119 - 153 km/h</td></tr><tr><td>CAT 2</td><td>154 - 177 km/h</td></tr><tr><td>CAT 3</td><td>178 - 208 km/h</td></tr><tr><td>CAT 4</td><td>209 - 251 km/h</td></tr><tr><td>CAT 5</td><td>≥ 252 km/h</td></tr></tbody></table>

TC events are classified by these seven denominations.

* **Central pressure (CP)**, in hPa, is the minimal pressure reached over the event area.

<figure><img src="image_11.png" /><figcaption>FIGURE 1.8 – Central pressure - hPa</figcaption></figure>

The *MWS* and the *CP* of the cyclone have two opposite but similar distributions. The minimal central pressure levels recorded are 879 hPa. This level was reached by Faxai (2001), Dianmu (2004) and Chaba (2004).

* **Radius of maximum sustained wind speed (RMAX)**, in km, is the distance between the location of the *CP* and the location of the *MWS*. It captures the size of the TC event.

<figure><img src="image_12.png" /><figcaption>FIGURE 1.9 – Radius of maximum sustained wind speed - km</figcaption></figure>

The *RMAX* obtained in the IBTRACS database is mostly concentrated around 10 to 50 km, with peaks
up to 200 km for some cyclones. The maximal RMAX is 200 km and was reached by Mindulle (2016).

<figure><img src="image_13.png" /><figcaption>FIGURE 1.10 – RMAX of Sandy (2012)</figcaption></figure>

The picture below shows that the location of the maximum winds of Sandy (2012), in red, are reached in New York, while the track crosses close to Philadelphia. The distance between both cities indicates a RMAX around 150 km.

### 1.5.3 Correlation matrix between variables

We visualize the relationship between each pair of variables X and Y by the Pearson correlation, which can be summarized as follows :

$$p_{X,Y} = \frac{\mathbb{E}[(X - \mu_X)(Y - \mu_Y)]}{\sigma_X \sigma_Y}$$

where :

- $\mu_X$ and $\mu_Y$ are the mean of X and Y

- $\sigma_X$ and $\sigma_Y$ are the standard deviation of X and Y
<figure><img src="image_14.png" /><figcaption>FIGURE 1.11 – Main IBTRACS parameters’ correlation matrix</figcaption></figure>

The maximum sustained wind speed *MWS* variable is negatively correlated to the cyclone central pressure *CP*. A strong correlation relation appears between the *RMAX* and *MWS* and also between *RMAX* and *CP*.

Coriolis effect is highlighted by the positive correlation between the change in longitude and the latitude. This shows that the trajectory of cyclones in general is not rectilinear and that the movement towards the west is accompanied by an increase in latitude.

## 1.6 Limits

In this section, focus is on the data quality of the main database used in this study, IBTRACS [20] and the limitations in the approach taken so far.

### 1.6.1 Limitation in the data

#### Completion of missing data

Some variables are rarely present in the trajectory characteristics. The following table shows the missing values rate per variable.
<figure><img src="image_15.png" /><figcaption>FIGURE 1.12 – Missing data for variables in the dataset</figcaption></figure>

67% of *RMAX* data is missing and labelled as « NA ». *MWS* and *CP* are missing for a small part of the dataset (resp. 15% and 2%), the other variables (*BEAR*, *F*) are well represented with 0% of missing data. Missing data for *RMAX*, *MWS* and *CP* is completed by linear interpolation between the two consecutive non missing time steps.

## Wind speed granularity

<figure><img src="image_16.png" /><figcaption>FIGURE 1.13 – *MWS* profile of TC Vicente 2012</figcaption></figure>

The wind profile of TC Vicente shows that the *MWS* has identical values over several timesteps especially when the cyclone is not yet very powerful. However, the IBTRACS data adequately captures the acceleration and deceleration phases of the *MWS*.

All in all, only 129 *MWS* unique values are taken by the 72,000 cyclones timesteps instances in the dataset. For instance, there are only 15 different values taken above 250 km/h.
1.6.2 Limitations in the approach

Agencies judgmental choice

In IBTRACS [20], Data for any given historical TC event is available for different meteorological agencies of different countries from the NWP basin. The choice of the agency has been consistent across all the database, but relying on judgemental criteria. These criteria include the distance of the agency to the sub-basin of the occurring events. It is expected that IBTRACS provides metadata for this choice in the years to come. This help will matter as results can significantly vary across agencies, for a same TC event.

Length of historical data

The historical data extracted from IBTRACS [20] for this study is short : about 40 years. Since the modelling
heavily relies on this historical data, there is a limitation in using this for estimating the intensity of rare and
extreme events. Insurers are specifically interested in modelling events that cause loss with a Return Period (RP)
over 100 years.

On one hand, such extreme losses originate in their location relative to the insurers’ exposure. For this
dimension, 40 years of historical track data permit to cover highly insured regions. On the other hand, 40 years
of historical windspeed data may not be sufficient to capture the extreme windspeed cause by such events, even
with appropriate distribution.

Missing parameters

Building on the previously mentioned limitation, extreme windspeed cannot be captured adequately by
statistical distributions when data is too scarce. Most of recent litterature on extreme TC windspeed modelling
relies on parametric relationships including the Sea Surface Temperature, see [29] [4].

In these references, TC events have a Maximal Potential Intensity (MPI) [15], which is measure of the maximum reachable windspeed based on environment and atmospheric factors. Sea Surface Temperature (SST) is one of these important factors.

SST is available in ERA5 [27] and can be joined to the IBTRACS data. Its inclusion can be considered as a next step for improvement of this study.

**1.7 Key takeaways**

The main output provided by catastrophe model is the YLT (Year Loss Table), which represents the distri-
bution of one year of catastrophe losses over an insured portfolio.

The first step in the creation of a YLT for Tropical cyclone risk is to identify the phases (origin, propagation, termination) and the modelled subperl : wind.

Afterwards, the study explores two data sources serving as a historical basis for the modelling. These two
sources complement each other : IBTRACS (NOAA) is a large historical set of TC observation data, while
INCYC (REASK) stores physical model data for a small number of cyclones. We only consider IBTRACS data
for the time span 1980-2018 and the spatial scope of North West Pacific Basin.

We then examine parameters that will be of interest in **Chapters 2 and 3**. They are of two types : location and intensity. We noticed strong correlations between intensity parameters such as Maximum sustained wind speed (MWS), Central pressure (CP), and Radius of maximum sustained wind speed (RMAX).

Finally, we take a step back and focus on the limitations in the data and in our approach. They include the lack of granularity for *MWS*, the low availability of *RMAX* in the IBTRACS database, and the absence of Sea Surface Temperature in the modelled parameters. They also include the modelling assumptions we made : Only 39 years of historical data are considered, and our choice of agencies in IBTRACS dataset is arbitrary. While
these limitations are notable, they are not detrimental to the construction of the next elements : the stochastic tracks catalog in **Chapter 2**, the footprints in **Chapter 3** and the YLT in **Chapter 4**.
## 2. Track Generation

After the presentation of location and intensity parameters in **Chapter 1, Chapter 2** focuses on the simulation of these parameters to form a set of data called « stochastic track catalogue ». This set represents the possible trajectories of a large number of simulated TC events. This is the first step towards the building of the YLT.

### 2.1 TC events frequency

Modeling the yearly number of TC events is the first step in building a stochastic track catalog. While some models rely on probability laws (e.g Poisson distribution) to model the frequency of TCs occurrence, we use a different approach to include the year to year correlation of the number of events. This choice is based on inter-annual climate phenomena such as the ENSO (El Nino - Northern oscillation) detailed in the Appendix H.

#### 2.1.1 Block bootstrap methodology

The methodology retained to describe the frequency of TC events is the block bootstrap. The block bootstrap consists in randomly sampling in IBTRACS, with replacement, various blocks of a given size *n* from historical occurrences of TC, and use them to define two main parameters :

- The number of TC events per year

- The Saffir Simpson (SS) category for each of these events.

For instance, if *n*=3, three consecutive historical years are drawn and will determine the number of events for each Saffir Simpson category, for three simulation years. One example is shown below :

Historical years drawn :

- 1980 : three TC events of which one CAT 2, two CAT 4

- 1981 : two TC events of which one CAT 3, one CAT 5

- 1982 : four TC events of which one CAT 2, two CAT 4, one CAT 5

Simulated frequency for a block of three simulations :

- **Simulation 1** : Event 1 - CAT 2; Event 2 - CAT 4; Event 3 - CAT 4

- **Simulation 2** : Event 1 - CAT 5; Event 2 - CAT 3

- **Simulation 3** : Event 1 - CAT 4; Event 2 - CAT 5; Event 3 - CAT 4; Event 4 - CAT 2

This process is repeated until reaching a number of 10,000 simulations. These simulations eventually represent 10,000 simulations of one year of TC occurrence.

The frequency of TCs can be represented by statistical laws (Poisson or Binomial distribution) such as in [4] but the block bootstrap approach has two distinct advantages :

1. The intra-year correlation between events is captured and replicates the historical one, if any.

2. The correlation between simulation years is captured and replicates the historical one, if any.
### 2.1.2 Time dependence structure and size of the bootstrap blocks

In order to assess the dependence in the observation data and to set the size n of the blocks, we use a sample autocorrelation functions of the data (using the acf R package) :

Let $x_1..., x_n$ be observations of a time series, for $h \in [-n; n]$, the sample auto-covariance function is :

$$ \hat{\gamma}(h) := \frac{1}{n} \sum_{t=1}^{n-|h|} (x_{t+|h|} - \bar{x})(x_t - \bar{x}) $$

where $\bar{x}$ is the sample mean of the observations of the time series.

Then, for $h \in [-n; n]$, the sample autocorrelation function is :

$$ \hat{\rho}(h) := \frac{\hat{\gamma}(h)}{\hat{\gamma}(0)} $$

This framework can be applied where each $x_i$ represents the number of TCs in year for a given Saffir Simpson category. Therefore, for each $CAT_i$ with $i \in [1; 5]$, $X_{1980}^{CAT_i}$, ..., $X_{2019}^{CAT_i}$ are the observations of the number of TCs at a given year.

For each $CAT_i$ with $i \in [1; 5]$, the autocorrelation function sample plot below show the value of $\hat{\rho}_{CAT_i}(h)$ for $h \in [1; 15]$.

The horizontal lines on the graph are the bounds $\pm \frac{1.96}{\sqrt{n}}$. In case of an independent year structure (i.e white noise time series), approximately 95% of the sample auto-correlations should fall between the bounds (since 1.96 is the .975 quantile of the standard normal distribution).

Category 1

<figure><img src="image_17.png" /></figure>

Category 2

<figure><img src="image_18.png" /></figure>

Category 3

<figure><img src="image_19.png" /></figure>

Category 4

<figure><img src="image_20.png" /></figure>
Category 5

<figure><img src="image_21.png" /><figcaption>FIGURE 2.1 – Auto-correlation plot of TC frequency depending on Saffir Simpson category</figcaption></figure>

We would expect roughly 15 × (1 − 0.95) = 0.75 values to fall outside the bonds for each plot. Except for CAT 5 TCs, at least 1 value (3 for CAT 2) fall outside the bonds, which means the number of TCs per category between years seems dependent.

Moreover, no value exceeds the bonds for a lag superior to 5. It implies that no time-dependent structure
superior to five years is highlighted by the plots.

Based on these two arguments, this confirms the use of the bootstrap block approach to take into account
the dependence of the TCs between years. In addition, the size n of the bootstrap blocks used in section 2.1 is
set at 5 years.

**2.1.3 Frequency validation**

The validation section focuses on the frequency of TCs by Saffir Simpson categories catalogue of tracks
at two levels : at basin level and at gate level. Validation at these two levels allows for a comprehensive and
fine-grained understanding of the tracks.

In this section and in the following validation selection, the validation of the stochastic catalogue of tracks
will be based on three sets of data :

- A subsample of the catalogue of stochastic tracks generated in this chapter representing 40 years of simulation (referred as source STO).

- The open source catalogue of stochastic tracks STORM published in [4]. This database is described in the Appendix H and is used for information mainly (referred as source STORM)

- The historical tracks from IBTRACS, (referred as source HISTO)

The criteria of success / or not, is the consistency between the historical data (HISTO) and the catalogue
of stochastic tracks (STO).
Frequency of TCs by Saffir-Simpson categories at basin level

The simulated annual frequency of TC (in green) of the stochastic catalogue is compared to historical annual
frequencies from sources HISTO (in red) and STORM (in blue).

<figure><img src="image_22.png" /><figcaption>FIGURE 2.2 – TCs frequency by Saffir-Simpson categories</figcaption></figure>

The bootstrap methodology for TC frequency simulation ensures a good consistency between the annual
frequencies by categories between sources HISTO and STO.

The overall frequency difference between the two data sources does not exceed 0.5 « event per year », which is quite small considering the differences with STORM. However, STORM does not include events with *MWS* below 85 km/h which explains the difference of the first set of bars.

**2.2 TC Track generation : generation of location parameters**

The TC track generation follows the life cycle of a TC described in **Chapter 1**. It relies on various statistical methods. Some of them were already implemented internally for other geographical scopes. In that case, this study helped to show their applicability to the NWP basin. On the contrary, some were newly developed for this basin and applied to other geographies. For that case, the present study helped the company to upgrade its internal modelling standards.

**2.2.1 Initiation**

Three main parameters are generated before starting the track generation :

1. $(lon_0, lat_0)$ is the location of the center of the TC event at t=0.

2. $F_0$ and $dF_0$ are respectively the forward speed and the change in the forward speed of the TC event at $t=0$.
3. $BEAR_0$ and $dBEAR_0$ are respectively the bearing and the change in the bearing of the TC event at $t=0$.

### Generation of starting points ($lon_0, lat_0$)

<table><caption>FIGURE 2.3 - Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Kernel density</td><td>IBTrACS - Western Pacific Basin</td><td>ks, MASS</td></tr></tbody></table>

Kernel is a powerful method to estimate the probability density function of a random variable. We prefer this approach to a random draw within the empirical density function approach. The former brings indeed more variation to the simulated coordinates compared to the historical ones.

We define $f$ a probability density over $\mathbb{R}^d$ with corresponding distribution $F$ with $f$ bounded (i.e $||f||_\infty < \infty$). Let $X = (X_1, ..., X_n)$ be $n$ i.i.d samples drawn from it.

We also define a kernel function $K: \mathbb{R}^d \to \mathbb{R}^+$ such that:

$$ \int_{\mathbb{R}^d} K(u)du = 1 $$

We define, $\hat{f}_n$ the Kernel density estimation of $f$ such that:

$$ \hat{f}_n(x; H) = n^{-1} \sum_{i=1}^{n} K_H(x - X_i) $$

where $H$ is a positive definite and symmetric matrix of dimension $d \times d$.

An interesting property of the function $\hat{f}_n(x; H)$ is its convergence speed to $f$:

$$ |\hat{f}_n(x; H) - f(x)| \le \epsilon_x + \sqrt{\log(n/(nh^d))} $$

where $h$ is the scalar bandwidth of $H$ and $\epsilon_x$ the function of the kernel and the smoothness of $f$ at $x$.

In the next subsection, this approach is used to generate the starting points ($lon_0, lat_0$) in two steps.

**Step 1: Fitting** The density of starting points is estimated through a non parametric estimator of the form:

$$ \hat{f}_n(x; H) = n^{-1} \sum_{i=1}^{n} K_H(x - X_i) $$

where:

- $n$ equals the number of starting points ($lon_0, lat_0$) in IBTRACS.

- $X_i = (X_{i1}; X_{i2})$ with $i \in (1: n)$ are the historical couples of starting points from IBTRACS.

- $H$ is the bandwidth matrix

- $K_H(x) = |H|^{-1/2}K(H^{-1/2}x)$ [12] with $K(x) = (2\pi)^{-1}\exp(-\frac{1}{2}x^T x)$ the Gaussian Kernel.

The bandwidth of the kernel changes its shape. A lower bandwidth limits the scope of the function and leads to the estimate curve looking rough and jagged. Thus, as by [11] the choice of the bandwidth matrix $H$ is crucial
in determining the performance of $\hat{f}_n$. For sake of simplicity, the chosen bandwidth matrix here corresponds to the following Pearson correlation matrix (positive definite and symmetric) :

$$H = \begin{bmatrix} \text{var}(X_1) & \text{cov}_{X_1,X_2} \\ \text{cov}_{X_1,X_2} & \text{var}(X_2) \end{bmatrix}$$

where $X_1$ is the random variable generating the $(X_{i1})$ points, and $X_2$ generating the $(X_{i2})$ points.

**Step 2 : Sampling** ($lon_0, lat_0$) are then sampled from the estimated density above.

This sampling requires :

- Estimating the correlation matrix $H$ (dimension 2x2)

- Sampling a couple of historical starting points ($lon_x$; $lat_x$) from IBTRACS

- Simulating $K_H$, i.e a multivariate normal law with mean = ($lon_x$; $lat_x$) and covariance matrix = $H$

## Validation of starting points

The simulated starting points (SIMU) of the stochastic track catalogue are compared to the starting points of source HISTO and simulated events from STORM.

<figure><img src="image_23.png" /><figcaption>FIGURE 2.4 - Starting points</figcaption></figure>

This graph shows a good overall consistency of the simulated starting points in comparison with the historical
ones. Detailed comments are given in the Appendix H.

### Generation of initial bearing change ((dBEAR₀) point)

<table><caption>FIGURE 2.5 – Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Empirical density</td><td>IBTrACS - Western Pacific Basin</td><td>MASS</td></tr></tbody></table>

Unlike the generation of starting points, it is not necessary to adopt a Kernel approach since historicals dBEAR₀ take a limited number of values.

In this subsection, an empirical distribution based approach is used to generate the bearing change in two steps:

**Step 1: Fitting** The cumulative distribution function (CDF) F of bearing change (dBEAR₀ point) is estimated through a non-parametric estimator of the form:

$$ \forall x \in \mathbb{R}, \forall \omega \in \Omega, F_n(x, \omega) = \frac{1}{n} \sum_{i=1}^{n} 1_{(X_i(\omega) \le x)} $$

where:

- $X_i$ are the historical bearing changes of starting points from IBTRACS. They are differences between the first two bearings of each historical events from IBTRACS, denoted by $BEAR_0$ and $BEAR_1$.

- n is the number of historical $BEAR_0$ and $BEAR_1$ taken in the IBTRACS dataset

**Step 2: Sampling** dBEAR₀ is then sampled from the empirical CDF above.

### Generation of the remaining parameters

Initial bearing ($BEAR_0$), Forward speed ($F_0$) and Forward speed change ($dF_0$) are sampled following the same methodology as $dBEAR_0$.

### 2.2.2 Track propagation

#### Transition probabilities fitting

<table><caption>FIGURE 2.6 – Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Hidden Markov Chain</td><td>IBTrACS - Western Pacific Basin</td><td>n/a</td></tr></tbody></table>

The track propagation model uses the Markov Chain methodology and follows [18]. The main advantage of this model is the smoothness of the simulated tracks. This is due to the chosen states (dF and dBEAR), rather than the coordinates directly.
**Markov Chain framework** A discrete-time Markov Chain is a sequence of discrete random variables $X_1, X_2, X_3, ...$ with the Markov property, namely that the probability of moving to the next state depends only on the present state and not on the previous states.

Let $(X_t), t \in \mathbb{N}$ be a sequence of discrete random variables with value in $E$. $(X_t)_{t \in \mathbb{N}}$ is a discrete Markov Chain of matrix of transition $P$ of general term $(p_{i,j})$ $(i,j) \in E^2$ if $\forall (x_0, ..., x_t) \in E^{t+1}$:

$$ \mathbb{P}(X_t = x_t | X_1 = x_1, X_2 = x_2, ..., X_{t-1} = x_{t-1}) = \mathbb{P}(X_t = x_t | X_{t-1} = x_{t-1}) = p_{x_{t-1}, x_t} $$

Moreover, if the transition mechanism remains the same over time, the Markov Chain is homogeneous and we have $\forall t \in \mathbb{N} \forall (x_0, ..., x_t) \in E^{t+1}$:

$$ \mathbb{P}(X_t = x_t | X_1 = x_1, X_2 = x_2, ..., X_{t-1} = x_{t-1}) = \mathbb{P}(X_1 = x_1 | X_{t-1} = x_{t-1}) = p_{x_{t-1}, x_t} $$

In that case, we obtain two properties of a homogeneous Markov Chain:

1. $P$ is a stochastic matrix i.e

$$ (1) - \forall i \in E, \sum_{j \in E} p_{i,j} = 1 $$

2. The law of a Markov Chain $(X_n)_{n \in \mathbb{N}}$ is entirely determined by its transition matrix $P$ and the law of $X_0$ i.e $\forall t \in \mathbb{N} \forall (x_0, ..., x_t) \in E^{t+1}$:

$$ (2) - \mathbb{P}(X_1 = x_1, X_2 = x_2, ..., X_{t-1} = x_{t-1}, X_t = x_t) = P[X_0 = 0] \times p_{x_0, x_1} \times \dots \times p_{x_{t-1}, x_t} $$

**Application** In our modeling, we will assume that the evolution of $dBEAR$ (resp. $dF$) is a discrete and homogeneous Markov Chain denoted $(dBEAR_t)_{t \in \mathbb{N}}$ (resp. $(dF_t)_{t \in \mathbb{N}}$). The choice of such a model is based on the fact that the future trajectory of a TC is based solely on its present location. Moreover, we discretize the variable $dBEAR$ (resp. $dF$) which is continuous by rounding its values.

According to property (2) above, $(dBEAR_t)_{t \in \mathbb{N}}$ (resp. $(dF_t)_{t \in \mathbb{N}}$) is defined by its initial law $dBEAR_0$ (resp. $dBEAR_0$) and by the transition matrix $P_{dBEAR}$ (resp. $P_{dF}$).

The initial law $dBEAR_0$ (resp. $dF_0$) is already defined in the initiation section.

**Transition Matrix** We determine the transition matrix $P_{dBEAR}$ as follows:

Given $E$, a space composed of discretized values in $[-180; 180]$, $\forall (i, j) \in E^2$, $\forall t \in \mathbb{N}$:

$$ P_{dBEAR_{i,j}} = P(dBEAR_{t+1} = j | dBEAR_t = i) = \frac{\text{Card} (\{x = j, x \in dBEAR_{dBEAR=i}^{J+1}\})}{\text{Card} (\{x = i, x \in dBEAR^J\})} $$

where:

- $dBEAR^J$ represents the set of $dBEAR$ values in the historical dataset.

- $dBEAR_{dBEAR=i}^{J+1}$ represents the set of $dBEAR$ values given that previous $dBEAR$ was equal to $i$ in the historical dataset.

However, as seen in chapter 1, extreme values of $dBEAR$ (i.e. $> 35$ or $< -35$) are scarce and $dBEAR$ does not take value over all the $[-180; 180]$ interval. To limit the number of null values in the transition matrix, we
judgmentally group the dBEAR into six bins ($S_k$)$_{k \in [1; 6]}$. Then, $\forall x \in E \ \forall t \in \mathbb{N}$:

$$P_{dBEAR_{i,j}} = P(dBEAR_{t+1} = j | dBEAR_t = i) = \sum_{k=1}^{6} P(dBEAR_{t+1} = j | dBEAR_t \in S_k) \times \mathbb{1}_{dBEAR_t \in S_k}$$

with:

$$P(dBEAR_{t+1} = j | dBEAR_t \in S_k) = \frac{Card(\{x = j, x \in dBEAR_{dBEAR \in S_k}^{J+1}\})}{Card(\{x \in S_k, x \in dBEAR^J\})}$$

where:

- $dBEAR_{dBEAR \in S_k}^{J+1}$ represents the set of dBEAR values given that previous dBEAR was in the bin $S_k$ in the historical dataset.

- $S_1 = [-180; -35[, S_2 = [-35; -5[, S_3 = [-5; -0[, S_4 = [0; 5[, S_5 = [5; 35[, S_6 = [35; 180]$

**Transition Matrix** $dF$ $P_{dF}$ is defined in a similar way than $P_{dBEAR}$ with intervals ($F_k$)$_{k \in [1; 5]}$ such as:

$$F_1 = [-100; -5[, F_2 = [-5; 0[, F_3 = [0; 5[, F_4 = [5; 15[, F_5 = [15; 100]$$

## Propagation of $BEAR_t$ and $F_t$

For $t \ge 1$, $BEAR_t$ and $F_t$ are given by:

$$\left\{ \begin{array}{l} BEAR_t = BEAR_{t-1} + dBEAR_t dt \\ F_t = F_{t-1} + dF_t dt \end{array} \right.$$

where:

- $BEAR_0$ is sampled following section 2.2.1 and $BEAR_{t-1}$, $t \ge 2$ is calculated at the previous recurrence step.

- $dBEAR_t$ is sampled from the conditional distributions defined previously.

- $F_0$ is sampled following section 2.2.1 and $F_{t-1}$, $t \ge 2$ is calculated at the previous recurrence step.

- $dF_t$ is sampled from the conditional distributions defined previously.

- $dt = 6$ (in hours)

## Propagation of $lat_t$ and $lon_t$

For $t \ge 0$, $lat_t$ and $lon_t$ are given by, according to [8]:

$$\left\{ \begin{array}{l} lat_{t+1} = lat_t + \frac{180}{\pi R} \times F_t \times dt \times \sin \frac{BEAR_t \times \pi}{180} \\ lon_{t+1} = lon_t + \frac{F_t \times 180 \times dt \times \cos \frac{BEAR_t \times \pi}{180}}{\cos \frac{lat_t \times R \times \pi}{180}} \end{array} \right.$$

where:
- $lat_0$ and $lon_0$ are sampled following section 2.2.1

- $lat_t$ and $lon_t$ are calculated at the previous recurrence step for $t \ge 1$

- $F_t$ and $BEAR_t$ are calculated at the previous step

- R is the earth radius (in km)

- $dt = 6$ (in hours)

### 2.2.3 Termination

The track should be terminated to represent TC's dissipation. To model track termination, a survival probability based on the historical events of IBTRACS database[20] is computed at each timestep,

**Fitting** For a modelled TC, given its last timestep coordinates ($lon$ and $lat$) the probability of survival is calculated as follows at each new timestep:

$$P[TC_{survival} | (lon, lat)] = 1 - \frac{TC_{100kmradius | (lon, lat)}^{end}}{TC_{100kmradius | (lon, lat)}}$$

where:

- $TC_{survival}$ is the TC's survival event.

- $TC_{100kmradius | (lon, lat)}^{end}$ accounts for the number of historical TCs at less than 100 km from last step coordinates($lon$ and $lat$) which ended at the following timestep.

- $TC_{100kmradius | (lon, lat)}$ accounts for the number of historical TCs in the 100 km - radius from last step coordinates ($lon$ and $lat$).

100 km is a good trade-off in order to obtain sufficient number of points without considering too large and heterogeneous (land vs. sea) geographical areas.

**Sampling** At each timestep, starting from the tenth timestep, we compute $TC_{survival}$ and generate a random value $s$ from a uniform law $U[0; 1]$.

If $s > TC_{survival}$: the TC goes on to the next propagation timestep. Otherwise the TC stops.

**Comments** The process starts from the tenth timestep to match historical TCs duration. Insights to improve this approach are rare in the literature. Some publications like [8] base the track termination on the MWS value (if MWS goes below 17 km/h then the track stops).

### 2.2.4 Validation

Track validation is performed at NWP basin scale and at « gate » scale. This allows us to assess the consistency of track behavior: i.e. over- or under-representation of tracks in one area compared to the historical record, at two different scales.

Particular attention is paid to the landfall variable since TCs cause damages mostly on lands.
## Basin level validation

### Frequency of landfalling TCs by Saffir Simpson categories, at basin level

<figure><img src="image_24.png" /><figcaption>FIGURE 2.7 – Frequency of landfalling hurricanes by categories</figcaption></figure>

The simulated (STO) and historical (HISTO) annual frequencies are similar, while STORM shows a truncated frequency for tropical storms (TS), and lower annual frequencies for CAT 5 events.

Another validation graph which gives a focus on the landfall frequency at a country level is given in the appendix H.

## Country Gate level validation

**Frequency of TCs at gate level in Philippines** This validation step focuses on the Philippines, which is a country with high exposure to TCs.

We compare the simulated annual frequencies from source STO and source HISTO at a gate level where gates are a sequence of 100 km segments that board the east coastline.

They are shown here:
<figure><img src="image_25.png" /><figcaption>FIGURE 2.8 – Gates plot Philippines</figcaption></figure>

The size of 100 km is judgmentally chosen.

In order to compare the historical gates frequencies to the simulated ones, we choose a bootstrapping approach to build subsets of simulated data.

More precisely, 100 blocks of 39 years of simulations are randomly drawn into the stochastic catalogue of tracks, with replacement, a certain number of time. For all these blocks, the number of TCs which intersects each gate every year is calculated. From this count, we obtain annual frequency from each gate and block. For each gate, a distribution is fitted on the blocks annual frequencies.

Consequently, we build a ridge plot which consists in showing for each gate :

- The historical annual frequency of TC events crossing the gate, represented by a black line.

- The density of the simulated annual frequencies fitted distribution from blocks

<figure><img src="image_26.png" /><figcaption>FIGURE 2.9 – Gate level validation Philippines</figcaption></figure>
There is a slight trend of having historical black line on the right-part of the distributions in the south of the country and the opposite effect in the north. Considering the historical black line as a draw in the distribution, no gate shows a line in the tails of the distribution beyond the 99th quantile.

Thus, the graph shows an overall good consistency between STO and HISTO data sources. Two other validation graphs for Japan and China are shown in the Appendix H.

## 2.3 TC Track generation : generation of intensity parameters

The previous section showed the generation of a sequence of $(lon_t, lat_t)$ that represents the position of TC event center over time.

This section shows how to assign four intensity parameters to these sequence :

- The maximum sustained wind speed reached at time t, $MWS_t$

- The radius of maximum speed at time t, $RMAX_t$

- The angle of maximum winds at time t, $AMAX_t$

- The minimal pressure reached at time t, $CP_t$

The models used to generate these parameters are machine Learning and deep Learning methods given the number of potential predictive values and possible relationships between explanatory variables and target values. The use of interpretability methods, especially when manipulating random forest, allows to evaluate the functioning of the models to get out of the black box issue.

### 2.3.1 Maximum sustained wind speed : MWS

<table><caption>FIGURE 2.10 - Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Deep learning - LSTM</td><td>IBTrACS - Western Pacific Basin</td><td>Keras, Tensorflow</td></tr></tbody></table>

There are two main challenges in estimating this parameter :

1. The maximum sustained wind speed is a key intensity parameter for building a TC model, as it directly relates to economic damage.

2. The parameter to estimate is a vector, $\forall t, MWS_t$, because we want to take into account the correlation between each timesteps.

This parameter is determined using a deep learning approach, and more precisely, LSTM [32]. Other methods such as K-Nearest Neighbours regression, Gated Recurrent Unit, Recurrent Neural Network and regression
Forest, that are shown in the validation part, are described in the Appendix H.

The problem that is solved through this methodology is the estimation of a conditional law $Y|X$ where:

* $Y$ is the target vector of maximum sustained wind speed $MWS$ of dimension $T$.

$$Y = (MWS_1, ..., MWS_T)$$

where:

- $MWS_t$ is the maximum sustained wind at time $t$.

- $T$ is the last timestep of the simulated TC.

* $X$ is the predictor matrix of dimensions $T \times p$ with $p = 7$, the number of predictors.

$$X = \begin{bmatrix} \text{lat}_0 & \text{lon}_0 & F_0 & \text{BEAR}_0 & dF_0 & \text{dBEAR}_0 & \text{dist2coast}_0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ \text{lat}_T & \text{lon}_T & F_T & \text{BEAR}_T & dF_T & \text{dBEAR}_T & \text{dist2coast}_T \end{bmatrix}$$

These predictors are chosen based on the fact that there are generated before $MWS$ in the model simulation process.

## LSTM: Long short-term memory

The underlying principles and application steps of the LSTM algorithm are described thereafter:

**Main principles** LSTM is an extension of neural networks which are designed to model chronological sequences and their long-range dependencies. Similarly to a neural network, it operates by updating the weights of each network cell through back-propagation mechanism led by a gradient descent algorithm. Structure of a LSTM and algorithm mechanism are shown in Appendix H.

**Rationale of use** Due to its time-dependency mechanism, where the value at time $t-1$ impacts the value at time $t$, the LSTM is the best candidate to model the $MWS$ evolution through times of TC.

**Main application steps** Several techniques are used to optimize the LSTM algorithm.

**Loss function** The performance and fitting of the algorithm is measured by the MSE (mean squared error) given by:

$$MSE = \frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{Y}_i)^2$$

where:

- $Y$ is a variable representing $MWS$ in the validation dataset

- $n$ the number of observations in the validation dataset

- $\hat{Y}$ is the predicted $MWS$

## Model structure
<figure><img src="image_27.png" /><figcaption>FIGURE 2.11 Neural network structure with LSTM layer</figcaption></figure>

The graph summarizes the general structure of the network. It consists of an input layer, a masking layer, a (bidirectional) LSTM layer, a Dense time distributed layer and the output layer based on [32] (Appendix H).

We will then compare the choice of the model hyper parameters by finding the values which minimize the MSE over all the epochs (i.e. the number of times that the learning base is entirely taken as an input by the model). Then we will choose the number of epochs that minimizes the overall model (Appendix H).

## Optimizer - Gradient descent algorithm

We compare the *MSE* for different gradient descent algorithms of the following list *RMSprop*, *AdaGrad*, *Adam* on the LSTM model ran over several number of epochs on the table below :

<figure><img src="image_28.png" /><figcaption>FIGURE 2.12 - Gradient descent algorithm - performance comparison</figcaption></figure>

Based on the graph, the chosen optimizer for the LSTM is **Adam**.

## Learning rate

Due to lack of computational power, we have restricted ourselves to three values representing distinct orders of magnitude of learning rate. 0.1 corresponds to a prevalence of what the algorithm has learned recently, 0.01 corresponds to a balance between recent and old learning and 0.001 gives a premium to the first learning.

We compare the *MSE* of the three learning rates on the LSTM ran over several number of epochs on the graph below :
<figure><img src="image_29.png" /><figcaption>FIGURE 2.13 – Learning rate - performance comparison</figcaption></figure>

Based on the graph, the learning rate chosen for the LSTM is 0.01.

## Number of epochs

We compare the *MSE* of the LSTM for different number of epochs on the learning curve below :

<figure><img src="image_30.png" /><figcaption>FIGURE 2.14 – Number of epochs - learning curve</figcaption></figure>

The MSE gradually decreases on the training data, but we do not rely on it to avoid over fitting.

Regarding the validation data, the MSE of the LSTM model decreases very quickly up to the twentieth epoch. Then the mean error seems to stabilize with less and less MSE volatility between each epoch. The number of epochs chosen for the model is 125. This choice is based on the fact that the MSE of the validation database stabilizes at the lowest level.

## Validation of MWS

We show the consistency between the historical *MWS* and the simulated *MWS* for historical events.

This criteria checks whether the generation methodology [section 2.4] can reproduce the MWS of the following historical events :

- CAT 4 TC MIREILLE (1991), landfalling in Japan.

- CAT 5 TC HAIYAN (2013), landfalling in Philippines.

- CAT 3 TC HATO (2017), landfalling in China.
The graphs below show the modelled and historical MWS of these events, based on different methodologies

<figure><img src="image_31.png" /><figcaption>FIGURE 2.15 – Maximum sustained winds evolution for three major events : MIREILLE, HATO and HAIYAN</figcaption></figure>

These graphs show that the Gated Recurrent Unit method is the best methodology to capture the *MWS* variability of the historical events. In particular, it manages to reproduce the acceleration and deceleration phases of the *MWS* of historical cyclones.

Overall, the relative evolution of the MWS is well reproduced even if the absolute values of MWS are slightly underestimated for all the three events, and at all timesteps. This is an issue since this intensity variable is the cornerstone for estimating the level of destruction of affected areas and thus economic losses.

<table><caption>FIGURE 2.16 – Metrics comparison for the different maximum sustained wind prediction models</caption><thead><tr><th></th><th></th><th>KNN Regression</th><th>RNN</th><th>GRU</th><th>LSTM</th></tr></thead><tbody><tr><td>MSE</td><td></td><td>38.90</td><td>38.96</td><td>33.69</td><td>32.57</td></tr><tr><td>MSE by TC</td><td></td><td>30.32</td><td>31.31</td><td>27.33</td><td>26.62</td></tr><tr><td rowspan="5">MSE by TC Saffir Simpson Categories</td><td>CAT 1</td><td>30.80</td><td>32.61</td><td>28.12</td><td>26.62</td></tr><tr><td>CAT 2</td><td>35.88</td><td>36.06</td><td>33.79</td><td>29.85</td></tr><tr><td>CAT 3</td><td>46.49</td><td>40.08</td><td>35.35</td><td>36.76</td></tr><tr><td>CAT 4</td><td>49.72</td><td>55.59</td><td>45.07</td><td>35.92</td></tr><tr><td>CAT 5</td><td>59.76</td><td>56.86</td><td>50.89</td><td>50.88</td></tr></tbody></table>

Looking at the general evaluation measures of the models used, LSTM has the most convincing error results in general but also in particular for the CAT 3,4 and 5 TC.
### 2.3.2 Central pressure : CP

<table><caption>FIGURE 2.17 – Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Power law regression</td><td>IBTrACS - Western Pacific Basin</td><td>n/a</td></tr></tbody></table>

The central pressure parameter at time t, denoted by $CP_t$, is estimated through a parametric relationship from [36] given by :

$$CP_t = CP_{env} - \exp\left(\log\left(\frac{MWS_t}{6.7} * \frac{1}{0.644}\right)\right)$$

where :

* 6.7 and 0.644 are regression parameters from [36]

* $CP_{env}$ is the mean sea level pressure; here a constant value of 1013 hPA is taken as it is the mean sea-level atmospheric pressure on Earth.

Most Pressure Wind relationships are of the form $V_m = \alpha \Delta p^\beta$ ([31]) where $\Delta p^\beta$ is the pressure deficit (in hPa) between the cyclone center and the environmental pressure. The relationship used here is also used in [2] or in the Dvorak pressure-wind model ([13]).

### 2.3.3 Radius of maximum winds RMAX

<table><caption>FIGURE 2.18 – Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Machine learning - Regression Forest</td><td>IBTrACS - North America Basin</td><td>QRF</td></tr></tbody></table>

As seen in chapter 1, RMAX variable is not always available in IBTRACS. To overcome this problem, the learning base for the RMAX study is based on IBTRACS data from the year 2000. This parameter is predicted using a regression forest approach [5] as it offers a good balance between interpretability of the predictors and performance on continuous variable prediction.

The problem to solve here is the estimation of a conditional law $Y|X$.

This estimation to the conditional law is represented by the regression function $m(x) = E[Y|X = x]$, given a training sample $D_n = ((X_1, Y_1), ..., (X_n, Y_n))$, where :

* $Y_t$ is the value of RMAX at time t

$$Y_t = RMAX_t$$

* $X_t$ is the predictor vector of dimension $p = 7$.

$$X_t = (lon_t, lat_t, dlon_t, dlat_t, MWS_t, CP_t, BEAR_t, dMWS_t)$$

These predictors are chosen according to the process described afterwards.
## Regression Forest

**Main principles** A regression forest (resp. random forest) [5] [28] is a supervised learning method which uses ensemble learning method such as bagging and bootstrapping for the estimation of conditional probability laws. It operates by constructing multiple decision trees at training and outputting the mean prediction (resp. majority prediction) over the individual trees.

<figure><img src="image_32.png" /><figcaption>FIGURE 2.19 – Random Forest</figcaption></figure>

The regression forest builds an estimate $m_n : \mathcal{X} \to \mathbb{R}$ of the regression function $m$.

A regression forest is a predictor consisting of a collection of $M$ randomized regression trees. For the $j$-th tree in the family, the predicted value at the query point $x$ is denoted by $m_n(x; \Theta_j, D_n)$ where $\Theta_1, ..., \Theta_M$ are independent random variables, distributed the same as a generic random variable $\Theta$ and independent of $D_n$

In practice, the variable $\Theta$ is used to resample the training set prior to the growing of individual trees and to select the successive directions for splitting—more precise definitions will be given later. In mathematical terms, the $j$-th tree estimate takes the form

$$m_n(x; \Theta_j, D_n) = \sum_{i \in D_n^*(\Theta_j)} \frac{1_{X_i \in A_n(x; \Theta_j, D_n)} Y_i}{N_n(x; \Theta_j, D_n)} = \sum_{i \in D_n^*(\Theta_j)} w_i(x, \Theta_j) Y_i$$

where $D_n^*(\Theta_j)$ is the set of data points selected prior to the tree construction, $A_n(x; \Theta_j, D_n)$ is the cell containing $x$, and $N_n(x; \Theta_j, D_n)$ is the number of points that fall into $A_n(x; \Theta_j, D_n)$.

Then, threes are combined to form the forest estimate :

$$m_{M,n}(x; \Theta_1, ..., \Theta_M, D_n) = \frac{1}{M} \sum_{j=1}^{M} m_n(x; \Theta_j, D_n)$$

which can also be expressed as a weighted sum over all observations :

$$m_{M,n}(x; \Theta_1, ..., \Theta_M, D_n) = \sum_{j=1}^{M} w_j(x) Y_j$$

where $w_j = \sum_{i \in D_n^*(\Theta_j)} w_i(x, \Theta_j)$.

**Convergence of a regression forest** As the number $n$ of observation grows, it does not lead to overfitting as shown in [5] :

$$\lim_{n \to \infty} \mathbb{E}[m_{M,n}(X; m_{M,n}(x; \Theta_1, ..., \Theta_M) - m(X))^2] = \mathbb{E}[m_{\infty,n}(X) - m(X)]^2$$
**Rationale of use** Even if the regression forest does not take into account the temporal dependence between the different *RMAX* values of a TC, it is still useful because it offers a granular level of interpretation: we can capture the importance of each of the explanatory variables for the prediction. Moreover, the properties of the regression forests ensure that overfitting is avoided as the number of trees increases [5].

We will use the interpretability of the algorithm to explain how we reduced the explanatory variables to the *X* vector above by doing feature selection.

**Step 1 : selection of features** Feature selection is done through the Accumulated Local Effect (ALE) method [44]. The selection of the features is done according to the non flatness of the feature to predictor curve. The ALE plot is a more efficient alternative to Partial Dependence Plot and M-plot. The three techniques are described in the Appendix H.

The way the accumulated local effect plot proceeds is shown using the example of the house price. The effect on the price associated with an area of 100*m*² is analysed. The ALE technique uses all houses with about 100*m*², gets the model predictions pretending those houses were 102*m*² minus the prediction pretending they were 98*m*². This gives us the pure feature effect of the area variable, without taking into account the effect of the correlated features. In summary it shows how the model predictions change in a small "window" of the feature around the target variable for data instances in that window.

The graph below details the ALE plot of each variable studied in a preliminary way so as to integrate the best variables as predictors of the *RMAX*.

<figure><img src="image_33.png" /><figcaption>FIGURE 2.20 – Accumulated local effect plot for each candidate variable to be a component of *X*</figcaption></figure>

Flat graphs show no correlation between *dBEAR*, *F*, *dF* and the predictor *RMAX*, while irregular graphs indicate a correlation structure.

For instance, the variable *MWS* shows a strong monotonous decrease in the ALE plot and hence appears as strongly anti correlated to *RMAX*. On the contrary, the graph of the variable *CP* shows a major positive
correlation with RMAX

Based on this plot, *lon*, *lat*, *dlon*, *dlat*, *MWS*, *CP*, *BEAR* and *dMWS* (change of Maximum sustained wind speeds) are selected to predict RMAX. Another method (Boruta Algorithm) has been tested and is shown in the Appendix H.

**Step 2: Trees fitting and prediction** The performance and fitting of the algorithm is measured by the MSE (mean squared error).

Parameters are fitted to minimize the MSE and their final values are as follows :

— Number of trees : 97 as it the value that minimises the mean square error

<figure><img src="image_34.png" /><figcaption>FIGURE 2.21 – RMAX - Mean square error by number of trees</figcaption></figure>

— Number of tested iterations at each division : 8. It is based on a search grid tuning of hyper-parameters [Appendix H].

## Validation of RMAX

The final MSE obtained on the test dataset arrives at 211, i.e the root mean square error of the RMAX predicted (in comparison with the validation data) is 14.5.

The breakdown of the MSE by TC categories shows that the MSE is much larger for intense TCs (CAT 3,4 and 5). This is likely due to (i) the scarcity of historical data for such events (ii) the higher variability of RMAX for these events, which the algorithm fails to capture.

<table><caption>FIGURE 2.22 – MSE by Saffir Simpson category</caption><thead><tr><th>Category Saffir Simpson</th><th>MSE</th><th>Mean error</th></tr></thead><tbody><tr><td>0</td><td>111</td><td>11</td></tr><tr><td>1</td><td>178</td><td>13</td></tr><tr><td>2</td><td>225</td><td>15</td></tr><tr><td>3</td><td>314</td><td>18</td></tr><tr><td>4</td><td>647</td><td>25</td></tr><tr><td>5</td><td>1144</td><td>34</td></tr></tbody></table>

We focus on the consistency between the historical RMAX and the simulated RMAX for historical events and the confidence interval around the simulated RMAX.

This criteria checks whether the RMAX prediction methodology (hyperlinkrfsection 2.3.3) can reproduce the RMAX of the following historical events :
- CAT 5 TC Haiyan (2013)

- CAT 5 TC Jebi (2018), landfalling in Japan

- CAT 4 TC Vicente (2012), landfalling in multiple countries

These three TCs are in the validation base.

Additionally, the criterion examines the uncertainty in the RMAX prediction by building 95% - confidence intervals around the regression forests predicted values. Then it checks whether the historical values lie in this interval. Thus, in the case of a consistent model, we should see 95% of the predicted values within this interval. The 95% - confidence intervals around the predicted value are generated using a quantile regression forest ([43]) approach described in Appendix H.

The graphs below show 95% - confidence intervals around the prediction and the historical values of RMAX :

<figure><img src="image_35.png" /><figcaption>FIGURE 2.23 – RMAX Comparison for major events</figcaption></figure>

- The red line represents the prediction of RMAX

- The blue lines represents the low and upper bounds of the 95% - confidence intervals.

- The black line represents the historical values of RMAX

The behavior of the simulated *RMAX* is much more stable than the historical values especially for HAIYAN and JEBI. This is due to the fact that unlike the LSTM model for *MWS*, the regression forest model for *RMAX* does not capture the variability between each timestep since it is not time dependent.

Of these three events, only four *RMAX* values fall outside the confidence interval over 150 timesteps (≈ 2.7%). If we extrapolate this data over the whole validation set, about 2.4% < 5% of the predicted values fall outside the confidence interval.

Despite a less satisfactory evolution profile of *RMAX* than *MWS*, the prediction model shows a good consistency.
### 2.3.4 Angle of maximum winds : AMAX

<table><caption>FIGURE 2.24 – Method, data and R packages used</caption><thead><tr><th>Methodology</th><th>Data used</th><th>R Package</th></tr></thead><tbody><tr><td>Machine learning - Circular Forest</td><td>INCYC</td><td>cirttree</td></tr></tbody></table>

This parameter is predicted using a circular regression forest approach [22], which is well adapted to the angular values of AMAX.

The problem to solve here is the estimation of a conditional law $Y|X$, where:

* $Y_t$ is the value of AMAX at time $t$

$$Y_t = AMAX_t$$

* $X_t$ is the predictor vector of dimension $p = 9$ at time $t$

$$X_t = (lon_t, lat_t, RMAX_t, F_t, BEAR_t, dlon_t, dlat_t, CP_t, MWS_t)$$

The selection of features and the simulation methodology is similar to the RMAX and explained in the Appendix H.

## 2.4 Tracks simulation procedure

### 2.4.1 Procedure

After having described:

- In section [2.1], how to simulate the occurrence of TC events for a given number of simulation years.

- In sections [2.2 to 2.4], for all of these events how to generate track parameters, the following simulation procedure is implemented:

- **Step 1**: draw a block of five consecutive historical years in the historical database from the IBTRACS database as built in chapter 1.

- **Step 2**: Select for each of these years, the number of TC events for each category.

- **Step 3**: Replicate this data to initiate a block of n consecutive simulated years, i.e with the same number of TC events for each category.

- **Step 4**: For each of these events, generate location parameters. [section 2.2]

- **Step 5**: For each of these events, generate intensity parameters. [section 2.3]

- **Step 6**: For each of these events, compute survival probability at each timestep to terminate the track.

Repeat this procedure until simulating 10,000 years. At the end of this procedure, we get 10,000 simulated years, associated to TC events, and for each of these events, a track. This set of data is called a « stochastic
A comparison between a set of simulated tracks and a set of historical tracks is shown below :

<figure><img src="image_36.png" /><figcaption>FIGURE 2.25 – A slot of 10 years of simulated tracks (left) and of historical tracks (right)</figcaption></figure>

- The simulated tracks and the historical tracks have similar form : they are convex, going first northwest and heading then to the northeast as getting close from the coast. The Coriolis effect is then well reproduced by the model.

- When hitting the land, the modelled TC tracks have a dissipation behavior that is similar to the historical tracks.

## 2.4.2 Remarks on the procedure

This procedure generates each track independently. Parallel computation mechanisms can then be applied to run the simulation procedure very efficiently.

These parallel mechanisms can be :

- Open source parallel libraries to run the procedure on a laptop in a reasonable time (R package : parallel)

- A cluster on the cloud to run the procedure in a very short amount of time.

# 2.5 Limits

## 2.5.1 Historical based approach for the track generation

The track generation methodology heavily relies on historical data. When historical data is scarce, for instance for latitudes or categories where too few events are recorded, the simulated locations of the TC center do not compare well with the historical ones. For instance, for high latitudes where no events are recorded the Markov Chain algorithm fails to simulate events there as shown in the north part of Japan in the Appendix H. However, it has been checked that AXA Group has no significant exposure in these regions.

## 2.5.2 Underestimation of the *MWS*

The LSTM methodology slightly underestimates the *MWS* over time, especially for rapid intensification events. This underestimation is observed by reproducing historical events and comparing the predicted values with the historical values of *MWS* over time. To overcome this issue, adjustments have been proposed and can be implemented in the next versions. One of this adjustment is in two steps :

- First, simulate the maximum *MWS* reached over time for the TC event.

- Second, predict the *MWS* for all time steps relatively to the maximum *MWS*.
## 2.6 Key takeaways

**Chapter 2** provided all the elements to build a stochastic catalog of tracks: the frequency of the TCs with the origin, propagation and termination of their tracks, i.e. the location parameters, and the intensity parameters as well.

First, a block bootstrap methodology is used to model the annual frequency of TCs. This method has the advantage of preserving the year-to-year correlation of the number of TC events.

Second, the approach used to model the location parameters through the TC's life cycle relies on several statistical methods. Tracks origin parameters (*F*, *dF*, *BEAR*, *dBEAR*) are modeled using empirical CDFs, while the starting coordinates points (*lon*₀, *lat*₀) are modeled using a Kernel density approach. While the first method ensures consistency with IBTRACS data, the second method provides more variability in the simulated starting points.

Moreover, Markov Chains based on *dF* and *dBEAR* are used to simulate the track propagation following the work of [18], whereas the modeling of the track termination is achieved through fitting empirical probabilities of survival.

While the generation of the location parameters relied on classical statistical methods, the intensity parameters were generated through machine learning methods: LSTM for *MWS*, Random Forest for *RMAX*, and Angular Random Forest for *AMAX*. These methods make it possible to take into account complex relationships and temporal correlations between variables.

Third, all the previous modelling steps helped to build a simulation procedure. This procedure enables the simulation of 10,000 years of TC event activity, with a track associated to each of the events. This output is named « stochastic track catalogue ».

Lastly, the main weakness of this track generation method is the exclusive reliance of the model on historical data, which makes it impossible to simulate tracks where there is no historical data.
# 3. Footprint Generation

After **Chapter 2** shows the building of a « stochastic track catalogue », **Chapter 3** provides the methodology to create a set of wind footprints. These footprints represent wind intensity outside of the points from the stochastic track catalogue, and account for terrain variability effects. This is the second step towards the building of the YLT.

## 3.1 Prediction of windspeed fields : framework and challenges

### 3.1.1 Problem framework

An example of track and windfield is shown below for TC Mangkhut (2018) :

<figure><img src="image_37.png" /><figcaption>FIGURE 3.1 – Visual representation of a TC windfield</figcaption></figure>

On this graph, the track is represented by the points in red. The windfields enable the generation of the inland shades of red, shown in the north of the Philippines and in China. Chapter 3 undertakes to calculate
these windfields, based on the track information provided in Chapter 2.

The problem to solve is the estimation, at each timestep, of a conditional law $Y|X$, where:

* $Y_t$ is the windfield at time $t$, with spatial resolution of 1 km and of dimension $d = 4RMAX \times 4RMAX$.

$$ Y_t = \begin{bmatrix} MWS_{t,1,1} & \cdots & MWS_{t,1,4RMAX} \\ \vdots & \ddots & \vdots \\ MWS_{t,4RMAX,1} & \cdots & MWS_{t,4RMAX,4RMAX} \end{bmatrix} $$

* $X_t$ is the track vector at time $t$, with intensity and locations parameters simulated at Chapter 2.

$$ X_t = (MWS_t, RMAX_t, AMAX_t, CP_t, lon_t, lat_t, F_t, dF_t, BEAR_t, dBEAR_t) $$

### 3.1.2 Challenges

The previously defined estimation problem brings a major challenge:

Prediction values of a matrix of $4RMAX \times 4RMAX$ pixels based on a vector of dimension 10 can lead to overfitting. A PCA approach is then applied to reduce the dimension of the windfield to predict, following [42].

We believe this approach is superior to the parametric radial wind profiles [50], [33] used in most of catastrophe commercial software. The implementation of this approach significantly upgraded the modelling of TC wind hazard at AXA Group.

### 3.1.3 Data used

The database used for the estimation of both $Y$ and $X$, as previously defined, is INCYC. It contains, for a certain number of historical events, tracks and associated windfields, at 1km resolution, at sea and at 10m above sea level. The track parameters provided by INCYC are consistent with the track parameters simulated at Chapter 2.

## 3.2 Dimension reduction

We focus on this section, in the reduction of windspeed fields into three principal components using a Principal component analysis.

<table><caption>FIGURE 3.2 – Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Principal Component Analysis</td><td>INCYC</td><td>stats (prcomp)</td></tr></tbody></table>

The problem to solve is the transformation of the windfield at time $t$, denoted by $Y$, of dimension $d = 4RMAX \times 4RMAX$, into a vector $\bar{Y}$ of reduced dimension, $p$ where $p << d$.

This transformation is performed using a PCA approach. The underlying principles and application steps of the PCA are described below:
## Main principle

PCA method consists in transforming data from a high-dimensional space to new data from a lower-dimensional space, with as few loss of information as possible.

<figure><img src="image_38.png" /><figcaption>FIGURE 3.3 – Visual representation of a principal component analysis (PCA)</figcaption></figure>

Formally, PCA solves the following optimisation problem :

$$ (P_1): \min_{W \in \mathbb{R}^{p \times d}, U \in \mathbb{R}^{d \times p}} \mathbb{E}(\|Y - UWY\|_2^2) $$

where :

- $Y \in \mathbb{R}^d$ is a sample of the windfield data

- $W$ is the compression matrix

- $y \in \mathbb{R}^d \rightarrow Wy \in \mathbb{R}^p$ produces a compressed version of the data, of dimension $p << d$

- $U$ is the recovery matrix

$(P_1)$ is equivalent to :

$$ (P_2): \max_{U \in \mathbb{R}^{d \times p}, U^T U = I_p} \mathrm{tr}(U^T \mathbb{E}(YY^T)U) $$

A solution to this problem is given by :

$$ U = [V_{d-p+1}, \dots, V_p] $$

where $[V_{d-p+1}, \dots, V_p]$ are the eigenvectors associated to the $p$ highest eigenvalues $\lambda_i$ ($i \in [1 : p]$) of $\mathbb{E}[YY^T]$.

In addition to that :

$$ \underbrace{\mathrm{tr}(U^T \mathbb{E}(YY^T)U)}_{\Delta} = \sum_{i=d-p+1}^{d} \lambda_i \le \sum_{i=1}^{d} \lambda_i $$

where $\Delta$ is an indicator of the quality of the approximation made when reducing the dimension of the data.
## Calculation steps

1. On our data, p is set at 3, such that the approximation quality ratio, $\frac{\sum_{i=d-p+1}^{d} \lambda_i}{\sum_{i=1}^{d} \lambda_i}$, arrives at 60%

2. For any given windfield data at time t, denoted by Y and of dimension d, the application of the compression function f is given by :

$$f(Y) = U^T Y = \begin{pmatrix} \sum_{i=1}^{d} U_{i1} y_i \\ \sum_{i=1}^{d} U_{i2} y_i \\ \sum_{i=1}^{d} U_{i3} y_i \end{pmatrix} = \begin{pmatrix} \langle PC_1, y \rangle \\ \langle PC_2, y \rangle \\ \langle PC_3, y \rangle \end{pmatrix} = \bar{Y}$$

where $(PC_i)$ ($i \in [1 : n]$) are called the principal components and $\bar{Y}$ the reduced windfield data, of dimension $p = 3$. $\bar{Y}$ is entirely defined by the input data Y and the principal components

Even if the PCA effectively reduces the number of dimensions (without losing much information, i.e. quality ratio of 60%), the flaw of this technique lies in the lack of physical interpretability of some of the new variables.

## 3.3 Prediction of reconstructed windspeed fields

After having described the reduction of windfield fields over 3 principal components, the next goal is to build windspeed fields for the simulated events based on the PCA reconstruction technique.

### 3.3.1 Reduced problem framework

<table><caption>FIGURE 3.4 – Method, data and R packages used</caption><thead><tr><td>Methodology</td><td>Data used</td><td>R Package</td></tr></thead><tbody><tr><td>Principal Component Analysis - Reconstruction</td><td>INCYC</td><td>stats (prcomp)</td></tr></tbody></table>

The reduced problem to solve is the estimation of a conditional law $Y|X$, where :

* $\bar{Y}_t$ is the compressed windfield data at time t such as

$$\bar{Y}_t = (\bar{Y}1_t, \bar{Y}2_t, \bar{Y}3_t)$$

* $X_t$ is the track vector at time t, with intensity and locations parameters simulated at Chapter 2

$$X_t = (MWS_t, RMAX_t, AMAX_t, CP_t, lon_t, lat_t, F_t, dF_t, BEAR_t, dBEAR_t)$$

The problem is solved by a quantile random forest algorithm, with detailed principles described in the Appendix H.

## Main principles

Quantile Regression Forest (QRF) [43] share similar principles than the regression forest. Both build a series of regression trees and aggregate their predictions to provide an overall prediction value. However, the quantile regression forest, instead of memorizing only the mean prediction value, retains the prediction value for each
tree. Then, a conditional distribution is obtained and quantiles of these distributions could be extracted from it.

The main advantage of quantile regression forest is that it is not restricted to predict a mean value. In the context of predicting Principal Component Analysis (PCA) weights and the variation of these values, it is useful to adjust the regression forest predictor by adding noise by randomly drawing quantiles from its distribution [42]. This impacts the variability of the simulated windfields.

## Main steps

**Trees fitting and prediction** The performance of the algorithm is measured by the MSE (mean squared error) method seen in section 2 [MSE].

Parameters are fitted to minimize the MSE and their final values are as follows :

- number of trees : 648

- maximum depth of the individual trees : 4

The final MSE obtained on the test dataset arrives at 0.013, which means the average error on the $\bar{Y}$ (weights) predicted in comparison with the windfields weights is 0.114.

### 3.3.2 Windfields reconstruction

For each timestep, the reconstruct windfield, denoted by $Y^*$ and of dimension $d$, can be reconstructed as follows :

$$Y \approx Y^* = U \underbrace{U^T Y}_{\text{compression}} = \sum_{i=1}^{p} \bar{Y}_i PC_i$$

where :

- $Y$ is the integral windfield of dimension $d$ at timestep $t$

- $Y^*$ is the reconstructed windfield of dimension $d$ at timestep $t$.

- $(Y_i)_{1 \le i \le 3}$ is the compressed windfield of dimension $p=3$, predicted at section 3.2.

- $(PC)_{1 \le i \le 3}$ are the principal components, each $PC_i$ has a dimension equal to $d$

- $U$ is defined in section 3.2

Thanks to this methodology, applied from [42], the final windfield simulated are "reconstructed" windfields. These reconstructions are based on simulated PCA compressed windfields.

Here is an illustrative example of the PCA reconstruction process. It highlights the fact that photography
or images can be considered as a matrix with coefficients represented by the pixel values.

<figure><img src="image_39.png" /><figcaption>FIGURE 3.5 – Reconstructed image</figcaption></figure>

### 3.3.3 From windfield to wind footprint generation

The wind footprint is the maximum wind over the duration of a TC event. It is calculated by taking the maximum windspeed value of each resolution squares of all the reconstructed windfields $Y_t^*$.

The wind footprint, also referred as MAXGRID, is then given by :

$$ (MAXGRID(k) = \max_t[Y_t^*(k)])_{k \in D} $$

where $D$ is the set of geographic coordinates of the NWP basin (at 1km resolution, restricted at 4RMAX around each point of the track) and $Y_t^*(k)$ the windspeed value of the reconstruction windfield at timestep $t$ and coordinates $k$.

<figure><img src="image_40.png" /><figcaption>FIGURE 3.6 – Representation of $Y^*$ transformation into MAXGRID</figcaption></figure>

On the left figure, we see three $Y^*$ at different timestep $t$ as well as three 1km resolution squares (orange-reddish colors) located at the same place belonging respectively to the three $Y_t^*$.

The figure on the right represents MAXGRID. The 1km resolution square (in red) shows the maximum of
the values taken by Y at all timestep t at this square of center (x, y). Each point of MAXGRID is obtained in this way.

## 3.4 Adjustment of wind footprints

So far, footprints do not reflect the irregularities of the land.

This section shows how wind footprints shall be adjusted to include:

- The terrain effect, included through the roughness effect which is the variable taking into account the land surface.

- The gust effect, which converts the previously calculated wind into another wind representation that correlates better with damage. The gust effect amplifies the MWS by simulating a wind speed calculated over 3 seconds (instead of 1 minute). This more conservative measure is based on the assumption that the damage to buildings correlates better with gust wind than with average winds [25].

### 3.4.1 Roughness effect

The wind speed varies according to certain parameters such as the altitude and the ground roughness (including the physical obstacles). Wind footprints until now are modeled represent wind at 10 meters above sea level ([42]). It is then necessary to adjust the values to take into account the terrain effect.

The wind footprints MAXGRID are modified to include the roughness effect by the following formula:
For every spatial pixel k:

$$ \text{Roughness Adjusted MAXGRID}[k] = \text{MAXGRID}[k] * \text{Roughnessfactor}[k] $$

where roughnessfactor[k] is given by:

$$ roughnessfactor[k] = \frac{MWS(10, l[k])}{MWS(10, l_{water})} \approx \frac{ln(10/l[k])}{ln(150/l[k])} \times \frac{ln(150/l_{water})}{ln(10/l_{water})} $$

where:

- MWS(z, l): is the MWS at height z in meters over a surface represented by the parameter l

- l[k] is the roughness length in meters, at a spatial pixel k. It is a representation of the terrain roughness. This parameter is available for any point on Earth from ERA5 database [27].

- lwater is the roughness length of water, uniformly equal to 0.002.

Note that (1) shows that roughness factor can be applied to MWS in a consistent way, and approximation
(2) is mainly based on the log wind speed formula see [49].

<figure><img src="image_41.png" /><figcaption>FIGURE 3.7 – Roughness factors over NWP area</figcaption></figure>

### 3.4.2 Gust effect

The so-called gust effect is defined as a sudden increase in wind speed above the average, which is likely to create building damage. It is represented by the maximum of 3 second-moving averages of wind observations. A Gust factor is hence used to converse the 1 min averaged MWS that was used so far, to gust wind.

$$ \text{Final Adjusted Windspeed Grid} = \text{Roughness Adjusted Windspeed Grid [k]} * \text{Gust factor} $$

where *Gust factor* = 1.49.

This factor is based on WMO report [47].
## 3.5 Footprints validation

The simulated footprints, named Final Adjusted Winsdpeed Grid [k], look like this :

<figure><img src="image_42.png" /><figcaption>FIGURE 3.8 – Modelled footprint - TC Hato</figcaption></figure>

In order to validate the Final Adjusted Winspeed Grid [k], we compare the modelled windspeed values to observed values, for a set of historical events.

The first step is to simulate wind footprints from historical events. The parameters are from the IBTRACS database.

The second step is to collect wind speed data, recorded on land at relatively low altitudes and reported by government meteorological authorities. This information is collected from two sources.

1. The NCEI (National Centres for Environmental Information) database includes globally observed winds.

2. The annual reports of the Hong Kong Observatory (HKO) have detailed wind observations associated with major neighbouring TC events.

The following QQ plots compare the distribution of observed and simulated winds, for three major events :

- CAT 3 - HATO 2017 - (26 observation points)

- CAT 4 - VICENTE 2012 - (31 observation points)

- CAT 5 - JEBI 2018 - (67 observation points)
<figure><img src="image_43.png" /><figcaption>FIGURE 3.9 – HATO, VICENTE and JEBI TCs - Observed gust versus modelled gust</figcaption></figure>

The X axis shows the observed gust values at the station location. The Y axis shows the simulated gust values at the same location. The red line is the graph of the identity function. If points are above this line, it means observed values are higher than modelled values. The color of the points is the distance in *km* of the observation location from the center of the TC, from the closest in dark color, to the farthest in light color.

It should be noted that for TC Jebi and TC Vicente, the most accurate modelled gusts (i.e the closest values to the red line) are on average located close to the center of TC (less than 100 km). This phenomenon shows that the model is more reliable when it is close to the track and seems to underestimate gust values when locations are far to the TC center. (i.e overall, light points are on the left of the graph).

We do not catch well the wind footprint variability far from the center of the TC. One possible explanation is that distant locations are more sensitive to the predicted *RMAX*. As shown in section 1, *RMAX* values are generally less than 50km. If the predicted *RMAX* is less than 50km then the wind speed decreases as the points are located further away than the predicted *RMAX*. Points close to the center (less than 100 km) are in general closer to the *RMAX* and experience less of this decrease.

Since the average *RMAX* of HATO is much smaller than Vicente and Jebi (16km vs. 30km and 26 km), the extent of the windfield footprint is also smaller. Thus the phenomenon of error amplification with respect to the distance of the cyclone is less true since the observation points are closer to the center.

Overall, even if we did not find in the literature models to compare our results with the simulated values seem to match the observed values at least for the closest locations.
## 3.6 Limits

### 3.6.1 Lack of observation data for wind footprint validation

First, observed gust data is scarce for North West Pacific geography. The networks of wind observation devices are less numerous than in the US for instance, and the data is also less shared.

Second, this data is valid for the specific locations of the wind stations, usually airports or open urban areas (like a tower summit). There is almost no observation data for mountainous areas, that are crucial for geographies like Japan or Taiwan. Then it is difficult to draw strong conclusions on the quality of the modelled gust based on such few data.

### 3.6.2 Methodology for gust and roughness effect evaluation

First, more sophisticated methods exist for better evaluation of gust and roughness effect, be it non parametrical [30] or parametrical [23]. These methods are implemented in most of commercial models and rapidly evolving. For instance, Reask provided significant improvements in that field in 2020 [41].

Second, these methodologies account for three parameters that are not included in this study : the topography, the land cover and the direction of winds.

(i) The topography effect impacts both roughness and gust evaluation. With equal roughness characteristics, the wind speed is not the same upstream and downstream of a valley.

(ii) The land cover strongly impacts the gust factor, as forestry areas or urban dense areas affect the wind differently.

(iii) The direction of the wind strongly impacts both the gust and roughness effects.

These models represent a significant additional computational cost, as the topography and land cover need to be known at a very fine resolution (≤ 1km), as well as a knowledge of wind behavior at the local scale (for the wind direction) which is difficult to achieve on a continental scale model.

## 3.7 Key takeaways

In **Chapter 3** the tracks generated in **Chapter 2** are used to predict windfields (i.e. windspeed values outside of the track points for every timestep).

The first step of **Chapter 3** is the dimension reduction of the historical windfields, as their dimension were too high compared to the one of the predictors. To do so, principal component analysis (PCA) was applied to the windfield at each timestep to reduce them to only tree values, which translates into a 40% loss of information. The result from this PCA is called the « reduced windfield ».

Second PCA reconstruction was applied to reduced windfields, to rebuild them at their initial dimension. The wind footprints are then defined as the maximum windspeed values of the reconstructed windfields, over the TC event duration.

The next step was to adjust simulated wind footprints to take into account two factors : the ground roughness and the gust effect. The first element is included using parametric formulas and ERA5 data. The second element, which is based on WMO methodology, adjusts the wind to account for the maximum of 3 seconds moving averages winds, leading to better estimation of insurance losses.

Finally, by simulating the windfields on each simulated track from previous chapter, we obtain a stochastic catalog of 10,000 years of simulated TC events characterized by their wind footprints.

Lastly, the approach used in **Chapter 3** has two main limits. First, the observations used to validate the adjusted wind footprints are scarce, and in particular not representative of mountainous areas. Second, the gust and roughness factors fail to include key parameters such as wind direction, land cover and topography
# 4. Insured Losses Generation

After **Chapter 3** shows the building of a « stochastic footprint catalogue », **Chapter 4** provides the methodology to translate each of these footprints into losses. These losses represent the economic impact on the insurer portfolio. This is the third and last step towards the building of the YLT.

## 4.1 Exposure data

The exposure data used in this study is an extraction from a P&C insurance portfolio in Philippines. We added some stochastic noise on the variables, to obtain a fictional but representative insurance portfolio.

The data contains the following information :

- Site Id, which refers to the ID of a given insured building;

- Latitude and Longitude, which refers to the geographical coordinates of the location of the insured building;

- Sum insured, which is the value of the building retained by the insurance policy;

- Coverage type is Building only. Content and Business Interruption coverages are filtered out;

- Country of the insured risks, here only Philippines;

- Line of Business, which refers to the type of business insured, here commercial, residential and industrial.

The geography of the exposure portfolio is based on the high frequency of TC in the Philippines (on average, since 1980, 3.8 landfalling TC per year, and 1.8 major landfalling TC per year) and the exposure of AXA Group.

The breakdown of the portfolio by business is shown below :

<figure><img src="image_44.png" /><figcaption>FIGURE 4.1 – Exposure portfolio breakdown by Line of Business (total insured value)</figcaption></figure>

This graph shows that insured risks attached to a residential line of business account for more than two-thirds of the sum insured values. Although industry accounts for 34% of the Philippine GDP, this is not reflected in the exposure portfolio.

## 4.2 Damage state and damage ratio

After having described the exposure portfolio on which TC risk is evaluated, we will now examine the parameters reflecting the effect of wind intensity on buildings : the damage states and the damage ratios.
### 4.2.1 Damage states

Based on Federal Emergency Management Agency (FEMA) Hazus report [25], we define five damage states that characterize the physical destruction of a building by TC wind. These five states are called DS0, DS1, DS2, DS3, DS4 for the different degrees of damage inflicted.

<table><caption>Damage states definition</caption><thead><tr><td>DS</td><td>Definition</td><td>Description</td></tr></thead><tbody><tr><td>DS0</td><td>No Damage</td><td>Little or no visible damage from the outside, no broken windows, minimal loss of roof cover with minimal water intrusion</td></tr><tr><td>DS1</td><td>Minor Damage</td><td>Maximum of one broken window or door. Moderate roof cover loss that can be covered to prevent additional water entering the building. Marks or dents of wall requiring painting or patching for repair</td></tr><tr><td>DS2</td><td>Moderate Damage</td><td>Major roof cover damage, moderate window breakage, minor roof deck failure, come resulting to interior damage from water</td></tr><tr><td>DS3</td><td>Severe Damage</td><td>Major window damage or roof sheathing loss, major roof cover loss, extensive damage to interior from water</td></tr><tr><td>DS4</td><td>Destruction</td><td>Complete roof failure: significant amount of the wall envelope opened through window failure. Extensive damage to interiors</td></tr></tbody></table>

<caption>FIGURE 4.2 – Damage states definition</caption>

These damage states are effectively observed after a TC event and reflect the difference levels of resistance of the building components. Usually, the first damage to appear at low windspeed is the loss of roof shingles or a limited number of light components (door or windows), which corresponds to damage states 0 and 1. Then, with windspeed increasing and according to the level of design of the building, the loss of more light components (multiple windows, flashing elements) can lead to partial structural damage and interior impacts due to rain intrusion. This corresponds to damage states 2 and 3. Damage state 4 usually happens when the TC wind exceeds the design windspeed for non-engineered or non-fortified buildings. More information on this can be found at [34].

### 4.2.2 Damage ratios

To each damage state DS(i) ($i \in [0:4]$), is associated a damage ratio DR(i). This damage ratio represents the ratio of the repairment cost paid by the insurer over the total building value.

<table><caption>FIGURE 4.3 – Correspondence table of damage ratio</caption><thead><tr><td>LOB</td><td>Country</td><td>Peril</td><td>Damage_state</td><td>Coverage</td><td>Damage_Ratio</td></tr></thead><tbody><tr><td>6</td><td>ALL</td><td>WS</td><td>DS0</td><td>B</td><td>0.000</td></tr><tr><td>6</td><td>ALL</td><td>WS</td><td>DS1</td><td>B</td><td>0.250</td></tr><tr><td>6</td><td>ALL</td><td>WS</td><td>DS2</td><td>B</td><td>0.500</td></tr><tr><td>6</td><td>ALL</td><td>WS</td><td>DS3</td><td>B</td><td>0.701</td></tr><tr><td>6</td><td>ALL</td><td>WS</td><td>DS4</td><td>B</td><td>1.000</td></tr></tbody></table>

While damage states relate to a physical damage, damage ratios relate to a cost. It explains that these two variables are not always proportional. Moreover, this distinction enables to model complex relationships between physical damage and repairment cost. [17]

## 4.3 Fragility curves

Fragility curves are functions that associate the gust wind speed calculated in Chapter 3 to a probability of being in a given $DS_i$, with $i \in [0:4]$.

In this study fragility functions are provided by open source Hazus model MH2.1 from FEMA [25]. They are based on building codes and Monte Carlo analyses, at component level. More information on the methodology
is given in [26].

These functions represent the so called « secondary uncertainty » in natural catastrophe models, which is the uncertainty in insured losses resulting from a TC.

An example of these fragility curves, is provided below.

<figure><img src="image_45.png" /><figcaption>FIGURE 4.4 – Fragility curves</figcaption></figure>

Above is presented a fragility function under a discrete probabilities form. For instance, for a gust between 269.68 and 280.14 km/h, the probability of reaching the damage state DS3 is 0.7094.

The probability of belonging to a damage state as a function of gust speed is not proportional. For example the probability of belonging to damage state 3 or above is 0% below 165 km/h and 99% at 250 km/h.

This can be explained by two phenomena, on the one hand the construction and prevention standards of the buildings make that they are likely to resist to average winds [10]. On the other hand, high winds generate turbulent debris clusters that increase the chances of affecting buildings [25].

## 4.4 Adapation of fragility curves

Fragility curves from FEMA [25] are based on building codes and building typologies from the United States exclusively. In Philippines, buildings have different vulnerabilities and characteristics compared to the ones from the US.

In particular, the insurance portfolio used in this study includes a large number of engineered buildings, such as office towers or condominiums that are not well represented by FEMA fragility curves. Hence, an optimisation technique is used to adapt the fragility curves to AXA claim data from wind events in the Philippines.

### Gradient descent overview

Gradient descent is an iterative optimization algorithm for finding the local minimum of a function.

The goal of the gradient descent algorithm is to minimize a cost function *J*. To achieve this goal, it performs two steps iteratively:
— Compute the gradient (slope), the first order derivative of the function at that point

— Make a step (move) in the direction opposite to the gradient by $\gamma$ times the gradient at that point.

$$a_{n+1} = a_n - \gamma \times \frac{\delta}{\delta a_n} J(a_n)$$

where:

— $J$ is the cost function

— $\gamma$ is the learning, it shows the change of direction importance

— $a_n$ is the $n$-th point of the series (a) that we want to minimise eventually.

The method used here tends to replicate the gradient descent logic by optimizing the algorithm on historical building destruction level data in Philippines from AXA.

## Application of Gradient Descent

**Data used** The data used include the fragility curves and AXA claim data mapped with wind footprints of historical events. This data called AXA vulnerability historical observations is composed of four parameters : lon, lat, MWS<sub>(lon,lat)</sub>, DS<sub>(lon,lat)</sub>.

**Fragility curve adjustment based on MWS<sub>(lon,lat)</sub>** We denote historical observations $X = (X_1, ..., X_T)$ where $T$ is the number of historical observations. $\forall j \in [1; T]$ :

$$X_j = (lon_j, lat_j, MWS_{(lon,lat)}(j), DS_{(lon,lat)}(j))$$

For each $X_j$ we follow these steps:

1. Find the two wind intensity bounds that frame $MWS_{(lon,lat)}(j)$ in the fragility curve table 4.4. We denote them by $B_1$ and $B_2$.

2. Compute weights $w_1$ and $w_2$ associated to the interval bounds. For example, if the wind gusts have reached 177 km/h, the two gust wind speed bounds of 174 and 182 km/h ($B_1$ and $B_2$) have the following weights:

$$\left\{ \begin{array}{l} w_1 = \frac{|177-174|}{182-174} = 0.375 \\ w_2 = \frac{|177-182|}{182-174} = 0.625 \end{array} \right.$$

3. The probabilities in the fragility curve table 4.4 are then updated for the two intensity intervals as follows. $\forall i \in [0, 4], \forall j \in [1; T]$ :

$$\mathbb{P}(DS_i(j) | MWS_{(lon,lat)}(j) \in B_1) = \mathbb{P}(DS_i(j-1) | MWS_{(lon,lat)}(j) \in B_1) - \gamma \times w_1 \times [\mathbb{P}(DS_i(j-1) | MWS_{(lon,lat)}(j) \in B_1) - 1_{(DS_{(lon,lat)}(j)=DS_i(j-1))}]$$

$$\mathbb{P}(DS_i(j) | MWS_{(lon,lat)}(j) \in B_2) = \mathbb{P}(DS_i(j-1) | MWS_{(lon,lat)}(j) \in B_2) - \gamma \times w_2 \times [\mathbb{P}(DS_i(j-1) | MWS_{(lon,lat)}(j) \in B_2) - 1_{(DS_{(lon,lat)}(j)=DS_i(j-1))}]$$

Finally, we keep the values of $(\mathbb{P}(DS_i(T)|B))_{i \in [1;4]}$ for each bound $B$, to obtain the optimized fragility curve.

This corresponds to an application of the gradient descent with:
— J, the cost function would equals the functions : ∀B, ∀i ∈ [1; 4], ∀j ∈ [1; T] :

$$MWS_{(lon,lat)}(j) \in X :\rightarrow \frac{\left( \mathbb{P} (DS_i(j-1) | MWS_{(lon,lat)}(j) \in B) - 1_{MWS_{(lon,lat)}(j)=DS_i(j-1)} \right)^2}{2}$$

The main difference is that J is not differentiable given the finite interval of T points.

— ∀k ∈ [1; 2], the learning rate equal to γ * w_k

## 4.5 Simulation procedure for Year Loss Table (YLT) and Occurrence Exceedance Probability (OEP) derivation

We will now group these simulated economic losses for the stochastic event catalog into a Year Loss Table and then an Occurrence Exceedance Probability table and curve, which have been described in Chapter 1.

In order to derive a Year Loss Table, the simulation procedure is as follows :

— Step 1 : For a given simulation year in the stochastic catalogue footprint derived in Chapter 3, identify all the events pertaining to this year and get their respective wind footprints.

E.g : for Simulation Year 1 and Footprint with ID = 1

<figure><img src="image_46.png" /></figure>

— Step 2 : For each of these footprints, associate each insured risk of the exposure data to a gust windspeed thanks to its latitude longitude information.

E.g : for Simulation Year 1 and Footprint with ID = 1

<table><thead><tr><th>Site ID</th><th>Sum Insured</th><th>Latitude</th><th>Longitude</th><th>Gust Windspeed</th></tr></thead><tbody><tr><td>1</td><td>100 000</td><td>0</td><td>100</td><td>150</td></tr><tr><td>2</td><td>200 000</td><td>1</td><td>102</td><td>120</td></tr></tbody></table>

— Step 3 : For each risk, given the gust windspeed derived at Step 2, draw a damage state DS in the fragility matrices and take the associated Damage Ratio (DR). Multiply this DR by the sum insured of the risk.

E.g : for Simulation Year 1, Footprint with ID = 1 and sites ID 1 and 2 from the exposure data

<table><thead><tr><th>Site ID</th><th>Sum Insured</th><th>Latitude</th><th>Longitude</th><th>Gust Windspeed</th><th>DS selected</th><th>DR associated</th><th>Loss</th></tr></thead><tbody><tr><td>1</td><td>100000</td><td>0</td><td>100</td><td>150</td><td>DS2</td><td>0.50</td><td>50000</td></tr><tr><td>2</td><td>200000</td><td>1</td><td>102</td><td>120</td><td>DS1</td><td>0.25</td><td>50000</td></tr></tbody></table>
— Step 4 : Sum all the losses of each risk, event by events.

E.g : for Simulation Year 1, Footprint with ID = 1 and sites ID 1 and 2 from the exposure data, the loss is 50000 + 50000 = 100000 and the YLT includes the following information :

<table><thead><tr><td>Year</td><td>Event ID</td><td>Loss</td></tr></thead><tbody><tr><td>1</td><td>1</td><td>100000</td></tr></tbody></table>

Based on the YLT, we derive an OEP curve with the corresponding destruction rates, for various return periods. This is shown below :

<figure><img src="image_47.png" /><figcaption>FIGURE 4.5 – OEP curve for TC risk - Philippines</figcaption></figure>

The destruction rate of the portfolio at the 99.5 quantile is 0.61%. This means that if we have an exposure portfolio of EUR 1 Bn, the modeled loss on a TC whose frequency of occurrence is (on average) once every 200 years is EUR 61 Mn.

It is difficult to compare this destruction rate with those of insurers because it is considered sensitive information and cannot be found publicly. Moreover, the destruction rate at the 99.5 quantile highly varies depending on the concentration of the exposure portfolio and the heterogeneity of the covered area.

## 4.6 Towards a framework for climate change impact assessment

The previous chapters showed how the simulation of key parameters such as forward speed (F) or maximum sustained windspeed (MWS) were used to eventually generate economic losses over an insured portfolio. Recent scientific literature demonstrated that climate change impacts some of these key parameters in a measurable way. Hence, we present here some framework to calculate the financial implication of climate change, by adjusting these parameters. The choice for this parameters, and the quantification of this adjustment are explained below :

### Frequency of high categories TC events

Two primary factors dictate how strong a TC wind event can be : Ocean surface temperature and moisture state of the atmosphere. Climate change increases these two factors with a high confidence and we can therefore expect a shift towards stronger storms.

Some researchers [40] support this expectation and mention a global increase of 25% of the probability for
a TC event to reach a major category (between 3 and 5).

**Proposed adjustment** : In section 2, adjust the simulation table by substituting CAT 0, CAT 1, CAT 2 TC events by CAT 3, CAT 4 and CAT 5 TC events.

The substitution of minor TCs into major TCs is done according to a frequency increase assumption of 25%. For instance, with this assumption, and with 132 major cyclones in a simulation table, we can transform 33 (0.25 × 132) minor TCs into major TCs.

Thus, in section 2.1 we obtain modified stochastic event blocks. Instead of having a three-year block like the one below :

— **Simulation 1** : Event 1 - CAT 2; Event 2 - **CAT 1**; Event 3 - CAT 4

— **Simulation 2** : Event 1 - CAT 5; Event 2 - CAT 3

— **Simulation 3** : Event 1 - CAT 4; Event 2 - CAT 5; Event 3 - CAT 4; Event 4 - CAT 2

You get a three-year block with potentially more minor TCs.

— **Simulation 1** : Event 1 - CAT 2; Event 2 - **CAT 4**; Event 3 - CAT 4

— **Simulation 2** : Event 1 - CAT 5; Event 2 - CAT 3

— **Simulation 3** : Event 1 - CAT 4; Event 2 - CAT 5; Event 3 - CAT 4; Event 4 - CAT 2

## Global frequency of TC events

There is no clear evidence that global TC genesis rates will increase due to climate change. The consensus is even frequency of all TC events will not change or decrease.

The method described above to model the increase of major cyclones, takes into account an unchanged frequency of events.

## Forward speed

Some literature points out a possible decrease in forward speed of TC events over land. Sally and Eta in 2020, are recent examples of TC events stalling on the American coast. [39] mentions a -20% change in the mean forward speed for North West Pacific basin.

**Proposed adjustment** : In section 2.2, we adjust by a 20% factor the Forward speed values from IBTRACS database.

## Track shift

TC event tracks are expected by some authors to move polewards, this shift being especially strong in the North West Pacific according to [37] and [38].

On average, more TC events tend to reach their *MWS* (maximum wind speed) at latitude over 30 degree,
which puts more risk on high latitude countries such as Japan, Korea or North East China.

**Proposed adjustment:**

- In section 2.3, we adjust 50% of the *dBEAR* values from IBTRACS database by a -5 radians value.

- In section 2.1, we resample a base of starting origin points with coordinates further north.

## Rapid intensification

[19] shows that the chance of having a TC event intensification by 35 mph or more in a 24h period has been multiply by 5 in the North Atlantic basin, partially due to climate change.

**Proposed adjustment:** Restrict the LSTM learning base to the most recent 12 years [19] where the phenomenon is visible.

A level of confidence on the expected change on these variables is given by recent literature review [24].

## 4.7 Limits

**Uncertainty around damage ratios** Damage ratios are based on US repairment cost statistics and do not reflect local practices in the Philippines. Moreover, damage states are associated to a unique damage ratios. Hence they don't account for uncertainties in the repairment costs, for a given physical damage situation.

**Non modelled insured losses** The modelling approach does not include key elements of insured losses caused by TC wind:

* Additional Living expenses (ALE) and Business interruption (BI) that account for the cost of interrupting the business of companies and relocating people. These two ones can happen, even at very low windspeed, due to downtimes of critical infrastructure or due to administrative closure. An example of BI modelling can be found at [35].

* The inflation of reconstruction cost due the high demand on labour and raw material markets after the event. This is often referred as « demand surge ».

* The interruption of supply chains, often referred as CBI.

## 4.8 Key takeaways

First of all, we presented a portfolio of insured buildings in Philippines as the exposure basis for the economic loss calculations. This geography is indeed highly prone to TC events and AXA group has insured interests there.

We then introduced the damages states (DS) and defined them as 5 consecutive steps of the physical destruction of buildings due to the impact by wind. They range from DS0, no damage, to DS4, complete destruction. To each damage state is associated a DR (Damage Ratio). The DR is the cost of reparation paid by the insurer divided by the sum insured of the building.

Afterwards, we defined the fragility curves as the probability for a building to be in a given damage state,
given a certain wind intensity. Hence, fragility curves provide a mapping between a windspeed and a couple of (damage state, damage ratio).

Since the fragility curves are suitable for American buildings only, we adjusted these curves to the Philippines exposure using a gradient descent method and local claim data.

Based on the intersection between the exposure database, the stochastic catalogue of footprints and the adjusted fragility curves, we produced a YLT. This YLT was used to derive OEP curves and Destruction rates (DR) at various return periods. For the 200 years return period, the DR arrives at 0.61%

Following these results, we investigated how the model could be adapted to account for climate change.

We noticed that adjustments of IBTRACS enables the inclusion of some predicted effects such as :

- The increase in frequency of extreme TC events

- The decrease in Forward Speed (F)

- The poleward movement of tracks
# 5. Conclusion and Perspectives

This study has presented a complete modelling framework for the evaluation of TC Wind risk over an insurance portfolio, with an application to the North West Pacific basin.

**Chapter 1** described the data used and the key parameters used to represent this risk.

**Chapter 2** and **Chapter 3** examined the modelling and provided two central outputs being tracks and footprints:

Those outputs are the result of a four step process. First, the track generation using Markov Chains, second the association to this trajectory of the main characteristics of the cyclone (RMAX, MWS, CP, AMAX) through machine learning algorithms, third the generation of wind footprints around these tracks, fourth the adjustment of these footprints to include roughness and gust effects.

This approach has shown some limitations as it heavily relies on historical data and leaves out some key variables such as the Sea Surface Temperature. Nonetheless, this process has shown great effectiveness and speed to produce a hazard catalogue of 10,000 years of simulation, with a quality equalling commercial software at some stages.

**Chapter 4** shows the intersection of this catalogue with a representative insurance portfolio in the Philippines, leading to the derivation of standard outputs of catastrophe models such as YLT, OEP and destruction rates.

This shows and let us conclude that developing in-house catastrophe models for TC wind risk is possible and can constitute an alternative to catastrophe commercial software. It brings many advantages, such as licence cost decreases and transparency in the modelling. It also brings challenges, such as IT maintenance or code versioning management. But some say you shall take the rough with the smooth. And with the rise of climate change impacts calculation requests, we conclude with the smooth:

The internal methodology presented here enables a flexible and transparent inclusion of climate change effects, described in **Chapter 4**. Interestingly enough, some other actors from the finance sector also followed this route [29] to evaluate these effects. We hope this study will convince others to join.
Bibliographie

[1] S. Athawes. Five years on : How haiyan shocked the world. www.rms.com, 2018.

[2] G. Atkinson and C. Holliday. Tropical cyclone minimum sea level pressure/maximum sustained wind relationship for the western north pacific. *Monthly Weather Review*, Volume 105, 1977.

[3] International best track archives IBTRACS. Tropical cyclones by year. www.ncei.noaa.gov, 2021.

[4] N. Bloemendaal, I. Haigh, and H. De Moel et al. Generation of a global synthetic tropical cyclone hazard dataset using storm. *Nature, Scientific Data*, Volume 7, 2020.

[5] L. Breiman. Random forest. *Machine Learning*, Volume 45, 2001.

[6] NASA's Goddard Space Flight Center. Nasa adds up heavy rains from typhoon jebi. phys.org/news/2018-09-nasa-heavy-typhoon-jebi.html, 2018.

[7] National center for atmospheric research(NCAR). Weather research and forecasting model (wrf). ral.ucar.edu, 2007.

[8] K.A Emanuel C.Risi. Statistical synthesis of tropical cyclone tracks in a risk evaluation perspective. www.lmd.jussieu.fr, 2006.

[9] C. Denise-Baillon. Modélisation statistique du risque de tempêtes tropicales dans le bassin atlantique nord. www.institutdesactuaires.com, 2013.

[10] Florida department of business and professional regulation. Florida building code. floridabuilding.org, 2020.

[11] T. Duong. Using ks for bivariate kernel density estimation. rdrr.io/cran/ks/f/inst/doc/kde.pdf, 2006.

[12] T. Duong. ks: Kernel density estimation for bivariate data. cran.r-project.org/web/packages/ks/index.html, 2021.

[13] V. Dvorak. Tropical cyclone intensity analysis and forecasting from satellite imagery. Monthly Weather Review, Volume 103, 1975.

[14] EM-DAT. 2020 - the non-covid year in disasters - global trends and perspectives. Centre for Research on the Epidemiology of Disasters, UN Office for Disaster Risk Reduction, 2021.

[15] K.A. Emanuel. The dependence of hurricane intensity on climate. Nature, Volume 326, 1987.

[16] K.A. Emanuel. Tropical cyclones. Annual Review of Earth and Planetary Sciences, Volume 31, 2003.

[17] A. Pothon et al. Assessing the performance of existing repair-cost relationships for buildings. 16ECEE, 2018.

[18] K.A. Emanuel et al. A statistical deterministic approach to hurricane risk assessment. Bulletin of the American Meteorological Society, Volume 87, 2006.

[19] K.Bahtia et al. Recent increases in tropical cyclone intensification rates. Nature Communication 10, 635, 2019.

[20] Knapp et al. The international best track archive for climate stewardship. Bulletin of the American Meteorological Society, Volume 91, 2010.

[21] L. Li et al. Field survey of typhoon hato (2017) and a comparison with storm surge modeling in macau. Natural Hazards Earth System Sciences, Volume 18, 2018.

[22] M. Lang et al. Circular regression trees and forests with an application to probabilistic wind direction forecasting. Journal of the Royal Statistical Society Series C Applied Statistics, Volume 69, 2020.

[23] S. Balica et al. Parametric and physically based modelling techniques for flood risk and vulnerability assessment: a comparison. Environmental Modelling and Software, Volume 41, 2013.

[24] T. Knutson et al. Tropical cyclones and climate change assessment. Bulletin of the American Meteorological Society, Volume 101, 2020.

[25] Federal Emergency Management Agency (FEMA). Understanding the hazards. www.fema.gov/pdf/library/ism2s1.pdf, 2000.
[26] Federal Emergency Management Agency (FEMA). Hurricane vulnerability. www.fema.gov/flood-maps/products-tools/hazus, 2005.

[27] European Centre for Medium-Range Weather Forecasts (ECMWF). www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5.

[28] G.Biau and E. Scornet. A random forest guided tour. arxiv.org/pdf/1511.05741.pdf, 2016.

[29] T. Le Guenedal, P.Dobrinski, and P.Tankov. Measuring and pricing cyclone-related physical risk under changing climate. Amundi - Working Paper 111, 2021.

[30] N. Hannah. A global tropical cyclone insights database. reask.earth, 2019.

[31] B. Harper. Tropical cyclone parameter estimation in the australian region : Wind-pressure relationships and related issues for engineering planning and design - a discussion paper. J0106-PR003E - Systems Engineering Australia Pty Ltd, 2002.

[32] S. Hochreiter and J.Schmidhuber. Long short-term memory. Neural computation, Volume 9, 1997.

[33] G.J Holland. An analytical model of the wind and pressures profiles in the hurricanes. Monthly Weather Review, Volume 108, 1980.

[34] IBHS. Website of insurance institute of business and home safety. ibhs.org.

[35] V. Jain and J. Guin. Modeling business interruption losses for insurance portfolios. www.researchgate.net, 2009.

[36] S. Kang, K. Jun, and K.Park. Empirical relationship between maximum sustained wind and central pressure for typhoons in northeast asian sea. Asian and Pacific coasts, editors : Y.Goda and W.Kikoa and K.Nadaoka, 2003.

[37] J. Kossin. The poleward migration of the location of tropical cyclone maximum intensity. Nature, Volume 509, 2014.

[38] J. Kossin. Past and projected changes in western north pacific tropical cyclone exposure. Journal of Climate, Volume 29, 2016.

[39] J. Kossin. A global slowdown of tropical-cyclone translation speed. Nature, Volume 558, 2019.

[40] J. Kossin, K. Knapp, T. Olander, and C. Velden. Global increase in major tropical cyclone exceedance probability over the past four decades. PNAS, Volume 117, 2020.

[41] T. Loridan. Reask's unified tc model : a global, climate-connected view of tropical cyclone risk. reask.earth/2020/10/14/unified-tropical-cyclone-model-a-global-climate-connected-view-of-tropical-cyclone-risk/, 2020.

[42] T. Loridan, R. Crompton, and E. Dubossarky. A machine learning approach to modeling tropical cyclone wind field uncertainty. Monthly Weather review, Volume 145, 2017.

[43] N. Meinshausen. Quantile regression forests. Journal of Machine Learning Research, Volume 7, 2006.

[44] C. Molnar. Machine learning interpretability. christophm.github.io/interpretable-ml-book/interaction.html, 2021.

[45] National Hurricane Center (NHC). Glossary of nhc terms. www.nhc.noaa.gov/aboutgloss.shtml, 2021.

[46] University of Rhode Island. Glossary of terms - wind shear. www.hurricanescience.org, 2020.

[47] World Meteorological Organization. Guidelines for converting between various wind averaging periods in tropical cyclone conditions. library.wmo.int, 2010.

[48] C. Planche. On the realism of the rain microphysics representation of a squall line in the wrf model. part ii : Sensitivity studies on the rain drop size distributions. Weather monthly review, Volume 147, 2019.

[49] H. Tennekes. The logarithmic wind profile. Journal of the Atmospheric Sciences, Volume 30, 1973.

[50] E. Willoughby, R. Darling, and M. Rahn. Parametric representation of the primary hurricane vortex. part ii : A new family of sectionally continuous profiles. Monthly Weather Review, Volume 134, 2006.