About the return period of a catastrophe

To understand catastrophes like earthquakes stochastically, When a natural hazard event like an earthquake affects 5 a region and generates a natural catastrophe (NatCat), the following questions arise: How often does such an event occurs? What is their return period (RP) should be quantified for the concerned region. Measures such as event magnitudes or indexes are less helpful for this purpose.)? We derive the combined return period (CRP) from a concept of extreme value statistics and theory the pseudo-polar coordinates of extreme value theory. The. A CRP is the (weighted) meanaverage of the local RPs and is again an RP; other metrics do not provide such testable reproductivity. We demonstrate CRP’s opportunities on 10 extratropical cyclones (winter storms) over Germany, including validation and bias correction of local RP estimates. Furthermore, we introduce new estimation methods for the RP of an event loss (risk curve) via CRP-scaling of historical storm fields. For high RP, the resulting event losses of the German insurance market are higher in the case of max-stable dependence. The latter means the same dependence level between local maxima of a year as of a decade. However, intensities. Since CRP’s reciprocal is its expected exceedance frequency, which applies to any RP per stochastic definition, the concept is testable. As 15 we show, the CRP is related to the spatial dependence is not stable but decreases by increasing period. Such control ofcharacteristic of the NatCat generating hazard event and their spatial dependence is not realized byof corresponding local block maxima (e.g., annual wind speed maximum). For this purpose, we extend previous risk modelsconstruction for maxstable random fields from science and industry. Our loss estimates for RP of 200 years areextreme value theory and consider a recent concept from NatCat research. Based on the CRP, we also significantly smaller than those of European regulation's 20 standard modeldevelop a new method to estimate the NatCat risk of a region via stochastic scaling of historical fields of local event intensities (represented by records of measuring stations) and averaging corresponding risk parameters such as the event loss with a defined RP. Our application example is winter storm (extratropical cyclones) over Germany. We analyze wind station data and estimate local hazard, CRP of historical events and the risk curve of insured event losses. The most destructive storm of our observation 25 period of 20 years is Kyrill 2002 with weighted CRP 16.97±1.75. The CRPs could be successfully tested statistically. We also state that our risk estimate is higher for the max-stable case than for the non-max-stable. Max-stable means that the dependence measure (e.g., Kendall’s ) for annual wind speed maxima of two wind stations has the same value as for maxima of higher block size such 10 or 100 years since the copula (the dependence structure) remains the same. However, the spatial dependence decreases with increasing block size; a new statistical indicator confirms this. Such control of spatial characteristic and 30 dependence is not realized by the previous risk models in science and industry. We compare our risk estimates to these.

Furthermore, we introduce new estimation methods for the RP of an event loss (risk curve) via CRP-scaling of historical storm fields. For high RP, the resulting event losses of the German insurance market are higher in the case of max-stable dependence.
The latter means the same dependence level between local maxima of a year as of a decade. However, intensities. Since CRP's reciprocal is its expected exceedance frequency, which applies to any RP per stochastic definition, the concept is testable. As 15 we show, the CRP is related to the spatial dependence is not stable but decreases by increasing period. Such control ofcharacteristic of the NatCat generating hazard event and their spatial dependence is not realized byof corresponding local block maxima (e.g., annual wind speed maximum). For this purpose, we extend previous risk modelsconstruction for maxstable random fields from science and industry. Our loss estimates for RP of 200 years areextreme value theory and consider a recent concept from NatCat research. Based on the CRP, we also significantly smaller than those of European regulation's 20 standard modeldevelop a new method to estimate the NatCat risk of a region via stochastic scaling of historical fields of local event intensities (represented by records of measuring stations) and averaging corresponding risk parameters such as the event loss with a defined RP.
Our application example is winter storm (extratropical cyclones) over Germany. We analyze wind station data and estimate local hazard, CRP of historical events and the risk curve of insured event losses. The most destructive storm of our observation 25 period of 20 years is Kyrill 2002 with weighted CRP 16.97±1.75. The CRPs could be successfully tested statistically. We also state that our risk estimate is higher for the max-stable case than for the non-max-stable. Max-stable means that the dependence measure (e.g., Kendall's ) for annual wind speed maxima of two wind stations has the same value as for maxima of higher block size such 10 or 100 years since the copula (the dependence structure) remains the same. However, the spatial dependence decreases with increasing block size; a new statistical indicator confirms this. Such control of spatial characteristic and 30 dependence is not realized by the previous risk models in science and industry. We compare our risk estimates to these.

Introduction
WhenAfter a natural catastrophe (NatCat)hazard event such as a large windstorm or an earthquake has occurred, in a defined region (e.g., in a country) and results in a natural catastrophe (NatCat), the question arises, how oftendoes such random eventsevent appear, i.e.,? What is the corresponding return period (RP, also called recurrence interval)? Before discussing this 35 issue, we underline that the extension of river floodsflood events or windstorms in time and space depend on the scientific and socio-economic event definition. These areThe definitions may vary by peril and is not our topic even though they influence our research objectthe RP of a hazard and NatCat event.
The RP of an event magnitude or index is so far frequently used as a stochastic measure of a catastrophe. For example, there are different magnitudes scales for earthquakes (Bormann and Saul, 2009). Unfortunately, corresponding RPs for the 40 source regionBut their RP may not correspond well with the catastrophiclocal consequences since the hypocentrehypocenter position also determines destructivenesslocal event intensities and effects. For floods, regional or global magnitude scales are not in use (Guse et al., 2020). For hurricanes, the Saffir-Simpson scale (National Hurricane Centre, 2020) is a magnitude measure; but the random storm track also influences the destruction's extent of the destruction. Extratropical cyclones hitting Europe, called winter storms, are measured by a storm severity index (SSI; Roberts et al., 2014) or extreme wind index (EWI; 45 Della-Marta et al., 2009). Their different definitions result in quite different RP for the same events. Also, their statistical models include assumptions and pitfalls. In rare scientific publications about risk modelingmodelling for the insurance industry, such as by MitechellMitchell-Wallace et al. (2017), better and universal approaches for the RP are not offered. In sum, previous approaches are not very fruitful. Nevertheless, RP is needed to understand and manage NatCat risk by decisionmakers in administration, policy, and industry and research the consequences of climate change , 50 Schwierz et al. 2010). In the end, the RP of losses and damage (the risk curve) is neededsuccessful regarding the stochastic quantification of a hazard or NatCat event what our motivation is to develop a new approach. We mathematically derive the concept of combined return period (CRP) being the average of RPs of local event intensities from an approach of extreme value theory and statistics. As we will show by a combination of previous and new approaches from stochastic and NatCat research, the concept of CRP is strongly related to the spatial association/dependence between the local event intensities, their 55 RPs, and corresponding block maxima such as annual maxima.
This spatial dependence is less or inappropriately considered in previous research about NatCat. The issue is only a marginal topic in the book about NatCat modelling for insurance industry by Mitchell-Wallace et al. (2017, Section 5.4.2.5). Jongman's et al. (2014) model for European flood risk considers such dependence explicitly. However, their assumptions and estimates are not appropriate according to Raschke (2015). In statistical journals, max-stable dependence models have been 60 applied to natural hazards without a systematic test of the stability assumption, such as the snow model by Blanchet and Davison (2011) for Switzerland and the river flood model by Asadi et al. (2015) for the Upper Danube River system. Maxstable dependence means that the copula (the dependence structure of a bi-or multivariate distribution) and corresponding value of dependence measures are the same for annual maxima as for ten-year maxima or these of a century (Dey et al., 2016). 4 trend in time series of a wind station in Potsdam (Germany) over several decades shows. Wichura (2009) assumes changed local roughness condition over the time as reason, Mudelsee (2020) the climate change.
The research of spatial dependence of natural hazards is not an end in itself, the final goal is an answer to the question about the NatCat risk. What is the RP of events with aggregate damage or losses in a region equal or higher to a defined level? By using CRP, we quantify the risk via stochastic scaling of fields of local intensity (e.g., river discharge or wind gust peak) 100 at a geographical site.intensities of historical events and averaging corresponding risk measures. This local intensity occurs (approximately) as point event of a Poisson process with points ; the number of independent new approach significantly extends the methods to calculate a NatCat risk curve. Previous opportunities for a risk estimate are the conventional statistical models that are fitted to observed or re-analyzed aggregated losses (also called as-if losses; ) of historical events as used by Donat et al. (2011) and by Pfeifer (2010) for annual sums. The advantages of such simple models are the controlled stochastic 105 assumptions and the small number of parameters; the disadvantages are high uncertainty for widely extrapolated values and limited opportunities to consider further knowledge. The NatCat models in (re)insurance industry combine different components/sub-models for hazard, exposure (building stock or insured portfolio) and corresponding vulnerability (Mitchell-Wallace et al. 2017, Section 1.8;Raschke, 2018) and offers better opportunities for knowledge transfer such as the differentiated projection of a market model on a single insurer. However, the corresponding standard error of the risk estimates 110 is frequently not quantified (and cannot be quantified). The numerical burden of such complex models is high. Tens of thousands of NatCat events must be simulated (Mitchell-Wallace eta al., 2017, Chapter 1). The question arises, what is the stochastic criterion for the simulation of a reasonable event set in NatCat modelling? As far as we know, scientific NatCat models for European winter storms (extratropical cyclones) are based on numerical simulations Schwierz et al., 2010;Osinski et al., with ≥ during a defined unit period (e.g., a season or a year) is a random integer 115 number that follows a Poisson 2016) and are not intensively validated regarding spatial dependence.
To answer our questions, we start with topics of extreme value statistics in the 2 nd Section and illuminate max-stability in the univariate sense, for the dependence structure (copula) of the bivariate case and max-stable random fields. We also extend Schlather's (2002) 1 st theorem with focus on spatial dependence. The more recent approaches of hazard event related area functions (Raschke, 2013) and survival functions (Jung and Schindler, 2019) of local event intensities within a region are 120 implemented therein to characterize spatiality. In the 3 rd Section, we derive the CRP from the concept of pseudo polar coordinates of extreme value statistics and explain its testability and scaling opportunity and corresponding risk estimate.
Subsequently, in Section 4, we apply the new approaches to winter storms (extratropical cyclones) over Germany to demonstrate their potentials. This application implies several elements of conventional statistics which are explained in Section 5. Finally, we summarize and discuss our results and give an outlook in Section 6. Some stochastic and statistical details are 125 presented in the Supplementary and Supplementary Data to remain clarity of the main paper and limit its extent. In the entire paper, we must consider several stochastic relations. Therefore, the same mathematical symbol can have different meanings in different subSections. We also expect that the reader is more familiar with statistics and stochastic than only with basics 5 about random variables. Statistical significance, goodness-of-fit tests, random fields, or a Poisson (point) process (Upton and Cook, 2008) should be familiar terms. 130 2 Max-stability in statistics and stochastic

