2.8. Catastrophe Modeling
Objectives: Applications of the Aggregate class to catastrophe risk evaluation and pricing using thick-tailed Poisson Pareto and lognormal models, including occurrence and aggregate PMLs (OEP, AEP) and layer loss costs. Covers material on CAS Parts 8 and 9.
Audience: Catastrophe modelers, reinsurance actuaries, and risk management professionals.
Prerequisites: Basics of catastrophe modeling, catastrophe insurance and reinsurance terminology, use of build.
See also: Capital Modeling and Risk Management, Strategy and Portfolio Management, Reinsurance Pricing, Individual Risk Pricing.
Contents:
2.8.1. Helpful References
2.8.2. Jewson’s US Wind PML Estimates
2.8.2.1. Model Description
Stephen Jewson Projections of Changes in U.S. Hurricane Damage Due to Projected Changes in Hurricane Frequencies (submitted, under peer review), Jewson [2022] reports the following frequency and severity statistics for US hurricane losses.
The dataframe jewson recreates the table and adds severity CVs.
In [1]: from aggregate import build, qd, mv
In [2]: import pandas as pd
In [3]: import matplotlib.pyplot as plt
In [4]: jewson = pd.DataFrame(
...: {'Num': [197, 84, 47, 43, 20, 3],
...: 'EN': [1.67, 0.71, 0.4, 0.36, 0.17, 0.025],
...: 'ES_W': [10.0, 2.28, 4.46, 13.0, 43.8, 46.5],
...: 'ES_M': [15.9, 2.96, 6.39, 17.9, 82.3, 55.2],
...: 'SD_W': [24.4, 8.63, 6.17, 21.9, 50.9, 51.5],
...: 'SD_M': [47.2, 9.62, 7.83, 29.9, 119.0, 60.1],
...: 'EX_W': [16.7, 1.63, 1.78, 4.73, 7.42, 1.18],
...: 'EX_M': [26.5, 2.11, 2.55, 6.52, 13.9, 1.4]},
...: index=pd.Index(['1-5', '1', '2', '3', '4', '5'],
...: dtype='object', name='Cat')
...: )
...:
In [5]: jewson['CV_W'] = jewson.SD_W / jewson.ES_W; \
...: jewson['CV_M'] = jewson.SD_M / jewson.ES_M
...:
In [6]: qd(jewson)
Num EN ES_W ES_M SD_W SD_M EX_W EX_M CV_W CV_M
Cat
1-5 197 1.67 10 15.9 24.4 47.2 16.7 26.5 2.44 2.9686
1 84 0.71 2.28 2.96 8.63 9.62 1.63 2.11 3.7851 3.25
2 47 0.4 4.46 6.39 6.17 7.83 1.78 2.55 1.3834 1.2254
3 43 0.36 13 17.9 21.9 29.9 4.73 6.52 1.6846 1.6704
4 20 0.17 43.8 82.3 50.9 119 7.42 13.9 1.1621 1.4459
5 3 0.025 46.5 55.2 51.5 60.1 1.18 1.4 1.1075 1.0888
Jewson models aggregate losses with Poisson frequency and lognormal severity assumptions.
Use build to create Aggregate models of the two implied distributions. Adjust bs from recommended 1/16 to 1/8 for thick tailed distributions.
In [7]: w = build('agg Cat:USWind:W '
...: f'{jewson.loc["1":"5", "EN"].to_numpy()} claims '
...: f'sev lognorm {jewson.loc["1":"5", "ES_W"].to_numpy()} '
...: f'cv {jewson.loc["1":"5", "CV_W"].to_numpy()}'
...: 'poisson'
...: , bs=1/8)
...:
In [8]: m = build('agg Cat:USWind:M: '
...: f'{jewson.loc["1":"5", "EN"].to_numpy()} claims '
...: f'sev lognorm {jewson.loc["1":"5", "ES_M"].to_numpy()} '
...: f'cv {jewson.loc["1":"5", "CV_M"].to_numpy()}'
...: 'poisson'
...: , bs=1/8 )
...:
In [9]: qd(w)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1.665 0.77498 0.77498
Sev 10.025 10.024 -4.5361e-05 2.4845 2.4841 9.9667 9.7215
Agg 16.691 16.691 -4.5361e-05 2.0756 2.0753 6.9539 6.8019
log2 = 16, bandwidth = 1/8, validation: fails sev skew, agg skew.
In [10]: mv(w)
mean = 16.6913
variance = 1200.188
std dev = 34.6437
In [11]: qd(m)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1.665 0.77498 0.77498
Sev 15.899 15.899 -4.9864e-05 3.0261 3.0218 15.065 14.399
Agg 26.473 26.471 -4.9871e-05 2.4699 2.4668 10.676 10.23
log2 = 16, bandwidth = 1/8, validation: fails sev cv, agg cv.
In [12]: mv(m)
mean = 26.4726
variance = 4275.28
std dev = 65.3856
Plots of the severity and aggregate distributions confirms they are very thick tailed.
In [13]: w.plot()
In [14]: m.plot()
2.8.2.2. Aggregate PML Estimates
It is easy to compute aggregate PML points (aggregate quantiles). The next table shows values at a range of return periods. The return period corresponding to a \(p\) quantile is \(1/(1-p)\). In a Poisson frequency model, the reciprocal of the frequency equals the average waiting time between events because of the relationship between the Poisson and exponential distributions. Amounts are in USD billions.
In [15]: agg_pmls = pd.DataFrame({'Return': [2, 5, 10, 20, 25, 50, 100, 200, 250, 1000, 10000]}, dtype=float)
In [16]: agg_pmls['p'] = 1 - 1/agg_pmls.Return
In [17]: agg_pmls['Weinkle'] = [w.q(i) for i in agg_pmls.p]
In [18]: agg_pmls['Martinez'] = [m.q(i) for i in agg_pmls.p]
In [19]: agg_pmls = agg_pmls.set_index(['Return'])
In [20]: qd(agg_pmls)
p Weinkle Martinez
Return
2.0 0.5 4.5 6.375
5.0 0.8 23.5 33
10.0 0.9 46 68.5
20.0 0.95 73.875 117.38
25.0 0.96 84.125 136.12
50.0 0.98 119.25 204.12
100.0 0.99 160.38 288.88
200.0 0.995 208.12 392.88
250.0 0.996 225.12 431
1000.0 0.999 350.5 727
10000.0 0.9999 657.88 1516
2.8.2.3. Occurrence PML Estimates
Occurrence PMLs, called occurrence exceeding probability (OEP) points, can be computed for a compound Poisson model as adjusted severity quantiles. The \(n\) year OEP is defined as the loss level \(\mathit{OEP}(n)\) so that
there is a \(1/n\) chance of one or more losses greater than \(\mathit{OEP}(n)\) per year.
The definition is valid only for \(n\ge 1\) since \(1/n\) is interpreted as a probability. If \(\lambda\) is the annual event frequency and \(S\) the severity survival function, then the annual frequency of losses greater than \(x\) equals \(\lambda S(x)\) and therefore chance of one or more losses greater than \(x\) equals \(1-\exp(-\lambda S(x))\) with Poisson frequency (one minus chance of no events). Rearranging gives
where \(q\) is the severity distribution quantile function.
Jewson [2022] considers the related notion of event exceedance frequency (EEF). Here, the \(n\) year EEF is defined as the loss level \(\mathit{EEF}(n)\) with a \(1/n\) annual frequency. Thus
This definition is valid for any \(n > 0\) since the result is a frequency rather than a probability. OEP and EEF are very similar for large \(n\) because \(\log(1+x)\approx x\) for small \(x\), but they diverge significantly for small \(n\) and only the latter makes sense for \(0 < n < 1\). Jewson shows EEFs in his Figure 2. See Aggregate and Occurrence Probable Maximal Loss and Catastrophe Model Output.
Jewson comments that OEP is useful for validating the frequency of annual maximum loss, which is affected by clustering. Thus OEP estimates are important for validating models that include clustering. EEF is useful for validating whether a model captures the mean frequency of events, including events with frequency greater than 1 per year.
The following table shows OEP and EEF points, comparing the two statistics.
In [21]: oep = pd.DataFrame({'Return': [2, 5, 10, 20, 25, 50, 100, 200, 250, 1000, 10000]}, dtype=float); \
....: oep['p'] = 1 - 1/oep.Return; \
....: oep['W OEP'] = [w.q_sev(1 + np.log(i) / w.n) for i in oep.p]; \
....: oep["W EEF"] = [w.q_sev(1 - 1 / i /w.n) for i in oep.Return]; \
....: oep['M OEP'] = [m.q_sev(1 + np.log(i) / m.n) for i in oep.p]; \
....: oep["M EEF"] = [m.q_sev(1 - 1 / i /m.n) for i in oep.Return]; \
....: oep = oep.set_index(['Return']); \
....: qd(oep)
....:
p W OEP W EEF M OEP M EEF
Return
2.0 0.5 3.625 6.375 5.125 8.875
5.0 0.8 18.75 21.125 26.125 29.75
10.0 0.9 37.75 39.375 56.75 59.375
20.0 0.95 62.125 63.125 100.5 102.38
25.0 0.96 71.25 72 117.62 119.25
50.0 0.98 103.12 103.62 181 182
100.0 0.99 141.62 141.88 261.75 262.38
200.0 0.995 187.38 187.62 362.75 363.12
250.0 0.996 203.88 204 400.12 400.38
1000.0 0.999 327.38 327.5 693.38 693.5
10000.0 0.9999 635.38 635.38 1482.8 1482.8
The next block of code shows the same information as Jewson’s Figure 2. It includes OEP and EEF for comparison. The dashed line shows a third alternative aggregate implementation using the exact continuous weighted severity survival function m.sev.sf, rather than the discrete approximation in m.density_df.
In [22]: fig, axs = plt.subplots(1, 2, figsize=(2 * 3.5, 2.45), constrained_layout=True)
In [23]: for ax, mw, title, xmax in zip(axs.flat[::-1], [m, w], ['Martinez estimates', 'Weinkel estimates'], [550, 240]):
....: bit = np.exp(-(1-mw.density_df.F_sev.loc[:1000]) * m.n)
....: bit = 1 / (1 - bit)
....: bit.plot(logy=True, ax=ax, label='OEP = Pr no events in year')
....: bit = 1 / ((1 - mw.density_df.F_sev.loc[:20000]) * m.n)
....: bit.plot(logy=True, ax=ax, label='EEF RP = 1/freq')
....: xs = np.linspace(0, 500, 501)
....: if title[0] == 'M':
....: rp = 1 / (m.n * m.sev.sf(xs))
....: else:
....: rp = 1 / (m.n * w.sev.sf(xs))
....: ax.plot(xs, rp, c='r', ls='--', lw=2, alpha=0.25, label='From sev sf')
....: ax.legend()
....: ax.set(xlim=[-10, xmax], ylim=[400, .5], title=title)
....: ax.set_yscale('log', base=10)
....: ticks = [1, 3, 10, 30, 100, 300]
....: ax.set_yticks(ticks)
....: ax.set_yticklabels([f'{x}' for x in ticks]);
....:
In [24]: fig;
2.8.2.4. Feller’s Relationship between AEP and OEP
For thick tailed distributions, AEP and OEP points are closely related by Feller’s theorem, which says that for \(A\sim \mathsf{CP}(\lambda, X)\) with severity \(X\) subexponential,
as \(x\to\infty\), see REF. The next plot confirms that Feller’s approximation is very good. Note the extreme return periods returned by aggregate that would be hard to estimate with simulation.
In [25]: fig, axs = plt.subplots(1, 2, figsize=(2 * 3.5, 2.55), constrained_layout=True)
In [26]: for ax, mw, lim, title in zip(axs.flat[::-1], [m, w], [5000, 5000], ['Martinez', 'Weinkle']):
....: bit = mw.density_df.loc[:5000, ['S', 'S_sev']]
....: bit['Feller'] = bit.S_sev * mw.n
....: bit = 1 / bit
....: bit.plot(xlim=[-10, lim], logy=True, ax=ax, ylim=[1000000, 0.5], lw=1)
....: ax.set_yscale('log', base=10)
....: ticks = [1,10,100,1000,10000, 1e5, 1e6]
....: ax.set_yticks(ticks)
....: ax.set_yticklabels([f'{x:.0g}' for x in ticks]);
....: ax.set(title=title);
....: if ax is axs[0]: ax.set(ylabel='Return period (years)');
....:
In [27]: fig.suptitle("Feller's approximation to aggregate PMLs")
Out[27]: Text(0.5, 0.98, "Feller's approximation to aggregate PMLs")
Note
These graphs demonstrate computational facility. I’m not suggesting one million year PML is a reliable estimate. But the figure is reliably computing what the specified statistical model implies. The losses shown range up to USD 5 trillion, about 20% of GDP.
2.8.3. Jewson’s US Wind Climate Change Estimates
Jewson Table 2 provides estimates for the impact of a 2 degree Celcius increase in global mean surface temperature (GMST) on event frequency by Safir-Simpson category. He also provides the standard deviation of the impact. These are added in the next dataframe.
In [28]: jewson['Freq Chg'] = [None, 1.011, 1.095, 1.134, 1.179, 1.236]
In [29]: jewson['Freq Chg SD'] = [None, 0.3179, .4176, .4638, .5174, .5830]
In [30]: qd(jewson.loc["1":"5", ['Freq Chg', 'Freq Chg SD']])
Freq Chg Freq Chg SD
Cat
1 1.011 0.3179
2 1.095 0.4176
3 1.134 0.4638
4 1.179 0.5174
5 1.236 0.583
He models the impact of climate change on PMLs by assuming the frequency of each category is perturbed using a lognormal with mean and standard deviation given by the last two columns of the above table. He assumes that the perturbations across categories are comonotonic. In actuarial terms, he is using comonotonic frequency mixing variables, to create a mixed compound Poisson.
We can create a similar effect using aggregate first by adjusting the baseline event frequencies by the Freq Chg column and then by applying shared mixing across all events together (resulting in comonotonic perturbations). We select a mix CV equal to Jewson’s estimate for Category 4 events. The categories are similar — in light of the overall uncertainty of the analysis.
In [31]: qd((jewson.iloc[1:, -1] / jewson.iloc[1:, -2]))
Cat
1 0.31444
2 0.38137
3 0.40899
4 0.43885
5 0.47168
In [32]: mix_cv = 0.5174 / 1.179
In [33]: mix_cv
Out[33]: 0.4388464800678541
The adjusted model is built using inverse Gaussian mixing variables (slightly thicker tail than gamma), rather than Jewson’s lognormals. Note that the standard deviations increase but the CVs decrease.
In [34]: wcc = build('agg Cat:USWind:Wcc '
....: f'{jewson.loc["1":"5", "EN"].to_numpy() * jewson.loc["1":"5", "Freq Chg"].to_numpy()} claims '
....: f'sev lognorm {jewson.loc["1":"5", "ES_W"].to_numpy()} '
....: f'cv {jewson.loc["1":"5", "CV_W"].to_numpy()}'
....: f'mixed ig {mix_cv}'
....: , bs=1/8)
....:
In [35]: mcc = build('agg Cat:USWind:Mcc '
....: f'{jewson.loc["1":"5", "EN"].to_numpy() * jewson.loc["1":"5", "Freq Chg"].to_numpy()} claims '
....: f'sev lognorm {jewson.loc["1":"5", "ES_M"].to_numpy()} '
....: f'cv {jewson.loc["1":"5", "CV_M"].to_numpy()}'
....: f'mixed ig {mix_cv}'
....: , bs=1/8 )
....:
In [36]: qd(wcc)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1.7954 0.86578 1.1454
Sev 10.646 10.645 -4.0062e-05 2.4249 2.4246 9.5485 9.3409
Agg 19.113 19.112 -4.0062e-05 2.0062 2.006 6.2354 6.1214
log2 = 16, bandwidth = 1/8, validation: fails sev skew, agg skew.
In [37]: mv(wcc)
mean = 19.1129
variance = 1470.276
std dev = 38.3442
In [38]: qd(mcc)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1.7954 0.86578 1.1454
Sev 16.95 16.949 -4.9432e-05 2.9533 2.9491 14.535 13.891
Agg 30.432 30.431 -4.9442e-05 2.368 2.3651 9.6255 9.2342
log2 = 16, bandwidth = 1/8, validation: fails sev cv, agg cv.
In [39]: mv(mcc)
mean = 30.4321
variance = 5193.212
std dev = 72.0639
The new models produce the following AALs, compare Jewson Figure 3.
In [40]: base = pd.concat((w.report_df.loc['agg_m'].T,
....: m.report_df.loc['agg_m'].T), axis=1,
....: keys=['Weinkle', 'Martinez']); \
....: cc = pd.concat((wcc.report_df.loc['agg_m'].T,
....: mcc.report_df.loc['agg_m'].T), axis=1,
....: keys=['Weinkle', 'Martinez']); \
....: df = pd.concat((base, cc), axis=1,
....: keys=[' Base', 'Adjusted']); \
....: df[('Change', 'Martinez')] = (df[('Adjusted', 'Martinez')] - df[(' Base', 'Martinez')]); \
....: df[('Change', 'Weinkle')] = (df[('Adjusted', 'Weinkle')] - df[(' Base', 'Weinkle')]); \
....: df[('Pct Change', 'Martinez')] = (df[('Adjusted', 'Martinez')] -
....: df[(' Base', 'Martinez')]) / jewson.iloc[0]['EX_M']; \
....: df[('Pct Change', 'Weinkle')] = (df[('Adjusted', 'Weinkle')] -
....: df[(' Base', 'Weinkle')]) / jewson.iloc[0]['EX_W']; \
....: df = df.iloc[[0,1,2,3,4,6]]; \
....: df.index = [1,2,3,4,5, 'Total']; \
....: df.index.name = 'Category'; \
....: df = df.swaplevel(axis=1); \
....: df = df.sort_index(axis=1, ascending=[False, True]); \
....: qd(df.stack(0).swaplevel(0).sort_index(ascending=[False,True]))
....:
Base Adjusted Change Pct Change
Category
Weinkle 1 1.6188 1.6366 0.017807 0.0010663
2 1.784 1.9535 0.16948 0.010149
3 4.68 5.3071 0.62712 0.037552
4 7.446 8.7788 1.3328 0.07981
5 1.1625 1.4368 0.27435 0.016428
Total 16.691 19.113 2.4216 0.14501
Martinez 1 2.1016 2.1247 0.023118 0.00087236
2 2.556 2.7988 0.24282 0.009163
3 6.444 7.3075 0.8635 0.032585
4 13.991 16.495 2.5044 0.094505
5 1.38 1.7057 0.32568 0.01229
Total 26.473 30.432 3.9595 0.14942
Here are plots of the base and adjusted AEP and OEP curves. Compare Jewson Figure 5 (a) and (b) for aggregate and Figure 6 (a) and (b) for occurrence.
In [41]: fig, axs = plt.subplots(2, 2, figsize=(2 * 3.5, 2 * 2.5), constrained_layout=True)
In [42]: axs = axs.flat[::-1]
In [43]: for axo, axa, (mw, mwcc), title in zip(axs.flat[0::2], axs.flat[1::2], [(m, mcc), (w, wcc)], ['Martinez', 'Weinkle']):
....: bit = 1 / ((1 - mw.density_df.F_sev.loc[:2000]) * mw.n)
....: bit.plot(logy=True, ax=axo, label='OEP');
....: bit = 1 / ((1 - mwcc.density_df.F_sev.loc[:2000]) * mw.n)
....: bit.plot(logy=True, ax=axo, label='OEP, climate chanage');
....: bit = 1 / (1 - mw.density_df.F.loc[:2000])
....: bit.plot(logy=True, ax=axa, label='AEP');
....: bit = 1 / (1 - mwcc.density_df.F.loc[:2000])
....: bit.plot(logy=True, ax=axa, label='AEP, climate change');
....: axo.set(title=f'{title} OEP');
....: axa.set(title=f'{title} AEP');
....:
In [44]: for ax in axs.flat:
....: ax.set(xlim=[-10, 325], ylim=[130, .5], xlabel='Loss');
....: if ax in [axs.flat[1], axs.flat[3]]:
....: ax.set(ylabel='Return period');
....: ax.set_yscale('log', base=10);
....: ticks = [1, 2, 5, 10, 20, 50, 100]
....: ax.set_yticks(ticks);
....: ax.set_yticklabels([f'{x}' for x in ticks]);
....: ax.legend()
....:
In [45]: fig.suptitle('Impact of climate change on AEP and OEP curves');
2.8.4. ILW Pricing
Industry Loss Warranties (ILW) are securities that pay an agreed amount if losses from a named peril exceed a threshold during the contract term. They are usually written on an occurrence basis and are triggered by losses from a single event. For example, a US hurricane $20 billion ILW pays 1 if there is a US hurricane causing $20 billion or more losses during the contract period. They are used by insurers to provide cat capacity. Because they are not written on an indemnity basis there is no underwriting, which simplifies their pricing.
Brokers publish price sheets for ILWs to give a view of market pricing. Price is expressed as a percentage of the face value. A recent sheet quoted prices for US hurricane as follows.
The next dataframe adds expected losses and compares them to the ILW pricing. The expected loss is given by the occurrence survival function — it is simply the probability of attaching the layer. The EL columns show Jewson’s expected losses across the four views discussed above. The impact on EL is only caused by greater event frequency. Its effect increases with attachment.
In [46]: views = ['Weinkle', 'Weinkle Adj', 'Martinez', 'Martinez Adj']
In [47]: ilw = pd.concat((x.density_df.loc[[15, 20, 25, 30, 40, 50, 60],
....: ['S_sev']].rename(columns={'S_sev': 'EL'})
....: for x in [w, wcc, m, mcc]),
....: axis=1, keys=views,
....: names=['View', 'Stat']); \
....: ilw['Price'] = [.47, .38, .33, .275, .175, .13, .1075]; \
....: ilw.index.name = 'Trigger'; \
....: ilw = ilw.iloc[:, [-1,0,1,2,3]]; \
....: qd(ilw, float_format=lambda x: f'{x:.4f}')
....:
View Price Weinkle Weinkle Adj Martinez Martinez Adj
Stat EL EL EL EL
Trigger
15.0 0.4700 0.1622 0.1733 0.2075 0.2207
20.0 0.3800 0.1261 0.1353 0.1666 0.1782
25.0 0.3300 0.1013 0.1091 0.1390 0.1492
30.0 0.2750 0.0832 0.0898 0.1188 0.1279
40.0 0.1750 0.0586 0.0634 0.0910 0.0983
50.0 0.1300 0.0431 0.0467 0.0725 0.0785
60.0 0.1075 0.0326 0.0354 0.0593 0.0643
Cat pricing is often expressed in terms of the implied multiple: the ratio of premium to EL (reciprocal of the loss ratio).
Cat pricing multiples are usually in the range of 2 to 5 times the standard commercial models’ estimates of expected loss. The base model pricing multiples, shown below, fall into this range. The climate adjusted multiples fall outside, which is not unexpected, since the pricing is for the coming period and not a future climate-impacted period.
In [48]: ilw[[(v, 'Multiple') for v in views]] = ilw[['Price']].values / ilw[[(v, 'EL') for v in views]]; \
....: qd(ilw.iloc[:, [0,5,6,7,8]])
....:
View Price Weinkle Weinkle Adj Martinez Martinez Adj
Stat Multiple Multiple Multiple Multiple
Trigger
15.0 0.47 2.8979 2.7117 2.2652 2.1299
20.0 0.38 3.0142 2.8079 2.2803 2.1327
25.0 0.33 3.2578 3.0257 2.3739 2.2117
30.0 0.275 3.3058 3.0637 2.3144 2.1502
40.0 0.175 2.9845 2.7581 1.9239 1.7806
50.0 0.13 3.0188 2.7849 1.7938 1.6562
60.0 0.1075 3.3001 3.0409 1.8132 1.6714
The next table shows implied distortion parameters calibrated to market pricing for the dual
and proportional hazard (PH)
parametric families (see Mildenhall and Major [2022]). In both cases, a higher parameter corresponds to a higher risk load. The dual is body-risk centric and the PH is tail-risk centric. The indicated parameters are quite high, consistent with the expense of bearing cat risk. (The parameters are incomparable between distortions.)
In [49]: params = pd.concat((np.log(1 - ilw[['Price']]).values / np.log(1 - ilw.xs('EL', axis=1, level=1)),
....: np.log(ilw[['Price']].values) / np.log(ilw.xs('EL', axis=1, level=1))),
....: axis=1, keys=['Dual', 'PH'])
....:
In [50]: qd(params.xs('Dual', axis=1, level=0))
View Weinkle Weinkle Adj Martinez Martinez Adj
Trigger
15.0 3.5877 3.3355 2.7301 2.5465
20.0 3.5474 3.2875 2.6224 2.436
25.0 3.7497 3.4678 2.6757 2.4785
30.0 3.7027 3.4194 2.5422 2.3499
40.0 3.1837 2.9347 2.0172 1.8596
50.0 3.1637 2.9131 1.8511 1.7036
60.0 3.4341 3.1598 1.8608 1.7107
In [51]: qd(params.xs('PH', axis=1, level=0))
View Weinkle Weinkle Adj Martinez Martinez Adj
Trigger
15.0 0.41507 0.4308 0.48009 0.49965
20.0 0.46722 0.48378 0.53997 0.56093
25.0 0.4842 0.50034 0.56186 0.58276
30.0 0.51916 0.53554 0.60606 0.62775
40.0 0.6145 0.63208 0.72704 0.7513
50.0 0.6487 0.66578 0.77736 0.80174
60.0 0.65132 0.66726 0.78937 0.8128
2.8.5. Secondary Uncertainty
Secondary uncertainty is the practice of expanding cat model simulated output by assuming that the results from each event form a distribution. It is usual to assume the distribution is a beta. The model output provides the beta’s mean and standard deviation. Given this output, modelers often need to compute statistics, such as a layer expected loss, reflecting the secondary uncertainty. This calculation can be performed in aggergate as follows.
Assumptions: Assume one location with a TIV of 2500 and simple cat model output with only three equally-likely events with mean losses 100, 200, and 1100 and secondary uncertainty standard deviation 100, 150, and 600. The overall event frequency is 1.6 with a Poisson distribution.
Question: What is the expected loss to a 1000 xs 1000 per risk cover with and without secondary uncertainty?
Solution: Start by building the answer without secondary uncertainty. It is convenient to put the assumptions in a dataframe.
In [52]: df = pd.DataFrame({'GroundUpLoss': [100, 200, 1100],
....: 'GroundUpSD': [100, 150, 600]})
....:
The model with no secondary uncertainty is a simple mixed severity.
In [53]: base = build('agg Cat:Base '
....: '1.6 claims '
....: f'dsev {df.GroundUpLoss.values} '
....: 'occurrence ceded to 1000 xs 1000 '
....: 'poisson'
....: , bs=1)
....:
In [54]: qd(base)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1.6 0.79057 0.79057
Sev 466.67 33.333 -0.92857 0.96362 1.4142 0.68097 0.70711
Agg 746.67 53.333 -0.92857 1.0979 1.3693 1.2973 1.3693
log2 = 15, bandwidth = 1, validation: n/a, reinsurance.
To incorporate the secondary uncertainty, we first compute the beta parameters using the method of moments. Then build the Aggregate model incorporating secondary uncertainty in each loss.
In [55]: tiv = 2500; \
....: m = df['GroundUpLoss'] / tiv; \
....: v = (df['GroundUpSD'] / tiv) ** 2; \
....: sev_a = m * (m * (1 - m) / v - 1); \
....: sev_b = (1 - m) * (m * (1 - m) / v - 1); \
....: sec = build(f'agg Cat:Secondary '
....: '1.6 claims '
....: f'sev {tiv} * beta {sev_a.values} {sev_b.values} wts=3 '
....: 'occurrence ceded to 1000 xs 1000 '
....: 'poisson')
....:
In [56]: qd(sec)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1.6 0.79057 0.79057
Sev 466.67 96.384 -0.79346 1.2367 2.5758 1.5399 2.6259
Agg 746.67 154.21 -0.79346 1.2573 2.1844 1.6706 2.4651
log2 = 16, bandwidth = 1/2, validation: n/a, reinsurance.
Including secondary uncertainty nearly triples the expected loss to the layer, from 53 to 154. Had the third loss been only 1000, there would be no loss at all to the layer without secondary uncertainty.
The next plot compares the gross severity and aggregate distributions.
In [57]: fig, axs = plt.subplots(1, 2, figsize=(2 * 3.5, 2.45), constrained_layout=True); \
....: ax0, ax1 = axs.flat
....:
In [58]: base.reinsurance_df['p_sev_gross'].cumsum().plot(xlim=[0, 2500], ax=ax0, label='Base'); \
....: sec.reinsurance_df['p_sev_gross'].cumsum().plot(xlim=[0, 2500], ax=ax0, label='Secondary'); \
....: base.reinsurance_df['p_agg_gross_occ'].cumsum().plot(xlim=[0, 2500], ax=ax1, label='Base'); \
....: sec.reinsurance_df['p_agg_gross_occ'].cumsum().plot(xlim=[0, 2500], ax=ax1, label='Secondary'); \
....: ax0.set(title='Occurrence', xlabel='Loss', ylabel='Distribution'); \
....: ax1.set(title='Aggregate', xlabel='Loss', ylabel='Distribution'); \
....: ax0.legend();
....:
In [59]: ax1.legend();
2.8.6. Summary of Objects Created by DecL
The following objects are created by build() in this guide.
In [60]: from aggregate import pprint_ex
In [61]: for n, r in build.qlist('^Cat:').iterrows():
....: pprint_ex(r.program, split=20)
....: