In [2]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
In [3]:
df = pd.read_csv(r'SigmaInsuredNatLosses.csv')
In [12]:
print(df.tail(10))
    Year  Loss
20  2015    38
21  2016    62
22  2017   185
23  2018   106
24  2019    66
25  2020   103
26  2021   130
27  2022   138
28  2023   115
29  2024   137
In [4]:
years = df.iloc[:, 0]
values = df.iloc[:, 1]

plt.figure(figsize=(12, 6))
plt.bar(years, values, alpha=0.7)

x_numeric = np.arange(len(years))
z = np.polyfit(x_numeric, np.log(values + 1), 1)
trend = np.exp(z[1]) * np.exp(z[0] * x_numeric)
plt.plot(years, trend, 'r-', linewidth=2, label='Exponential Trend')

plt.xlabel('Year')
plt.ylabel('Losses')
plt.title('Natural Catastrophe Losses by Year')
plt.legend()
plt.tight_layout()
plt.show()

We can see we've got a pretty good fit against the Swiss Re chart. We can also guess where the 2025 estimate comes from. Let's project the exponential trend line out and compare to their 2025 value

In [6]:
x_2025 = len(years)
value_2025 = np.exp(z[1]) * np.exp(z[0] * x_2025)
print(f"2025 prediction: {value_2025}")
2025 prediction: 146.05338915565198

When we hover over 2025 in the Sigma chart, we see their value is 145. So close enough!

Now let's calculate the trend still remaining in the data, and apply this to on-level the past years.

In [8]:
ln_values = np.log(values)
slope = np.polyfit(years, ln_values, 1)[0]
result = np.exp(slope) - 1
print(f"exp(slope) - 1 = {result}")
exp(slope) - 1 = 0.06544563714388674
In [9]:
df['Revalued'] = values * (1 + result) ** (2024 - years)
print(df)
    Year  Loss    Revalued
0   1995    30  188.594412
1   1996    20  118.006591
2   1997    12   66.454779
3   1998    31  161.129616
4   1999    56  273.193531
5   2000    16   73.260701
6   2001    20   85.950773
7   2002    26  104.872554
8   2003    31  117.359674
9   2004    79  280.706842
10  2005   166  553.608437
11  2006    21   65.732872
12  2007    37  108.701051
13  2008    65  179.231394
14  2009    34   87.993045
15  2010    67  162.746980
16  2011   170  387.575005
17  2012    90  192.583047
18  2013    51  102.426993
19  2014    39   73.515271
20  2015    38   67.230332
21  2016    62  102.953722
22  2017   185  288.330640
23  2018   106  155.057807
24  2019    66   90.615066
25  2020   103  132.727956
26  2021   130  157.230660
27  2022   138  156.654068
28  2023   115  122.526248
29  2024   137  137.000000
In [10]:
plt.figure(figsize=(12, 6))
plt.bar(years, df['Revalued'], alpha=0.7)

x_numeric = np.arange(len(years))
z_new = np.polyfit(x_numeric, np.log(df['Revalued'] + 1), 1)
trend_new = np.exp(z_new[1]) * np.exp(z_new[0] * x_numeric)
plt.plot(years, trend_new, 'r-', linewidth=2, label='Exponential Trend')

plt.xlabel('Year')
plt.ylabel('Revalued Losses')
plt.title('Revalued Natural Catastrophe Losses by Year')
plt.legend()
plt.tight_layout()
plt.show()

We can see that applying this trend factor has resulted in a dataset which has no exponential trend. A couple of interesting points - The big peak in 2005 is largely Katrina, which now sit's at a around a 550bn year. 2011 had a number of large earthquakes. Now let's do 2 things, calculate an empirical 1 in 10 OEP, and also fit a basic lognormal curve Since we've got 30 years of data, we can look at the 3rd worse year as our 1-in-10

In [12]:
third_worst = df['Revalued'].nlargest(3).iloc[2]
print(f"3rd worst year: {third_worst}")
3rd worst year: 288.3306399453474

So we see, that our value comes to 288, which corresponds to a revalued 2017 year.

Now let's build our lognormal model. We're just going to fit it using a method of moments.

In [14]:
mean_val = df['Revalued'].mean()
std_val = df['Revalued'].std()

mu = np.log(mean_val**2 / np.sqrt(std_val**2 + mean_val**2))
sigma = np.sqrt(np.log(1 + std_val**2 / mean_val**2))

return_periods = [2, 5, 10, 25, 50, 100, 200]
exceedance_probs = [1 - 1/rp for rp in return_periods]

percentiles = [stats.lognorm.ppf(p, sigma, scale=np.exp(mu)) for p in exceedance_probs]

results = pd.DataFrame({'Return Period': return_periods, 'Exceedance Prob': exceedance_probs, 'Loss': percentiles})
print(results)
   Return Period  Exceedance Prob        Loss
0              2            0.500  133.024872
1              5            0.800  221.457096
2             10            0.900  289.065790
3             25            0.960  384.047419
4             50            0.980  461.417552
5            100            0.990  544.240675
6            200            0.995  633.005794

Actually our 1-in-10 as per the lognormal matches pretty closely to the empirical percentile.