The univariate case
Before we formulate the CRP and discuss their opportunities, we must present, discuss, and extend a corresponding topicmax-stability in extreme value statistics especial of random process and fields. Max-stability has its origin in univariate statistics. The cumulative distribution functions (CDF) ( ) of maximum = ( 1 , … , ) of identical and 135 independently distributed (iid) random variables with CDF ( ) (for the non-exceedance probability ( ≤ )) is A CDF ( ) is max stable if the linear transformed maximum (with parameters and ) has the same distribution , Def. 3.1) The Fréchet distribution (Beirlant et al., 2004;Falk et al. , Tab. 2.1) is such a max-stable distribution, also called extreme value distribution, with CDF For the unit Fréchet distribution is = 1 and the transformation parameters are = 0 and = . The most distribution types are not max-stable, but their distribution of maxima (1) converges to an extreme value distribution by increasing sample 145 size , called the block size in this context (Beirlant et al., 2004, Chapter 3). These are well-known facts, and we can only refer to some of a very high number of corresponding publications (e.g., de Haan and Ferrira, 2007;Falk et al. 2011) with an ).  gives a good overview for practitioners.

Max-stable copulas
It is also well-known that a bivariate CDF ( , ) can replace by a copula ( , ) and the marginal CDFs ( ) and ( ) 150 The copula approach represents a universal distinction between the marginal distributions and the dependence structure and was introduced by Sklar (1959). As there are different univariate distributions (types) there are different copulas (types). Mari and Kotz (2001) presents a good overview about copulas, their construction principals, and different views on dependence. Max-stability is also a property of some copulas, called max-stable copula or extreme copula. A max-stable copula 155 remains the same for pairs of component wise maxima ( , ) as it was already for the underlying pairs ( , ); the copula parameters including dependence measure such as Kendall's rank correlation are equal. The formal definition is (Dey et al., 2016, (2.3)) ( , ) = ( 1/ , 1/ )

Max-stability of stochastic processes 160
The spatial extension of the bivariate situation and corresponding distribution is the random field ( ) at points in the space ℝ with dimensions (e.g., Schlather, 2001). In our application, ℝ 2 is the geographical space and is the corresponding coordinate vector. At one point/site in ℝ , ( ) is the marginal distribution of the local random variable . There are various differentiations and variants such as (non)stationarity or (non)homogeneity. A max-stable random field has max-stable marginal distributions and the copulas between to margins are also max-stable. Schlather (2002) has formulated and proofed 165 a construction of a max stable random field (we cite his 1 st theorem with the same notation) is a stationary max-stable process with unit Fréchet margins. 170 Extreme value statistics is interested in the max-stable dependence structure (copula) between the margins, the unit Fréchet distributed random variables at fixed points in space ℝ . From perspective of NatCat modelling in the geographical space ℝ 2 and with ( ) ≥ 0, the entire generating process is interesting. The Poisson (point) process represents all hazard events (e.g., storms) of a unit period such as a hazard season or a year and has two parts, and . The point events on (0, ∞) are a stochastic event magnitude and scale the field of local events intensity ( ), short called intensity 175 which represents all point events ( ) at sites . The random coordinate is a kind of epicenter in the meaning of NatCat with the (tendentiously) highest local event intensity such as maximum wind speed, maximum hail stone diameter, or peaks of earthquake ground accelerations. The copied random function ( ) determines the pattern of a single random event in the space ℝ . ( ) or its local expectation (converges to 0 or is 0 if magnitude ‖ ‖ of coordinate vector converges to infinity due 180 to the measurability condition in Theorem 1. This also applies to NatCat events with limited geographical extend. Schlather (2002) has demonstrated the flexibility of his construction by presenting realizations of maximum fields for different variants of ( ) . Its measurability condition is fulfilled by classical probability density functions (PDF, first derivative of the CDF, Coles, 2001, Section 2.2) of random variables. For instance, Smith (1990, an unpublished and frequently cited paper) used the PDF of the normal distribution. We present in the Supplementary, Section 4, some examples of the 185 random function ( ) to illustrate the universality of the approach. ( ) can also imply random parameters such as verities of standard deviation of applied PDF or is combined with a random field. 7 Both, and ( ) with fixed , are point events of Poisson processes with intensity −2 . This is the expected value) [… ] -point density and determines the exceedance frequency (EF) function (with . The latter is the expected number of point events ( ) > and > 190 The entire construction of Theorem 1 is also a kind of shot noise field according to the definitions of Dombry (2012); and Schlather (2002) has also published a construction of max-stable random field without a random function but with a stationary random field. The logarithmic variant of Theorem 1 (logarithm of (6,7)) also results in a max-stable random field, however, the marginal maxima are unit Gumbel distributed and (8) would be an exponential function. The Brown-Resnick process -195 well-known in stochastic (e.g., Engelke et al., 2011) -determines a max-stable random field with such unit Gumbel distributions and use a nonstationary random field. It is implicitly a construction according to Theorem 1 since for exponential transformation (inverse of logarithmic transformation) the nonstationary random field is the random function of Theorem 1.
The origin of a Brown-Resnick process in ℝ can be fixed but can also be a random coordinate as is in Theorem 1.
The construction of Theorem 1 is already used to model natural hazards in the geographical space. Smith (1990) has applied 200 the bivariate normal distribution as ( ) in a rainstorm modelling. The Brown-Resnick Process has been already applied to river flood (Asadi et al., 2015). Blanchet and Davison (2011) have applied a max-stable model for snow and Raschke et al. (2011) for winter storm, both in Switzerland. And there are similarities to conventional hazard models. Punge's et al. (2014) hail simulation includes maximum hail stone diameter that acts like ln ( ) in (6,7). Raschke (2013) already stated similarity between earthquake ground motion models and Schlather's construction. However, the earthquake magnitude can have a wider 205 influence on the geographical event pattern than a simple scaling. This was one of the motivations to extend and generalize the Schlather's construction (7)  As we show in the Supplementary, Section 2, the marginal Poisson processes ( ) in (9) have the same exceedance frequency (8) as (7). Correspondingly, ( ) in (10) is also unit Fréchet distributed as in (6). Schlather's construction is a special case of (9,10) with = 0; (9,10) only implies max stability of spatial dependence in this case, what we discuss in the following Section.

Spatial characteristics and dependence 215
We now illustrate spatial max-stability and its absence by examples of (9, 10) with standard normal PDF as random function ( ) in a one-dimensional parameter space ℝ =1 . For this purpose, we apply the simulation approach of Schlather (2002) and 8 generate random events within a range (-10,10) for local event intensities within the region/range (-4,4) in ℝ 1 by a Monte Carlo simulation. According to Schlather's procedure that processes a series of random numbers from a (pseudo) random generator, only the events for the large are simulated which implies incompleteness for smaller events. This does not 220 significantly affect the simulated field ( ) of maxima. However, we can only consider this simulation for ≥ 0 in (9,10) since the edge effects increase for increasing if < 0. In Figure 1 a, we show fields for one realization of Schlater's theorem ( = 1, equivalent to one year or one season in NatCat modelling) for the max-stable case with = 0 in (9,10). With the same series of random numbers, we generate fields of = 100 realizations of in Figure 1b. It has the same pattern = 1 and is the same when we linear transform the local intensities with division by = 100. The entire generating processes 225 are max-stable as the resulting marginals and dependence and association/dependence between marginals are. In contrast to this total max-stability, the example with = 0.1 results in different pattern for = 1 and = 100 in Figure 1 c and d. The shape of the event fields gets sharper for larger , only the marginals are max-stable, not their spatial relations.
To illustrates the effect on spatial dependence quantitatively, we have generated local maxima ( ) from (10) by Monte Carlo simulation with 100,000 repetitions and computed corresponding dependence measure Kendall's (Kendall, 1938;Mari 230 and Kotz  Section 6.2.6). As depicted in Figure 2 a and b, the functions are the same if = 0 and differs if = 0.1, the dependence is decreasing by increasing if > 0. In Figure 2c, the functions are shown for the limit cases full dependence with the same value of ( ) at each point and full independence with ( ) = 0 everywhere except one point.

240
Beside our extension of Schlather's theorem, we also consider a more recent approach from NatCat research to understand the spatial characteristic. Raschke (2013) described an earthquake event by its area function for the peak ground accelerations.
This is a cumulative function and measures the set of points in the geographical space (the area) with an event intensity higher than the argument of the function. The area function is limited here to a region and is normalized ( and symbolises the region's bounds, the integral in the denominator is the area of the region in ℝ 2 , is an indicator function ) 245 . (1) and is now like a survival function of a random variable (decreasing with value of functions between 0 and 1) which describes the exceedance probability in contrast to a CDF for non-exceedance probability (Upton and Cook, 2008). Jung and Schindler (2019) have already applied such aggregating functions to German winter storm events and call them explicitly survival function. However, not every normalized aggregating decreasing function is based on an actual random variable. And survival functions are not used in statistics to describe regions of random fields or random function as far as we know. Nonetheless, we 255 use the area function (11) to characterize and research the spatiality of the event field ( ) in a defined region. As an example, the area function for the strongest events in Figure 1 is shown in Figure 3 a. The differences between the variants = 1 versus = 100 and = 0 versus = 0.1 corresponds with the differences between these events in Figure According to (12), any scaling of by a factor > 0 results in proportional scaling of expectation and standard deviation 265 in (12), and the CV remains constant. Correspondingly, random magnitude in (9,10) only scales the field ( ) in the maxstable case with = 0 and influences the expectation of ( ) but not the CV. Thus, the CV is independent on the expectation. This does not apply to the non-max-stable case with ≠ 0 in (9,10). These different behaviors are detectable for the examples of Figure 1 b and d in Figure 3 c and d. For the max-stable case, the scale/slope parameter of the linearized regression function does not differ significantly from 0 according to the t-test (Fahrmeir et al., 2013, Section 3.3). For max-stable case, the 270 10 regression function is statistically significant with a p-value of 0.00. Linearization is provided by logarithm of CV and expectation/average. For completeness, the full dependence case of Figure 3 b corresponds with an CV 0.
In sum of Section 2, Schlather's 1 st Theorem has parallels to NatCat models, is used already in hazard models and was extended here to the non-max stable case regarding spatial dependence and characteristic. Statistical indication for maxstability is the independence of the spatial dependence measure from the block size (e.g., one versus ten years) and 275 independence between CV and expectation of the area function (11). Otherwise, non-max-stability is indicated.

The stochastic derivation
Let the point event , be the local intensity at site of a hazard event as a member of the set of all events of a defined unit period such as a year, hazard season, or half season. This local event intensity might be the maximum river discharge of a flood, the peak ground acceleration of an earthquake or the maximum wind gust of a windstorm event. The entire number of 285 events with , > during the unit period is = ∑ . is (at least approximately) a Poisson distributed (Upton and Cook, 2008) discrete random variable with an expectation -the expected exceedance frequency, that is the local hazard function in a NatCat model (this is not the hazard function/hazard rate of statistical survival analysis, Upton and Cook, 2008).
This is the bijective frequency function. The and the local hazard curve. Its reciprocal of expected local EF isdetermines the 290 localhazard curve for the RP ( ) = 1/ ( ). As is a point event it's RP ( ) is also a point event and the EFof a point process with frequency function of isaccording to (14) It might be called the RP process. In the Supplementary (Figure 1), the EF function is depicted beside realizations, one of a single process and a second of two associated processes. In extreme value theory and statistics Beirlant et al., 305 2004;Falk et al. 2011), the density of point events at the line of positive real numbers is expressed by the Poisson process's intensity, being the negative first derivative of (3).
Two-point processes can be associated.Since (15) is the same as (8), Schlather's theorem and our extensions directly apply to RP. For completeness, the marginal maxima have a CDF for unit periods (a unit Fréchet distribution for = 1 according to (3)) 310 This is applicable because the probability of non-exceedance for level of the block maxima is the same probability that no events occur with , > which is determined by the Poisson distribution; (6,8) also imply this link.
Schlather's theorem is also based on and implies the concept of pseudo polar coordinates. According to de Haan (1984),) and well explained by  and Falk et al. (2011, Section 8.3.2), two max-stable linked point processes with expected 315 exceedance frequency (15) and point events 1 and 2 are also represented by pseudo polar coordinates with radius and angle = 1 + 2 , = 1 1 + 2 , 1 = , 2 = (1 − ) = 1 1− . (4) As we describe in the Supplementary, Section 1, the expectation of the random element (1 − )/ is 1 and for the conditional expectation of unknown RP 2 with known 1 applies (association is provided) The interest in extreme value theory and statistics , Section 3.8; Beirlant et al., 2004, Section 8.2.3;Falk et al. 2011, Section 4.2) is focused on the distribution function of pseudo angle with cumulative distribution function (CDF) ( ).
It determines the dependence structure, called copula (Marie and Kotz, 2001), of the block maxima of the ). As  write "the angular spread of points of N [the entire point processes -here local maxima of all point events (or ) during 330 considered unit periods (e.g., half years, a season, years, or decades). The univariate CDF of maximum occurred RP of unit periods considers ( ) of (1) and also of (3) If the copula of a bivariate distribution ( 1 , 2 ) of two local maxima] is max-stable, the dependence between annual maxima determined by H, and is the same as between century maxima. It is also called extreme value copula (Marie and Kotz, 2001). The independence gives this max-stability of the dependence structure between pseudoindependent of radial distance [ ]", angle and pseudo radius occurs independently to each other and determine the copula between two marginal maxima ( ) in (4) (Theorem 1. 340 According to . The occurrence of , Section 3.8), the pseudo radius in (17) is once again a point event of a Poisson process with EFfrequency ( ) = 2/ -the double of (315). This means the average of two RP 1 and 2 results in a combined return period (CRP) with EFexceedance frequency function (3). This8,15). We do not have a mathematical proof that (18,19) also applies for the 345 non-max-stable associated point processes. However, max-stable case betweenand non-max-stable cases have the same limits -: full dependence ( 1 = 2 ) and no dependence/full independence ( 1 = 0 if 2 > 0 and vice versa, = = 0 represents the lack of a local event). These are also the limits of non-max-stable dependence structures. Therefore, CRP(19) should also apply forto the non-max-stable dependence, whichcase between these limits. This can be validated heuristically as shownwe demonstrate by an example in the Supplementary. The critical result is that the average of two RPs is again an RP with EF (, 350 Section 3). Such reproductivity does not apply to a metric of the EF or the logarithm of RP. RP has probably been .
More than one RP can be averaged in NatCat studies before, but it was not known that this average remains an RP.since the averaging of two RPs can be done in serial (and the pseudo polar coordinates are also applied to more than two marginal processes). Serial averaging (averaging the last result with a further RP) also implies a weighting; the first considered RPs would be smaller weighted than the last in the final CRP. The general formulation of averaging of RP with weight is 355 A weighting in the averaging is also possible because of the commutative property of numbers The corresponding continuous version via a defined region in a geographical within the region's bounds and in space with If ( ) = ( ) = 1, the integral applies in (21) then the denominator is the area of the region. and CRP is the expectation 370 of the area function (11). This also applies for other weightings if we consider it in the area function, here written for RP ( )

Opportunities and implications
with empirical version for measuring station in the analyzed region .
(23) 375 The weighting, especially the empirical one, can be used in hazard research to compensate an inhomogeneous geographical distribution of measurement stations or a different focus than the covered geographical area such as the inhomogeneous distribution of exposed values or facilities in NatCat research. It has the same effect on the area function as a distortion of the geographical space as used by Papalexiou et al. (2021). Weighted or not, CRP and CV are parameters of the area function.

Testability 380
Before the CRP is applied in stochastic NatCat modelling, it should be tested statistically to validate the appropriateness. A sample of CRPs can be tested by a comparison of its EFexceedance frequency function (315) and their empirical variant.
Therein the empirical EF of the largest CRP in the sample is the reciprocal of the length of the observation period. The 2 nd largest CRP is hence associated to twice the EFexceedance frequency of the largest CRP and so on. It is the same as for empirical exceedance frequency for EQ (e.g., the well-known Gutenberg -Richter relation in Seismology, Gutenberg-Richter, 385 1956). However, not all small events are recorded; the sample is thinned and incomplete in this range. This is. This completeness issue is well known for earthquakes and is here less important if only the distribution (616) of maximum CRPs 14 areis tested. There are a number goodness-of-fit testestests (Stephens, 1986, Section 4.4) for this situation with athe case of known distribution model. The Kolmogorov-Smirnov test is a popular variant.

Besides that, 3.2 The scaling opportunity 390
The CRP also offers the opportunity of stochastic scaling. The CRP and all local RPs in (7) or20) (and ( ) in (8), and by this, the pseudo radius in (4)21)) are scaled viaby a factor . The pseudo angle stays the same.
This means for the pseudo polar coordinates in (17), which applies to the max-stable case, = , = . 395 The pseudo angle is not changed as expected since pseudo radius and pseudo angle are independent in the pseudo polar coordinate for the max-stable case (Section 3.1). This also means that a scaling must be more complex if there is non-max-400 stability. We cannot offer a general scaling method for this situation; however, it must consider/reproduce the pattern of the relation CV versus CRP (example in Figure 3 d) adequately. Irrespective of this, the corresponding event field of local intensities (e.g., maximum wind gust speed) can be computed for the scaled RPs via the inverse of the EF function (1) of local intensity. If the association between the local point processes is max-stable, the value of ( ) is not changed, and nothing is changed in the sampling regarding dependence. This matches the scaling in Schlatter's (2002) two theorems about max-stable 405 random fields (by Schlatter's random element ). Such simple scalability does not apply if we lack max-stability in the spatial dependence; an extensive scaling without special considerations would be critical because of uncontrolled consequences.local RPs via the inverse of the local hazard function ( ) in (14) or ( ) in (13).

Risk estimates by scaling and averaging
The main goal of a NatCat risk analysis is the estimate of a risk curve (Mitchell-Wallace et al., 2017, Section 1), the bijective 410 functional of event loss in a region and corresponding RP, called here event loss return period (ELRP) . As aforementioned, there are two approaches for such estimates with corresponding pros and cons. We introduce an alternative method. According to (18), the expectation of an unknown ELRP is the CRP of the local event intensities; the CRP is an estimate of the ELRP (max-stability between ELRP and CRP provided). To get a good estimate of ELRP, we must average the of many events with the same event loss. We cannot observe such, but we can stochastically scale historical events respectively their 415 local intensity observations. The modelled event loss is the sum of the product of local loss ratio , determined by local event intensity , and local exposure value over all sites (Klawa and Ulbrich, 2003;Della-Marta et al. 2010 with the local vulnerability function , ( , ). To get the event loss for the scaled event, the observed , is replaced by ) (27) 420 with local hazard function ( ) (13) for exceedance frequency and its invers function −1 ( ). The formulation with RP (14) is equivalent. The scaling factor in (27) is the same for all sites/locations respectively as it is in (24) for the CRP and must be adjusted in an iteration until the defined event loss is the result of (26). The scheme in Figure 4 a includes all elements and relations of the scaling approach. Therein, the numerical determination in the scaling scheme has only one direction, from scaled CRP to the event loss. The idea of CRP averaging is also illustrated by Figure 4 b. The standard error of the averaging 425 is the same as for the estimates of an expectation by the sample mean (Upton and Cook, 2008, catchword central limit theorem). According to the well-known Delta method, well explained by Coles (2001, Section 2.6.4), statistical estimates and their 430 standard error can be transferred in other parameter estimation and corresponding standard error by the determined transfer function and its derivates. In its meaning, we can also average the event loss for a fixed/determined CRP respectively its scaled variant. This also applies to the exceedance frequency, the reciprocal of RP, for a determined/fixed event loss.
There is a further chain of thoughts as argument for averaging the exceedance frequency. The scaled intensity fields are like sub-sets in the set of all possible event fields. Each of these subsets implies a relation exceedance frequency to event loss 435 a risk curve. Furthermore, we assume an unknown probability that this sub-set generates a part of the entire risk curve. The latter is for a fixed event loss the aggregation of subset exceedance frequency multiplied with their probability. This is basically the same as the definition of the expectation (12). Therefore, we can apply estimator of the expectation in the estimation of a risk curve and average the exceedance frequencies of risk curves of the subsetsthe scaled event fields. CRP represents the expectation (or its estimate). It can also be derived from the moments of random variables that the coefficient of variation (CV; Upton and Cook, 2008) for (10) is not be concerned about scaling (9) for max-stable situations.
This implies that the CV is independent of CRP in the case of max stability. This underlines that the CRP and area function 450 corresponds with maxima's spatial dependence in the sense of max-stable random fields (Schlather, 2002).

Analysis of 28)
The right side of the equations in (28) implies actual values which can be and are be replaced by estimates. Corresponding uncertainties must be considered in the final error quantification.
We draw attention to the fact that the explained scaling does not change the CV of (23) which implies independence 455 between CRP and CV (Section 2.4). Therefore, the presented scaling only applies to the max-stable case of local hazard. For the non-max-stable case, the scaling factor in (27) must be replaced event wise by which reproduce the observed relation between CRP and CV. An example without max-stability was shown in Figure 3 d.

Overview about data and analysis 460
We have selected the peril winter storms (also called extratropical cyclones or wind winter windstorms) over Germany to demonstrate the application of CRP, we analyze opportunities of the CRP because of good data access and since we are familiar with this peril Raschke, 2015). Our analysis follows the scheme in Figure 4 a, important results are presented in the subsequent Sections, technical details are explained in Section 5. At first, we give an overview.
We analyzed 57 winter storms (also called winter windstorms) over Germany over 20 years, from autumn 1999 to spring 465 2019. Wind station data have been provided by the German meteorological service (DWD, 2020). Different completeness criteria were considered for station selection, and (Supplementary Data, Table 1 and 2) to validate the CRP approach. Different references (Klawa and Ulbrich, 2003;Gesamtverband Deutscher Versicherer, [GDV], 2019; Deutsche Rück, 2020) have been considered for selectingto select the time window per event. The lists are presented in the Supplementary (Excel file). The one and Our definition of a half-power of wind speed maximum are analyzed as local peril intensity. The reasonwinter storm 17 season is explained in Section 3.1 and the appropriateness of the Gumbel distribution for the from September to April of the subsequent year and accept a certain opportunity of contamination of the sample of block maxima of local event intensities and corresponding computation of RP per event with bias correction. To increase the accuracy of estimates by higherextremes from convective windstorm events and a certain opportunity of incompleteness since extratropical cyclones can also be observed outside our season definition (Deutsche Rück, 2020). The term winter storm is only based on the high frequency of 475 extratropical cyclones during the winter. The seasonal maximum is also the annual maximum of this peril.
The maxima per half season (bisected by turn of the year) are analyzed to double the sample size and to increase estimation precision. The appropriateness of this sampling is discussed in Section 5.1. We considered records of wind stations in Germany of DWD (2020; FX_MN003, a daily maximum of wind peaks [m/s], usually wind guest speed) that include minimum record completeness of 90% for analyzed storms, at least 90% completeness for the entire observation period and minimum 55% 480 completeness per half season. Therefore, we analyze the maxima of half-seasons per only consider 141 of 338 DWD wind stations (Supplementary Data, Table 3). We think this is a good balance between large sample size and high level of record completeness.
The intensity field per event is represented by the maximum wind gust for the corresponding time window of the event at each considered wind station. The local RP per event is computed by a hazard model per wind station. This is an implicit part 485 of the estimated extreme value distribution per station as explained in Section 5.1. The resulting CRPs per event and corresponding statistical tests are presented in the following Section 4.2. We have considered two weightings per station, capital, and area. Both are also computed per wind station by assigning the grid cells with capital data of the Global Assessment Report (GAR data; UNISDR, 2015) via the smallest distance to a wind station. We also use this capital data to spatially distribute our assumed total insured sum 15.23 Trillion € for property exposure (residential building, content, commercial,

The CRP of past events and validation
As announced, we have computed the CRP according to (20) with the wind gust peaks listed in the Supplementary. We considered simplified area weights and the grid's capital weights from Global Assessment Report (GAR; UNISDR, 2015), assigned to the wind station by shortest distance. The first variant considers the storm as a pure phenomenon in the geographical 505 space; the second is more focused on consequences. The Data, Table 2, and local hazard models according to (30). Our local hazard models are discussed in Section 5.1 and parameters are presented in the Supplementary Data, Table 4. We have considered two weightings for the CRP, a simple area weighting and a capital weighting (Supplementary Data, Table 3). In incompleteness of records.our record list. In the medium range, the differences between the model and empiricism are not statistically significant (for. In detail, we observe 27 storms with ≥ 1 within 20 years; expected were 20. According to the Poisson distribution the probability of 7.8% for the 27 exceedances or more for 20 years). The bias correction of local RP affects higher CRPs (is 7.8%. A two-sided test with = 5% would reject the model if this exceedance probability would be 2.5% or smaller. 520 The seasonal/annual maxima of CRP must follow a uniform Fréchet distribution ( = 1 in (3)) according to (16). We plot this and the empirical distribution in Figure 5 c). In Figure 1 d, the maximum CRP of a season follows the CDF (6) very well.. The Kolmogorov-Smirnov (KS) test (Stephens, 1986, Section 4.4) for the fully specified distribution model accepts theour model forat the extremevery high significance level of 25% (the pure KS statistic is 0.1422, sample size = 20).for the capital weighted variant. Usually, only level 5% is used; however. The KS testconsidered. This result should not be affected seriously 525 by the absence of one (probably the smallest) maximum due to an incomplete event listincompleteness issue. In summary, we state that the CRP offers a stable, testable, and robust method, to stochastically quantify winter storms over Germany.

a) comparison of area and capital weighted CRPs, b) Storm Kyrill, 2007, inverse distance weighting (by software Qgis with distance parameter 6) for local wind peaks as reproduction of the wind field and field of local RP (same colour -same percentile) , c) comparison of theoretical and observed EFexceedance frequency of capital weighted CRPs with influence of correction, dCRP, c) test of distribution of seasonal maxima of CRP with capital weighting of a seasons (equivalent with annual maxima)..
As aforementioned, the CRP corresponds with the dependence between local RPs. Therefore, the characteristic of spatial dependence between block maxima is also analyzed. The plot of the estimates of dependence measure Kendall's  (Upton and Cook, 2008) is depicted in Figure 2 b. It decreases by increasing distance. The scatter range of the half seasons is smaller than for two seasons due to different sample sizes. In Figure 2 c, we also compare Germany with results for Switzerland . For the latter, the spatial dependence is smaller, probably because of the alpine topography. The sample mean of 540 Kendall's  decreases by value 0.05 with increasing block size (half-season versus two seasons). This is statistically significant; the normally distributed confidence range of estimated expectation of the differences indicates probability 0.002 for values 0.
Max-stability is unlikely. The characteristics of spatial dependence also influence the area functions in Figure 2 c. The hypothetical case of fully spatial dependence is shown as a benchmark (CRP 100 years) and has no scattering with a CV of 0.
The scaling opportunity is demonstrated by the storm Kyrill of 2007 with = 100 years. The differences between area 545 functions of the area and capital weightings depend on the concrete storm. The CV should be independent of RP in case of max-stable dependence. This does not apply according to Figure 2  non-max-stable scaling are presented in Supplementary. In Figure 2e, the relation between scaled CRP and CV of the proxy is 550 shown. The regression function of Figure 2 d is reproduced.

Figure 2: Spatial characteristics: a) 4.3 Spatial characteristic and dependence
As discussed in Section 2.4, the spatial characteristic is an important aspect from stochastic perspective. Therefore, we have 555 analyzed the relation between distance and dependence measure. Here we have applied Kendall's  (Kendall, 1938;Mari and Kotz, 2001, Section 6.2.6) and show the dependence between half-season maxima and two season maxima for 9,870 pairs of stations in Figure 6 a. Since the sample sizes is relatively small, the spreading is strong that is caused by estimation error.
Furthermore, the differences between the estimates for one hazard season maxima and two hazard seasons maxima are almost perfectly normal distributed and should be centered to 0 in case of max-stability (CDF in Figure 6 b). This does not apply with 560 expectation 0.051 and standard deviation 0.182. According to a normally distributed confidence range for the estimated expectation with standard deviation 0.002, the probability, that the actual value is <0, is smaller than 0.00. This is in line with the differences of the non-max-stable example in Figure 2 b.
For completeness, we compare the current estimates of Kendall's  with these for Switzerland from Raschke et al. (2011) in Figure 6 c. The spatial dependence is higher for Germany. A reason might be differences in the topology. Original regression Proxy by modified scaling Max-stable scaling scaling that consider the global scaling factor and the ratio between local RP and CRP. In this way, we could reproduce the 570 observed pattern (Figure 6 f). Details of this workaround are presented in the Supplementary, Section 7. The differences between the scaling variants for storm Kyrill do not seem to be strong (Figure 6 d).

The risk estimates by averaging 580
The interest of NatCat risk analysis is in the risk curve, i.e., the relation between an event loss and its RP . It is a bijective Before we estimated risk curves according to the approach of Section 3.4, we must estimate a vulnerability function ( ) with its inverse variant ( ) and can be estimated by relation (5)  (11) 590 Further details of the concept are explained and illustrated in Section 4.
Our demonstration object remains Germany, and we estimate its risk curve for insured winter storm losses. The local losses are simply the product of local exposure value and loss ratio . This is determined by local maximum wind speed per event via the (local) vulnerability function (loss or damage function). Our loss function follows the approach of Klawa and Ulbrich (2003) and is fitted by the data of the General Association of German Insurer (GDV, 2019) for 16 historical events from 2002 595 to 2018 (Figure 3 a). These loss data are already inflated/scaled to the German insurance market 2018. An exemplary variant of our vulnerability functions is compared to previous models in Figure 3 b; different geographical resolutions might explain the differences. The processed model of the geographical distribution of exposed values is the same as for capital weighting in the previous section. The assumed total sum insured (as plotted in Figure 7  The differences between max-stable and non-max-stable scaling in the risk estimates are demonstrated in Figure 7 c. For smaller RP, no significant difference can be stated in contrast to higher RP. This corresponds with the differences between the CV in relation to the CRP for max-stable and non-max-stable case in Figure 6 f. These are also higher for higher CRP. 610 We also compare our results with previous estimates in Figure 7 c. For this purpose, we must scale these to provide comparability as good as possible. The relative risk curve of Donat et al. (2011) is scaled simply by our TSI) is 15.23 Trillion euros and based on the assumptions of assumption for exposure year 2018. The vendor models of Waisman (2015). His) are scaled by the average of ratios between modelled and observed event losses from storm Kyrill since a scaling via TSI for Germany is scaled/inflated from 2015 to 2018 by factor 1.12 under consideration of increasing construction prices 615 (Statistisches Bundesamt, 2020) and building stock (GDV, 2019).
The three risk curves, that have been estimated by (9) with wind fields scaled by capital weighted CRP are shown in Figure   3 c. Also, reported losses of historical events 0 and their empirical RP (observation period of 17 years) and the estimated CRP are plotted. The differences between the risk curves are minor. The standard error (SE) of estimates by loss averaging (also called RP scaling) is more minor than 8.5% of the estimated event loss for RP ≤ 400 years. The range of two SE is not a 620 confidence range for empirical RP. Its uncertainty is higher and explains the higher differences to the risk curves.
Since German winter storms do not imply max-stable dependence, we also estimate the risk by non-max-stable scaling.
The result is depicted in Figure 3 d and shows lower event losses for high RP than in the max-stable case. We also compare our risk curves with these vendor models, which have been presented by Waisman (2015). However, the exposure assumptions are not the same (exposure year,was not possible (uncertain market share, and split between residential, commercial, and 625 industry). To provide comparability, the vendor estimates are scaled by factor 1.34. This is the average of the ratios between reported and modeled loss for storm Kyrill. We also consider the risk estimate by Donat et al. (2011) and exposure). The result of the standard model of European Union (EU) regulations (European Commission, 2014), also known as Solvency II requirements. It considers the exposure per CRESTA zone (, is also based on our TSI assumption, split into the Cresta zones by the GAR data. The Cresta zones (www.cresta.org). For the standard model, we split our assumed TSI to the Cresta zones 630 by the GAR data. The same TSI is used to scale Donat's et al. (2011) estimated RP for event loss ratios (NCEP variant in Figure 7); our extrapolation complies with their distribution assumption. We highlight that we cannot fully ensure comparability to the previous risk estimates) are an international standard in insurance industry and corresponds to the two digit postcode zones in Germany.
The risk estimate of Donat et al. (2011) is based on a combination of frequency estimation and event loss distribution by 635 the generalized Pareto distribution which is fitted on a sample of modelled event losses for historical storms. The corresponding risk curve differs very much from other estimates and obviously overestimate the risk of winter storms over Germany. The standard model of EU only estimates the maximum event loss for RP 200 years, the estimated event loss is very high. The vendor models vary but have a similar course as our risk curves. The non-max-stable scaling is in the lower range of the vendor models, the unrealistic max-stable scaling is more in the middle. The concrete names of the vendors can be found in Waisman's 640 (2015) publication. The reader should be aware that the vendors might have updated their winter storm model for Germany in the meantime.
The major result of Section 5 is the successfully demonstration that the CRP can be applied to estimate reasonable risk curves under controlled stochastic conditions. We have also discovered the high influence of the underlying dependence model (max-stable or not) and corresponding spatial characteristic to loss estimates for higher ELRP. 645

24
5 Technical details of the application example 650

Modelling and estimation of local hazard
As aforementioned, the maximum wind gusts of half seasons of winter storms (extratropical cyclones), the block maxima, have been analyzed. Therein the generalized extreme value distribution (Beirlant, et al., 2004, (5.1)) is applied As discussed below, the Gumbel distribution (Gumbel, 1935(Gumbel, , 1941, as a special case in (29) with extreme value = 0, is 655 an appropriate model. The scale parameter is , location parameter is . The local hazard function (13,14) can be derived directly from estimated variant of (29) according to the link between extreme value distribution and exceedance frequency (16); the accent symbolizes the point estimation) We apply the Maximum likelihood method for the parameter estimation (Clarke, 1973, Section 2.6.3). The 660 incompleteness of wind records per half season have been considered in the ML estimates by a modification of the procedure as explained in the Supplementary, Section 5. A Monte Carlo simulation confirms the good performance of our modification.
The biased estimate of for our sample size = 40 was also detected which we considered ̂=̂/0.98 as corrected estimation. Landwehr et al. (1979) have already stated such bias. A further bias was discovered, the EF is well estimated by (30) in contrast to the RP ̂, this is strongly biased. We also corrected this as documented in the Supplementary, Section 6. 665 The analyzed half-season maxima, record completeness and parameter estimates are listed in Supplementary Data, Table 4, 5 ad 6.
We have validated the sampling of block maxima per half season. The opportunity of correlation between the first and second half-season maxima has been tested, for significant level  = 5% around 6% fails the test with Fisher's ztransformation (Upton and Cook, 2008). This corresponds to the error of the first kind and is interpreted as correlation not 670 being significant. Similarly, the Kolmogorov-Smirnov homogeneity test rejects 4% of the sample pairs first half to second season half season for a significant level of 5%. This 5% are the expected share of falsely rejected correct modelsthe first kind of error (type I error, e.g. Lindsey, 1996, Section 7.2

.3).
To optimize the intensity measure of the hazard model, we have considered the wind speed with power 1, 1.5, and 2 as the local event intensity in a first fit of the Gumbel distribution by the maximum likelihood method. According to these, power 675 1.5 offers the best fit of wind gust data to the Gumbel distribution. Such wind measure variants were already suggested by Cook (1986) and Harris (1996).
We do not apply the generalized extreme value distribution in (24) with extreme value index ≠ 0 but the Gumbel case with = 0 for the following reasons. At first, an extensive physical explanation would be required if some wind stations are concerned by a finite upper bound for < 0 and other stations not with ≥ 0 according to (29). Why should be local wind 680 hazard short tailed for some wind stations and heavy tailed for others? River discharges at different gauging stations could imply such physical differences since there variants with laminar and turbulent stream (catchword Reynolds number) or very different retention/storage capacities of catchment areas (e.g., Salazar et al. 2012). Such significant physical differences do not exist for wind stations which are placed and operated under consideration of rules of meteorology (World Meteorological Organization, 2008, Section 5.8.3) to provide homogeneous roughness condition due to generate comparable data. Besides, 685 we also found several statistical indications for our modelling. Information criteria AIC and BIC (Lindsey, 1996; here over all stations) indicate that the Gumbel distribution is the better model than the variant with a higher degree of freedom. Furthermore, the share of rejected Gumbel distributions of the Goodness-of-fit test (Stephens, 1986, Section 4.10) is with 6% around the defined significance level of 5% (the error of 1 st kindfalsely rejected correct models). We have also estimated for each station and got a sample of point estimates. The sample mean of is with 0.002 very close to = 0 which confirms our 690 assumption. Moreover, the sample variance is 0.018 which is around the same what we get for a large sample of estimates ̂ for samples of Monte Carlo simulated and Gumbel distributed random variables ( = 40). All statistics validate the Gumbel distribution.

Modelling and estimation of vulnerability
To quantify the loss ratio at location (wind station) and event in the loss aggregation (26), we use the approach of Klawa 695 and Ulbrich (2003) for Germany. The difference to the origin is not significant. The event intensity is the maximum wind gust speed. 98% is the upper 2% percentile from empirical distribution of all local wind records. The relation with vulnerability parameter is Donat et al. (2011) have used a similar formulation but with an additional location parameter. This is discarded here since the 700 loss ratio must be = 0 for local wind speed < 98% . This is also a reason, why a simple regression analysis (Fahrmeir et al., 2013) is not applied to estimate . We formulate and use estimator with historical storms, corresponding reported event losses , wind stations, and local exposure value , being assigned to the wind station. , be fixed for every station if there were wind records for every storm at each station . However, the 705 wind records are incomplete and the assumed TSI must be split and assigned to the stations a bit differently for some storms.
The exposure share of the remaining stations is simply adjusted so that the sum over all stations remains the TSI.
Our suggested estimator (32) has the advantage that it is less affected by the issue of incomplete data (smaller events with smaller losses are not listed in the data) than the ratio of sums over all events, and the corresponding standard error can be quantified (as for the estimation of an expectation). The current point estimate is ̂=9.59E-8±5.97E-9. 710 26 An example of our vulnerability function (with average of 98% over the wind stations) is depicted in Figure 8 and compared with previous estimates for Germany. It is in the range of previous models. Differences might be caused by different geographical resolutions of corresponding loss and exposure data. A power parameter of 2 in (31) might also be reasonable since the wind load of building design codes (European Union, 2005, Eurocode 1) is proportional to the squared wind speed.
The influence of deductibles (Munich Re, 2002) per insured object is not explicitly considered but smoothed in our approach. 715   Heneka et al. (2006) Munich Re (2002), homeowner, 1999, upper bound Munich Re (2002), homeowner, 1999 3 Secondary Methods in the Analysis of German Winter Storms

Modeling of local hazard
The local hazard curve for the EF of local event intensity must be estimated to compute the local RP by event and station 725 for our examplewinter storm over Germany. For this purpose, the relation between local event frequency and CDF of corresponding block maxima is used according to equation (3) and (6). We model the CDF of local block maxima of intensity by the Gumbel distribution. This is a special case of generalized extreme value distribution Beirlant, et al. 2004) with extreme value index = 0 in the CDF . 730 (12)

Numerical procedure of scaling
We briefly explain here the numerical procedure to calculate a risk curve via averaging the event loss. For any supporting point of a risk curve during an event loss averaging, the ELRP is defined and determine the scaled CRP for all historical 735 events. For each historical event, the scaling factor is = / according to (24) and is applied in (27)  The historical events are also scaled for a defined event loss and the corresponding scaled CRP is averaged. However, the Goal Seek function in MS Excel is applied to find the correct scaled CRP and corresponding scaling factor . For the averaging of the exceedance frequency, the reciprocal of is averaged. All these apply to max-stable scaling. For the non-745 max-stable scaling the scaling factor is adjusted to a local variant according to the description in the Supplementary, Section 7. Therein, the factor is adjusted for each station and depends on the relation of local RP to CRP of the historical event. This adjustment is made for each historical storm individually.

Error propagation and uncertainties
The uncertainty of the local hazard models influences the accuracy of the CRP since the CRP is an average of estimates of 750 local RP. The issue is that there is a certain correlation between the estimated hazard parameters of neighboring wind stations.
We consider this by application of the Jack-knife method (Effron and Stein, 1986). According to these, the mean squared error (MSE, which is the standard error if the estimate is bias free as we assume here) of the original estimated parameter ̂ is (accents symbolize estimations) with the estimates ̂− for the Jack-knife sample of observations, being the original sample but without one of the observations/realizations. Therefore, it is also called also called leave-one-out method. The estimator (33) implies a parameter sample of ̂− of size , with one estimated parameter or parameter vector for each Jack-knife sample of observations.
To consider any correlation in the error propagation of CRP estimate, the maximum of the same half-season is left out synchronously when the parameter sample is computed for each wind station. Without changing the order in the parameter 760 sample of each wind station, the CRP ̂− of the concrete historical event is computed with the hazard parameters ̂− of each station. Finally, for this storm, the standard error of point estimate ̂ is computed according to (33).
We use the same approach to consider the error propagation from local hazard models to the risk estimate for the maxstable case in Section 4.4. But the finally estimated parameter ̂ in (33) is the averaged event loss ̂( ) for scaled CRP. This only covers a part of uncertainties in risk estimate. We consider two further sources of uncertainty and assume that they 765 influence the risk estimate independently to each other. The uncertainty of loss averaging is the same as during an estimation of an expectation from a sample mean and is determined by sample variance and sample size (number of scaled events). The propagation of the uncertainty of the vulnerability parameter is computed via the delta method (Coles 2001, Section 2.6.4).
The aggregated standard error is the square root of the sum of squared errors. This implies a simple variance aggregation according to the convolution of independent random variables (Upton and Cook, 2008). 770 The computed standard errors in Figure 7 b are in the range 7.5 to 8.5% of estimated event loss per defined ELRP. The scale parameter is , location parameter is . For completeness, this distribution is max-stable in a univariate sense. Extreme value index is independent of the block size ( in equation (6)); is independent of the Gumbel case's block size. For a sampling of block maxima, we only consider wind speed maxima of the winter storm season that we define from September to April. A shorter definition is from October to March. The differences between the resulting block maxima are not significant. 775 To increase the sample size, we analyze the maxima of half seasons of Autumn 1999 to Spring 2019 and have 40 observations. An entire season corresponds to one year. The distinction between block size (number of unit periods in (6)) were considered in the post-processing of the results. There is the opportunity that other kinds of wind perils contaminate our sample. As a specific compensation, extratropical cyclones also (albeit rarely) occur outside our sampling periods.
We consider records of wind stations in Germany of DWD (2020; FX_MN003, a daily maximum of wind peaks, usually 780 wind guest speed [m/s]) that include minimum record completeness for one of the analyzed storms of 90%, at least 90% completeness for all seasons and 55% minimum completeness per (half) season. Therefore, we only use records of 141 stations.
The autocorrelation between the subsequent half-season maxima has been tested, for significant level  = 5% around 6% fails the test. This corresponds to the error of the first kind what is interpreted that autocorrelation is not significant. Similarly, the homogeneity test (two-sample test; Conover, 1971) rejects 5% of the pairs of samples first season half versus second season 785 half for a significant level of 5%. Besides, we have considered the wind speed with power 1, 1.5, and 2 as the intensity in a first fit of the Gumbel distribution by the maximum likelihood method (Coles, 2011). According to these, power 1.5 offers the best fit of wind gust data. Such wind measure adjustment was already suggested by Cook (1986) and Harris (1996).
We do not consider the generalized extreme value distribution with index ≠ 0 in (12) for the following reasons. It would require an extensive physical explanation if some wind stations are concerned by a finite upper bound for < 0 and other 790 stations not with ≥ 0. Besides, information criteria AIC and BIC (Lindsey, 1996) indicate that the Gumbel distribution is the better model. In addition, the share of rejected Gumbel models (one by station) of the Goodness-of-fit test (Stephens, 1986) is with 6% around the defined significance level of 5%, the error of 1st kindfalsely rejected correct models. Also, when we estimate the extreme value index , then the standard deviation of all estimates equals the SE for estimates for actual index = 0 for the same sample size. All statistics indicate the Gumbel distribution. 795 Even though we considered only stations with a high level of record completeness, not every sample is complete; some daily records are missed. This affects the sampling of block maxima and indirect likelihood (ML) estimation for their CDF. This is considered in the estimates as explained in the Supplementary. An extensive Monte Carlo simulation with 100,000 repetitions for sample size = 40 has been done to validate the correction's performance. We also realized that complete and incomplete samples result in a biased estimate of ; it is only 98% of the actual value. This ML method's bias for small sample 800 sizes is already stated by Landwehr et al. (1979) and is corrected for our application. Furthermore, EF ( ) estimates for a defined intensity level are (more or less) unbiased. This does not apply to the estimated RP as reciprocal of EF (parameters corresponds to (12)) Via the estimation results of the Monte Carlo simulation, a correction function for the local RP estimates has been derived and is presented in the Supplementary. The correction function is also applied to estimates for sample size = 39 because of the small difference. The smaller sample is applied to estimate the SE of CRPs by the leave-one-out method, also called 810 Jackknife method (Lindsey, 1996).

Modeling of vulnerability
To quantify the loss ratio at location (wind station) and event , we use the approach of Klawa and Ulbrich (2003) for Germany. The difference to the origin is not significant. However, our variant has a bit higher correlation between reported and modelled loss. The event intensity is the maximum wind gust speed. 98% is the upper 2% percentile from empirical 815 distribution of all local records. The relation is with parameter 820 Donat et al. (2011) have also used this approach but with an additional location parameter. This is discarded here since the loss ratio must be = 0 for local wind speed = 0. This is also a reason, why a simple regression analysis (Lindsey, 1996) is not applied to estimate . We use with historical storms, corresponding reported event losses , wind stations, and (here modelled) local exposure value , being assigned to the wind station. It would be the same for every station if there were wind records for every storm at each station. However, the wind records are incomplete and the assumed TSI must be split and assigned to the stations a bit differently for some storms. The normal share per station is pre-estimated by capital grid values of GAR data (normalized that 830 sum is 1). These have been assigned to the station by the shortest distance. Our estimate for the scale parameter by (15) is 0.0960 ± 0.0060. The power parameter (exponent) and upper 2% percentile in (14) are defined and influence the risk estimate.
A higher power parameter would result in higher event losses for higher RP. However, an exponent larger than 3 is not likely since an exponent of 2 would be also reasonable according to building codes (European Union, 2005) with proportional relation between squared wind peaks and structural load. 835

Error estimate and accuracy
The SE of the winter storm risk estimate for Germany is only computed for the RP scaling (loss averaging). Thereby we consider three components. The uncertainty from the local hazard estimates is considered by the Jackknife method, applied synchronously by half-season for all wind stations to consider the correlation between the SE of different stations. The SE for the vulnerability function's scaling parameter in (14) and (15) is computed as the standard error of expectation being estimated 840 31 by the sample mean. The same applies to the loss estimated by scaled CRP of 16 historical storms. The three SE are combined under the assumption of independence as the sum of their corresponding variances (square of SE).
The shares of uncertainty components on the error variance (squared SE) of our risk estimates depend on the RP. On average offor our supporting point, these are 15% for the limited sample of scaled historical events, 24% for the uncertainty of local hazard parameters, and 61% by the vulnerability model's parameter. We do not know a published error estimation for 845 a vendor model for risk from winter storm over Germany and can only compare our estimates with these of Donat's et al. (2011). Their confidence range indicate a smaller precision than ours.
Our two weighting variants for Winter Storms over Germany correlate much better with Kendall's = 0.814 and Spearman's = 0.946 than the RPs of Della-Marta's 0 et al. event measures for European windstorms. The accuracy of the current RP estimates is also slightly higher than these of Della-Marta et al. For instance, an estimated RP 15-20 years has a SE 850 smaller than two years currently; Della-Marta et al. report RP a range of 14 years for the 95% confidence interval. This is approximately equivalent to SE of more than three years. The current precision of the German storm risk estimate with a sample of 16 historical storms is also higher than the estimates by Donat et al. (2011) with a sample of 30 historical storms.

RP of vendor's risk estimate
We comparehave compared our results with vendor models in Section 4.4. These have estimated the risk curve for the 855 maximum event loss within a year. This is a random variable, and their pseudo-RP is the reciprocal of the exceedance probability. This pseudo RP and can never be smaller than 1. Under the assumption of a Poisson process, we transform the pseudo RP to an actual and event related one bywith the relation between EF and CFCDF in (6).13,16) With increasing event loss, the difference between its pseudo RP and the actual oneELRP converges to 0. The concrete name of the vendors can be found in Waisman's (2015) publication. The reader should be aware that the vendors might have updated their winter storm 860 modelrelative differences are around 5% for Germany in the meantimeELRP 10 years and 0.5% for 100 years.

Stochastic background and interpretation 870
The estimation is based on following stochastic relations and assumptions (or proxies) The origin is (5); the well-known delta method (Coles, 2011) for computation of propagation of errors is also a base. A more illustrative explanation is provided for the loss scaling by Figure 6 a. The sample mean (sample average) of corresponding 875 CRPs is the estimate for their expectation and thereby the estimate of RP of event loss. Figure 6 b offers a further stochastic interpretation. The basic set includes all possible events which match with the risk curve. In addition, every possible wind field has a CRP. This also provides a link between some of the possible event fieldsthe CRP scaling and pools wind fields to subsets. The actual probability density of an observation of a random variable is mostly not known. In the same way, each "observed" subset has an unknown probability (density). And each subset corresponds to a risk curve. Their expectation 880 determines the risk curve of the entire set of event fields and is estimated by the average of a sample of linked event fields since the expectation of a random variable is estimated by the corresponding sample mean. Three averaging variants are possible.

General
CRP opportunities prove once again that mathematical statistics and stochastic is the central technology for analysis, modeling, and understanding risks such as NatCat. The CRP is a simple, In the beginning, we asked the questions about the RP of a 890 hazard event in a region, the corresponding NatCat risk, and necessary conditions for a reasonable NatCat modelling. To answer our questions, we have mathematically derived the CRP of a NatCat generating hazard event from previous concepts of extreme value theory, the pseudo polar coordinates (17). This implies the important fact that the average of the RPs of random point events remains a RP with exceedance frequency (8, 15). Furthermore, we extended Schlatter's 1st theorem for max-stable random fields to non-max-stable spatial dependence and characteristic. We have also considered the normalized 895 variant of the area function of all local RP of the hazard event in a region with parameters CRP and CV. The absence of maxstability in the spatial dependence results in correlation between CRP and CV, which is a further indicator for non-max-stability beside changes of measures for spatial dependence by changed block size (e.g., annual maxima versus two years maxima).
The derived CRP is a simple, plausible, and testable stochastic measure for a catastrophe. In addition, hazard and NatCat event. The weighting of local RP in the computation of the CRP can be used to compensate an inhomogeneous distribution of 900 corresponding measuring stations if the physical-geographical hazard component of a NatCat, the event intensity field, is of interest. However, the concentration of human values in the geographical space could also be considered in the weighting to get a higher association of the CRP with the ELRP of a risk curve. This link implies the conditional expectation (18) under the assumption of max-stable association between CRP and ELRP and offers the new approach offers the opportunity to estimate risk curves, the bijective function event loses to ELRP, via a stochastic scaling of historical intensity fields and averaging of 905 corresponding risk parameters. The averaged parameters can be the scaled CRP for a defined event fields under stochastic control. By this and a simple exposure and vulnerability model, the risk curve forloss, corresponding exceedance frequency, or the event loss for a defined/scaled CRP.
The differences between the three estimators are small in our application example, insured losses from winter storm losses inover Germany could be derived.. In contrast to this, the influence of the stochastic assumptions regarding spatial dependence 910 and characteristic (max-stable or not) is significant in the range of higher ELRP. This highlights the importance of realistic consideration of spatial dependence and characteristic of the hazard in a NatCat model. Besides, our risk curves for Germany have a similar course as thesethose derived by vendors (Figure 7 d). The risk assumption by EU for Germany with RP 200 years is significantly higher than ours. The estimate by Donat et al. (2011) differs significantly and seems to be implausible.
for higher RP. A reason might be their statistical modelingmodelling by the generalized Pareto distribution as already applied 915 for wind losses by Pfeifer (20112001). The tapered Pareto distribution (Schoenberg and Patel, 2012)), also called tempered Pareto distribution (Albrecher et al., 2021), or a similar approach (Raschke, 2020) is aprovide more appropriate proxyproxies for our risk curve's tail.
According to our results, necessary conditions for an appropriate NatCat modelling are the realistic consideration of local hazard and their spatial dependence (max-stable or not?). Correspondingly, the spatial characteristic of NatCat events, 920 described here by relation CRP to CV, must be reproduced. In addition, the CRPs of a simulated set of hazard events in a NatCat model should have an empirical exceedance frequency that follows the theory (15). That the standard error of an estimate should be quantified, the sampling should be appropriate, and overfitting should be avoided (catchwords over parametrization and parsimony), applies to all scientific models with a statistical component (e.g., Lindsey, 1996). A reasonable weighting for the CRP might depend on the examined peril. The expected annual sum of local losses might be a universal weighting since the sum of local expectations determines the expectation of total annual loss. The various weighting opportunities also underline the fact that there is not an absolute view on hazard and risk.

Requirements of the new approachapproaches 935
Our approach to CRP is based on two assumptions. At first, the local and global events occur as a Poisson process. This is frequently useda common assumption or approximation in extreme value statistics and its application ) and, Chapter 7) and the corresponding Poisson distribution of number of events can be statistically tested (Stephens, 1986;Section