Extreme events and emergency scales.

An event is extreme if its magnitude exceeds the threshold. A choice of a threshold is subject to uncertainty caused by a method, the size of available data, a hypothesis on statistics, etc. We assess the degree of uncertainty by the Shannon's entropy calculated on the probability that the threshold changes at any given time. If the amount of data is not sufficient, an observer is in the state of Lewis Carroll's Red Queen who said "When you say hill, I could show you hills, in comparison with which you'd call that a valley". If we have enough data, the uncertainty curve peaks at two values clearly separating the magnitudes of events into three emergency scales: subcritical, critical, and extreme. Our approach to defining the emergency scale is validated by 39 years of Standard and Poor's 500 (S&P500) historical data. Published by Elsevier B.V.

Introduction

Not a single day passes by without hearing about extreme events which surround us almost everywhere. On one hand, climate change results in droughts, heat waves, tornadoes, storms etc; movement of the tectonic plates is responsible for earthquakes and volcano eruptions. On the other hand, certain aspects of human activities may crash stock markets, influence tensions not only among groups of people in a country, but also between countries sometimes leading to military confrontations and mass migrations, etc. While there is no single definition of extreme events [1], they are considered events that cause infrastructure failures, economic and property losses, risk to health and life.

In order to quantify extreme events, practitioners developed several scales. For instance, Modified Mercalli Intensity scale [2] (describes the intensity of visible damage of earthquakes), Beaufort Wind scale [3] (measures speed and observed effects the wind has on the sea and land), Saffir-Simpsons Hurricane scale [4] (measures wind speed), Fujita scale [5] (rates the intensity of a tornado after it has passed), US Homeland Security Terror Alert scale [6] (measures five color-coded terror alert levels), U.S. Climate Extremes Index [7], etc. Rohn Emergency Scale [8] unites emergency scales using three independent dimensions: (i) scope; (ii) topographical change (or lack thereof); and (iii) speed of change. The intersection of the three dimensions provides a detailed scale for defining any emergency [8]. In some papers, the threshold for an extreme event is related to the number of standard deviations from the average amplitude [9], [10]. However, existing empirical scales tend to describe the characteristics of the event itself rather that the consequences; such scales are ill-suited to describe emergencies in a way that is meaningful for response [11].

For instance, extreme events of different magnitudes in financial markets range from global recessions (defined by a global annual GPD growth rate of 3.0 percent of less) that happened in 1975, 1982, 1991, 2009 to “flash crashes” (for instance, on May 6, 2010, the S&P500 declined 7% in less than 15 min, and then quickly rebounded). Unlike the Richter magnitude scale, the severity of flash crashes in financial markets is defined by the measures that need to be taken to ease panic such as halting trading. Under 2012 rules, market-wide circuit breakers (or ‘curbs’) kick in when the S&P 500 index drops 7% for Level 1; 13% for Level 2; and 20% for Level 3 from the prior days close. A market decline that triggers a Level 1 or 2 circuit breaker before 3:25 p.m. Eastern Time will halt trading for 15 min, but will not halt trading at or after 3:25 p.m. Circuit breakers can also be imposed on single stocks as opposed to the whole market. Under current rules, a trading halt on an individual security is placed into effect if there is a 10% change in value of a security that is a member of the S&P 500 Index within a 5-minute time frame, 30% change in value of a security whose price is equal or greater than $1 per share, and 50% change in value of a security whose price is less than $1 per share [14]. Ironically, in August 2015, single stock circuit breakers produced unprecedented disruption as 327 exchange-traded funds experienced more than 1000 trading halts during a single day.

For the short-term reactions of stock markets, measured in terms of returns, there exist several approaches to defining their severity. For instance, some authors [12] define extreme events of stocks if prices change greater than 2.5%, for both positive and negative returns. However, others [13] suggest a threshold of 10%.

In general, statistics of extreme events show that the extreme events are found in the tails of probability distributions (i.e. the distributions extremities). Inference over tails is usually performed by fitting an appropriate limiting distribution over observations that exceed a fixed threshold. The choice of the threshold over which to fit a probability distribution is hard and arbitrary. Although tools to guide this choice exist, inference can greatly vary for different thresholds. In addition to that, different distributions may admit asymptotic tails of the same appearance (a power law). In our paper we apply two approaches existing in practical extreme value analysis to study S&P 500 time series in the period from January 2, 1980 till December 31, 2018.

The first one relies on deriving block maxima series as a preliminary step. We fit the actual data to the Generalized Extreme Value Distribution (GEVD). Our results (Section 4.1) show that the distribution of block maxima is a composition of several distributions (Fig. 6b).

Fig. 6

(a) A general extreme value QQ-plot with maximum likelihood estimation; (b) Density plot of empirical data where a dashed curve A is based on the empirical data, and a dashed curve B is modeled. and bandwidth is 135.9.

The second approach relies on extracting the peak values reached for any period during which values exceed a certain threshold (falls below a certain threshold). We will cover major methods of choosing the threshold in Section 4.2. Usually with the second approach the analysis may involve fitting two distributions: one for the number of events in a time period considered and a second for the size of the exceedances. But in practice a suitable threshold is unknown and must be determined for each analysis [15].

We demonstrate several methods of choosing a threshold. An example of an empirical method, the rule of thumb, is studied in Section 4.2.3. A choice of the threshold can be made using statistical analysis tools, such as graphical approaches (Section 4.2) and automatic methods (Section 4.2.2). Rather than characterize extreme events solely in terms of the distribution of weekly maxima, threshold methods take account of all events beyond a given high threshold (Fig. 1 ).

Fig. 1

Degree of uncertainty of different values of the thresholds based on the amount of trading days taken into account. Three shaded regions mark three scales of emergency: I - subcritical, II - critical, III - extreme and the solid curve represents uncertainty of the Red Queen State.

Every decision we make or approach we choose is made in the face of uncertainty. For instance, as of March 25, 2020, 414, 179 cases and 18, 440 deaths due to coronavirus disease 2019 (COVID-19), caused by the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), had been reported worldwide [16]. In attempts to slow down spread of the disease, the economies of many counties slammed brakes and financial markets plummeted in the light of uncertainty, but the above does not describe severity of this particular extreme event, but rather, a reaction of the mankind to COVID-19. One natural characteristic of severity of pandemic is death rate, which in case of SARS-CoV-2 ranges from 0.66% [17] to 4% [18]. Mortality rate is higher among people with underlying medical conditions, sometimes, deaths partly occur due to coronavirus. Therefore, some believe that this pandemic is poised to be more dramatic than anything we have seen in our lifetimes [19], but other compare the situation to an elephant being attacked by a house cat. Frustrated and trying to avoid the cat, the elephant accidentally jumps off a cliff and dies [20].

An uncertainty depends on amount of data considered, a method used, etc. Moreover, once new data becomes available, we cannot guarantee sustainability of the previously chosen threshold value. In addition to that, assessment of the severity of the extreme event is also ruled by consequences, which may be revealed ex post facto. In Section 4.3 we presented statistics of extreme events under threshold uncertainty. We analyzed the degree of uncertainty of the threshold value (chosen with the rule of thumb) based on the probability of its change at any given day. Fig. 16 demonstrates the degree of uncertainty that depends on both the amount of data available and the values of thresholds. We observed several cases:

Fig. 16

Degree of uncertainty of different values of the thresholds based on the amount of trading days taken into account. Three shaded regions mark three scales of emergency: I - subcritical, II - critical, III - extreme and the solid curve represents the Red Queen State.

If the amount of data is not sufficient (a solid line corresponding to an 18-day window of obervation in Fig. 16), then the uncertainty curve forms a skewed profile attaining a single maximum for some value of the threshold. In this situation, an observer’s perception of events reminds Red Queen from “Through The Looking-Glass and What Alice Found There” by Lewis Carroll [21] who said “When you say hill, I could show you hills, in comparison with which you’d call that a valley”. Our understanding of the events whether they are extreme or not is very limited and uncertainty is blurry. As events become more severe, our uncertainty that the events are extreme decreases. The observer realizes that the events are extreme, but a precise point at which the events turn to be severe cannot be determined. This case is called the Red Queen State.

As the window of observation becomes larger, the uncertainty curve exhibits two maxima indicating that the amount of data is sufficient. As we further extend the window the curve torrents into sharp peaks (Fig. 16). The latter ones clearly separate the threshold values into three regions: three levels of emergency.

In the region I (subcritical), the threshold values are small to raise concern about extreme events. Then we can observe a spike with the degree of uncertainty attaining its first maximum. This extremum indicates a transition to the next kind of uncertainty.

In the region II (critical), uncertainty is conceptualized with a question whether a magnitude of the event is already critical, extreme or not yet. Further, we see another jump of uncertainty. At this point we certain that events are not regular anymore.

We consider all events in this region extreme with our uncertainty decreasing. If these events are not extreme then what are they?

We consider the Red Queen state and three scales, which are based on the degree of uncertainty, our contribution to the discussion on the extreme events and emergency scales. We conclude in the last Section 5.

Data source and description

The analysis of extreme events in the S&P 500 time series of log-return has been performed based on the data collected during 9835 trading days (39 years), in the period from January 2, 1980, till December 31, 2018 (see Fig. 2 ) acquired using the publicly available source at Yahoo Finance (https://finance.yahoo.com/quote/%5EGSPC/).

Fig. 2

The S&P 500 Index of 500 large-cap U.S. stocks assessing market performance.

The data set contains S&P 500 index at market open, highest point reached in the day, lowest point in the day, index of the stock at market close, number of shares traded. We have used index at market close since this information was always present in the data set. Computations were made using Pythons numerical libraries, such as NumPy and Pandas, as well as R-language to perform statistical analysis.

Log returns: Between noise and random walks

Due to complexity of financial data, in order to analyze dynamics of the given financial time series, we computed log return, denoted R ln, according to the following formula [22]:where is the value of S&P 500 at market close (Fig. 3 ).

Fig. 3

The S&P 500 index log-return at market close.

The lowest drop of log-return in Fig. 3 happened Monday, October 19, 1987 known as Black Monday when the S&P500 shed value of nearly 23% [23]. Another significant drop occurred in 2008 when S&P 500 fell 38.49%, its worst yearly percentage loss [24]. In September 2008, Lehman Brothers collapsed as the financial crisis spread. However, on Oct 13, 2008: S&P 500 marks its best daily percentage gain, rising 11.58% and registers its largest single-day point increase of 104.13 points [24].

Distribution of R ln values (Fig. 4 ) is asymmetrically skewed, with fat heterogeneous right and left tails.

Fig. 4

Distribution of log-return values in the log 2-linear scale. Solid lines correspond to a Zipf’s Law (), dotted lines represent Power Law (), and dashed lines correlate to Gaussian distribution (). Curves are given for reference only.

The log-return time series has a scale invariant structure when the structure repeats itself on sub-intervals of the signal, where the Hurst exponent H characterizes the asymptotic behavior of the auto-correlation function of the time series [25], [26], [27]. The larger Hurst exponent is visually seen as more slow evolving variations (i.e., more persistent structure) of the time series [26], [28], [29]. Processes with 0 < H < 0.5 exhibit antipersistence, with an increase in the process is likely to be followed by a decrease in the next time interval resulting in sample paths with a very rough structure [26], [28], [29]. On the contrary, values 0.5 < H < 1 lead to long-range dependence (“long memory”) in the time series, with more likely the same sign of successive increments (persistence) and smoother sample trajectories. Finally, the time series constitutes random walks when H > 1 that have more apparent slow evolving fluctuations [26], [28], [29]. The q-order Hurst exponent H is only one of several types of scaling exponents used to parameterize the multifractal structure of time series [26], [30].

The log-return time series for S&P 500 exhibits local fluctuations with both extreme small and large magnitudes, as well as short- and long-range dependences on different time scales [31], [32]; it is not normal distributed and all q-order statistical moments should to be considered to describe the spatial and temporal variation that reveals a departure of the log-return time series from simple random walk behavior [26], [28]. The q-order weights the influence of segments with large and small fluctuations. The negative q′s are influenced by the time series segments with small fluctuations, and large fluctuations influence the time series segments for positive q’s.

In our work, we use the standard multifractal detrended fluctuation analysis (MFDFA) algorithm [26], [30] for estimating the q-order Hurst exponents and the multifractal spectra directly from the time series:

The original time series x, is aggregated by computing the cumulative sums where ⟨x⟩ denotes the sample mean;

The aggregated data is divided into ⌊N/s⌋ non-overlapping segments of length s;

The maximum likelihood estimator of the residual variance in segment ν, where y(i) is the degree polynomial fitted the aggregated observations in the segment;

For each segment of length s and for each positive or negative values of the moment order q, the q-order fluctuation function,

are calculated. The local fluctuations F with large and small magnitudes is graded by the magnitude of the negative or positive q-order, respectively;

A linear regression of ln F(s) on ln s for all s is performed, and the slope of the linear function ln F(s)∝Hln s is used as an estimator of the q-order Hurst exponent H for each q-order fluctuation function F.

The fractal structures of the positive and negative log-return time series and its deviations within time periods with large and small fluctuations are assessed by the q-order Hurst exponents (see Fig. 5 ).

Fig. 5

The q-order Hurst exponents H for the time series of positive (the dashed line) and negative (the bold line) log-returns.

The slopes H of the regression lines are q-dependent for the multifractal time series of positive (the dashed line) and negative (the bold line) log-returns. (see Fig. 5). Decreasing H with the q order indicates that the segments with small fluctuations have a random walk like structure whereas segments with large fluctuations have a noise like structure.

Tails, thresholds, and extreme events

There are two primary approaches to analyzing extreme values (the extreme deviations from the median of the probability distributions) in data:

The first and more classical approach reduces the data considerably by taking maxima of long blocks of data, e.g., annual maxima. The GEVD function has theoretical justification for fitting to block maxima of data [33].

The second approach is to analyze excesses over a high threshold. For this second approach the generalized Pareto distribution (GPD) function has similar justification for fitting to excesses over a high threshold [33].

Generalized extreme value distributions

The GEVD is a flexible three-parameter continuous probability distributions that was developed with extreme value theory to combine the Gumbel, Fréchet, and Weibull extreme values distributions into one single distribution [34], [35]. The GEV distribution has the following pdf [36]:whereand is the location parameter, σ > 0 is the scale parameter, and is the shape parameter. When the shape parameter ξ is equal to 0, greater than 0, and lower than 0 [33], the GEV distribution is equivalent to Gumbel [37], Fréchet [38] and “reversed” Weibull distributions [39], respectively.

The Gumbel distribution, also named as the Extreme Value Type I distribution, has the following pdf and cdf: where μ is the location parameter, β > 0 is the scale parameter. Specially, when and the distribution becomes the standard Gumbel distribution. Generalizations of the Gumbel distribution, which are of flexible skewness and kurtosis due to the addition of one more shape parameter are widely used for extreme value data as they better fit data [40]. The distribution in (4.1) has been employed as a model for extreme values [41], [42]. The distribution has a light right tail, which declines exponentially, since its skewness and kurtosis coefficients are constant.

The Fréchet distribution, also known as the Extreme Value Type II distribution, has the following pdf and cdf, respectively:where α > 0 is the shape parameter and β > 0 is the scale parameter.

The Weibull distribution is known as the Extreme Value Type III distribution. The pdf and cdf of a Weibull random variable are shown as follows, respectively: where λ > 0 is the scale parameter and k > 0 is the shape parameter.

Further we show the application of the GEV model to the stock market close price using the weekly-return data that was calculated byThe results of fitting the GEV distribution to (weekly) block maxima data is presented in Fig. 6 and Table 1 that present the Quantile-quantile plot (QQ-plot), quantiles from a sample drawn from the fitted GEV pdf against the empirical data quantiles with 95% confidence bands. The maximum likelihood estimators of the GEV distribution are the values of the three parameters (μ, σ, ξ) that maximize the log -likelihood. The magnitude along with positive sign of ξ indicates the fat-tailness of the weekly-return data, which is consistent with the quantile plot.

Table 1

Parameter estimates for the GEV fitted model with maximum likelihood estimator. The 95% confidence intervals for each estimate are included.

	Location μ^	Scale σ^	Shape ξ^
Estimated parameter	606.3260	511.1713	0.1215
95% a lower bound of the confidence interval	576.43	486.42	0.05
95% an upper bound of the confidence interval	636.22	535.92	0.19

Based on the statistical analysis presented above (Fig. 6a), we see that the distribution of the weekly-return data can be described by a combination of different distributions. The density plot (Fig. 6b) having two humps validates the idea of a mixture of distributions.

How to choose a threshold

The classical approach for modeling extreme events is based on the GPD. It was proved [43] that if a threshold u is chosen and are observations above u, then the limiting distribution for excess over threshold is indeed GPD. In applications, the GPD is used as a tail approximation [44] of values exceeding the threshold u. The GPD is determined by scale and shape parameters σ > 0 and ξ, respectively, or in terms of threshold excess producing the following formulawhere . When ξ > 0, it takes the form of the ordinary Pareto distribution. This case is the most relevant for financial time series, since it is heavy-tailed. For security returns or high-frequency foreign exchange returns, the estimates of ξ are usually less than 0.5. When the GPD corresponds to the exponential distribution [45].

There are several properties of GPD [43], such as, ‘threshold stability’ property: if X is GPD and u > 0, then provided X > u is also GPD. Therefore, a Poisson process of exceedance times with generalized Pareto excess implies the classical extreme value distributions [46]. The above suggests that generalized Pareto distribution is a practical tool for statistical estimation of the extreme values, given a sufficiently high threshold. The rest of this chapter is devoted to a question about how high a threshold should be.

Graphical approaches to estimate threshold

One of the most common ways to determine a suitable threshold is to graphically inspect data. This approach [44] requires substantial expertise, that can be subjective and time consuming. In some cases, when dealing with several data sets, a uniform threshold may be proposed and kept fixed making the entire evaluation even more subjective.

The most common graphical tools are: mean excess plot [47], threshold stability plot [44], QQ-plot [48], Hill plot [49], return level plots [44], etc.

The mean excess plot is a tool widely used in the study of risk, insurance and extreme values. One use is in validating a generalized Pareto model for the excess distribution. The distribution of the excess over a threshold u for a random variable X with distribution function F is defined as

This excess distribution is the foundation for peaks over threshold modeling which fits appropriate distributions to data on excesses and widespread with many application in hydrology [50], [51], actuarial science [52], [53], survival analysis [54]. This modeling is based on the GPD that is suitable for describing properties of excesses. The mean excess (ME) function is one of the most common tools to determine a suitable threshold u. The ME function of a random variable X is defined asprovided which is also known as mean residual life function. As Ghosh [47] noted for a random variable X ≈ G , if and only if ξ < 1 and in this case, the ME function of X is linear in u:where if ξ ∈ [0, 1) and if ξ < 0. The linearity of the ME function characterizes the GPD class [47]. Davison and Smith [46] developed a simple graphical tool that checks data against a GPD model. Let be the order statistics of the data, then ME plot depicts the points where is the empirical ME function defined asIf the ME plot is close to linear for sufficiently large values of the threshold then there is no evidence against use of a GPD model. Another problem is to obtain a natural estimation of ξ. There are several methods to estimate such as: (i) Least squares [44], (ii) Maximum likelihood estimation [46], (iii) the Hill estimator [49], (iv) the Pickands estimator [55], (v) quantile-quantile plot (QQ-plot) [44], (vi) the moment estimator [56].

For example, the QQ plot depicts the points

where m < n and ξ > 0. In case ξ < 0, QQ plot is the plot of points where is an estimate of ξ based on m upper order statistics.

Recently, new graphical diagnostic tools have been introduced: a new multiple-threshold GP model with a piece-wise constant shape parameter [57]; plots measuring a surprise at different candidates to select threshold, using results of Bayesian statistics [58], [59]; structure of maximum likelihood estimators have been studied to develop diagnostic plots with more direct interpretability [60].

With this choice for a threshold, we get 507 exceedances with the empirical life becoming close to linear above this choice of the threshold (Fig. 7 ). Similarly, we find a threshold for negative returns. In this case, all computations were repeated for absolute values of negative returns (Fig. 8 ). With this choice for the threshold there are 462 exceedances. One can see that empirical MRL becomes almost linear above .

Fig. 7

Mean residual life plot for the S&P 500 positive returns. Solid jagged line is empirical MRL with approximate pointwise Wald 95% confidence intervals as dashed lines. The threshold u is estimated at 0.016. A vertical dashed line marks this threshold.

Fig. 8

Mean residual life plot for the S&P 500 negative returns. Solid jagged line is empirical MRL with approximate point-wise Wald 95% confidence intervals as dashed lines. The threshold u is estimated at 0.017. A vertical dashed line marks this threshold.

Ghosh and Resnick [47] noted that despite graphical diagnostic is a tool commonly accepted by practitioners, there are some problems associated with the methods mentioned above, such as: (i) an analyst needs to be convicted that ξ < 1 since for ξ ≥ 1 random sets are the limits for the normalized ME plot. Such random limits lead to wrong impressions. Certain methods described above work with ξ defined on specific intervals; (ii) in this case distributions are not close to GPD can mislead the mean excess diagnostics.

Based on the graphical approach, the threshold for negative return of S&P 500 was chosen at and for positive return at . Distribution of exceedances over respective thresholds are shown in the Fig. 9 .

Fig. 9

Distribution of exceedances normalized by thresholds for positive and for negative returns, respectively. Dotted lines represent Zipf’s Law (), dash dot lines represent Gaussian distribution () and dashed lines represent power law (). The curves are given for reference only.

The QQ-plots shown in the Fig. 10 demonstrate that with a few exceptions, exceedances follow the GPD.

Fig. 10

(a) Quantile-Quantile plot with maximum likelihood estimation for the negative threshold; (b) QQ-plot with maximum likelihood estimation for the positive threshold.

Automatic methods to estimate thresholds

As was mentioned above, the graphical approaches as well as rules of thumb can be highly objective, time consuming, and require certain professional background. Thus some authors have proposed automatic selection methods that can treat chunks of Big Data: a pragmatic automated, simple and computationally inexpensive threshold selection method based on the distribution of the difference of parameter estimates when the threshold is changed [61]: it was shown that better performance is demonstrated by graphical methods and Goodness of Fit metrics that rely on pre-asymptotic properties of the GPD [62] using weighted least squares to fit linear models in the traditional mean residue life plot; the recently developed stopping rule ForwardStop [63], which transforms the results of ordered, sequentially tested hypotheses to control the false discovery rate [64] that provides reasonable error control [58].

A particular interest has a method that suggests a way to determine threshold automatically without time consuming and rather subjective visual approaches based on L-moments of GPD that summarize probability distributions, perform estimation of parameters and hypothesis testing [58].

Probability weighted moments, defined by Greenwood [65], are precursors of L-moments. Sample probability weighted moments computed from data values arranged in increasing order, are given by

L-moments are certain linear combinations of probability weighted moments that have simple interpretations as measure of location, dispersion and shape of the data sample. The first few L-moments are defined by(the coefficients are those of the shifted Legendre polynomials).

The first L-moment is the sample mean, a measure of location. The second L-moment is (a multiple of) Gini’s mean difference statistic, a measure of the dispersion of the data values about their mean. By dividing the higher-order L-moments by the dispersion measure, we obtain the L-moment ratios, .

These are dimensionless quantities, independent of the units of measurement of the data. τ 3 is a measure of skewness and τ 4 is a measure of kurtosis - these are respectively the L-skewness and L-kurtosis. They take values between and . For random variable with GPD with ξ < 1, the particular relationship between L-skewness and L-kurtosis is defined as

Given a sample the Automatic L-moment Ratio Selection Method (ALRSM) works as follows [58]:

Define the set of candidate thresholds as 20 sample quantiles, starting at 25% by steps of 3.7%.

Compute the sample L-skewness and L-kurtosis for the excess over each candidate threshold and determine - the Eucledian distance:for with

The threshold after which the behavior of the tail of the underlying distribution can be considered approximately GPD is then automatically selected asthat is, the level above which the corresponding L-statistics fall closest to the curve.

Using the L-moments method, we computed thresholds for S&P 500 log-return depending on an observation period, 100 trading days and 400 trading days. The Fig. 11 indicates that a threshold not only time dependent, but also depends on the size of data set used. Once again, we cannot choose one value of the threshold that can be absolutely accurate.

Fig. 11

Value of the thresholds for positive and negative log-return based on L-moments. The solid black and blue lines correspond to negative and positive log return thresholds, respectively, and based on a window of 100 trading days. The dotted black and blue lines correspond to negative and positive log return thresholds, respectively, and based on a window of 400 trading days. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Rules of thumb to choose a threshold

As earlier noted, the threshold sequence is a function of the properties of the GPD provided that a population is in the domain of attraction of the GPD. In case a distribution function F is known, derivation of the threshold selection is possible, however, in practice if F is unknown then there is no general form for the threshold sequence [44].

Practitioners often use so-called rules of thumb, many of them have little to know theoretical justification, for instance, a fixed quantile rule [66]: the upper 10% rule or its derivative 5% rule; the square root rule [67]. There is a procedure that tries to find a region of stability among the estimates of the extreme value index [68]. This method depends on a tuning parameter, whose choice is further analysed in Neves and Alves [69]. Unfortunately, no theoretical analysis exists for this approach [70]. More comprehensive reviews of the threshold selection methods can be found in Scarrott and MacDonald [44].

Even though most of methods mentioned above have no theoretical justification for an exact value of a threshold, we can find an approximate location for a threshold given a data set R that has M values of log-return, .

Let R be a list of n consecutive values of log-return, . Split R into two parts: and ; then sort them in an increasing order such thatwhere such that and n is the width of observation or a window of observation.

Next, we compute medians of call them and . These values will be a lower bound for a positive threshold and an upper bound for a negative threshold, respectively.

At this step, an upper bound for a positive threshold x and a lower bound for a negative threshold x are estimated using the rule of thumb, namely, the fixed quantile rule with upper 10% for positive log-return values and lower 10% for negative log-return values. With this we haveThe indices u, l can be found as

A threshold for negative log-return values ranges from x to and a threshold for positive values is within and x, based on n observations from R.

Further, we shifted R to to estimate new values of thresholds based on previous n observations. The process repeated until we exhausted the entire data set R.

We chose a window of 300 days that was moving over the entire dataset producing a domain for threshold existence as shown in the Fig. 12 . It is clear from the Fig. 12 that certain values of thresholds cannot sustain the entire period and must be updated from time to time.

Fig. 12

Possible threshold changing ranges from 03/11/1981 to 12/31/2018 based on 300 preceding trading days. A green strip represents positive log-return and an orange strip shows the threshold domain for negative log-return values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Similarly, the window of 300 trading days, the Fig. 13 demonstrates how ranges of thresholds change as we move a window of 600 trading days across available data. Once again, some values of thresholds can exist for almost entire period, while other can exist a few months and then must be replaced with an updated value.

Fig. 13

Possible threshold changing ranges from 05/18/1982 to 12/31/2018 based on 600 preceding trading days. A green strip represents positive log-returns and an orange strip shows the threshold domain for negative log-return values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Based on Figs. 12 and 13 we can see that a choice of the threshold depends on a size of the data set used. Moreover, the rules of thumbs alike graphical approaches require practitioners’ involvement in studying data before making a final choice for thresholds.

Equipped with these results, we can compute a validity period τ(u) for a threshold u. Let n be the window of observation. By moving the window over the data set containing M values of log-returns we obtain four lists with elements in eachcontaining medians of positive and negative log-returns, upper 10% cut-offs and lower 10% cut-offs for positive and negative log-returns, respectively, as computed previously.

First, we compute the k threshold candidates u for positive log-return as

. In a similar fashion, we compute a set of threshold candidates for negative log-returns,

.

For each threshold candidate u found in the previous step we compute its validity duration, i.e. a number of days a candidate would fall within admissible ranges. Set . For if then for positive threshold candidates, similarly, if then for negative threshold candidates.

The probability that a threshold candidate will be changed on any given day is .

The Fig. 14 shows probabilities that a particular choice of the threshold can be changed on any day depending on the size of the dataset used, which is important to know especially when new data becomes available and is included for consideration. The more data we use to estimate a value of a threshold, the more likely the threshold will stay unchanged, however, as we add data, the threshold should be reconsidered. It also brings another issue: in many cases, statistical analysis is performed on a historical dataset that does not reflect a phenomenon we study at the present time.

Fig. 14

The probability that a threshold of the log-return values will be changed on any given day calculated over the different data windows ranging from 25 to 6000 trading days.

Statistics of extreme events under threshold uncertainty

The statistics of extreme events under threshold uncertainty can be described with the help of a simple model, in which the log-return of the index and the value of threshold are treated as random variables that yet can change inconsistently. The model that we are going to adopt and modify had been put forward by us for the first time in Floriani et al. [71] to describe the behavior of systems close to a threshold of instability and used later in Volchenkov [72], 73] to model survival under uncertainty and the events of mass extinction. In the model of extreme events under threshold uncertainty, the current value of the log-return is quantified by a random number x ∈ [0, 1] drawn accordingly some probability distribution function . The threshold value that might change any time once the new data become available is another random number y ∈ [0, 1], which is drawn from another probability distribution function, We assume that the rate of daily variations of the log-return values is greater than or equal to that of the threshold values, ultimately determining whether the current log-return value is extreme or not. In fact, it is the relative rate of random updates of x and y described in our model by the probability of inconsistency η ≥ 0, that actually determines the statistics of extreme events.

At time the value of log-return x is chosen with respect to the probability distribution function F, and the value of the threshold y is chosen with respect to the probability distribution function G. If y ≥ x, the event is regular and the process keeps going to time . At time t ≥ 1, either, with probability η ≥ 0, the value of log-return x is drawn anew from the probability distribution function F, but the threshold keeps the value y it had at time or with probability the value of log-return x is updated anew from the probability distribution function F, and the level of supply y is updated either with respect to the probability distribution function G. As long as the value of threshold is not exceeded (x ≤ y), the event is classified as regular, but the event is extreme whenever x > y.

The value of probability η > 0 can be interpreted as the reciprocal characteristic time interval, during which the threshold level remains unchanged, and vice versa the probability that a threshold of the log-return values will stay unchanged on any given day can be calculated as the inverse of the time interval during which the threshold value stays put. The probability that a threshold of the log-return values will stay unchanged on any given day calculated over the different data windows ranging from 25 to 6000 trading days is shown in Fig. 14.

We are interested in the probability distribution P(t) of the interval duration t between the sequential extreme events for some probability distribution functions F and G and a given value of the probability η ≥ 0. A straightforward computation [71], [72], shows that, independently upon the value of η, the initial probability of choosing the level of demand below the supply level (to start the subsistence process) is The general formulas for P(t) can be found in Floriani et al. [71], Volchenkov [72].

When the values of log-return and the threshold value are updated coherently the resulting probability function,decays exponentially fast with t, for any choice of the probability distribution functions F and G. In particular, if the threshold level and log-returns are drawn uniformly at random, over the interval [0,1], the occurrence of an extreme event is statistically equivalent to simple flipping a fair coin, for which head and tail come up equiprobably,

On the contrary, when the threshold level is kept unchanged, the statistics of intervals between the sequential extreme events,decays asymptotically algebraically as t ≫ 1 [71], [72]. For example, in the special case of uniformly random updates of the threshold and log-return values, the probability function (4.10) decays algebraically as

For a general family of invariant measures of a map of the interval [0,1] with a fixed neutral point defined by the probability distributions F and G, absolutely continuous with respect to the Lebesgue measure, i.e., Eq. (4.10) gives the probability function that exhibits the power law asymptotic decay for t >  > 1 [71], [72], [73]:The asymptotic decay of (4.13) seems to be algebraic,for for any choice of the distributions F and G although it is mainly the character of the probability function G that determines the rate of decay of with time. In the limiting case when the support of the probability distribution G(x) determining the choice of the supply level is concentrated close to i.e., is zero everywhere in the interval [0,1], except for a small interval of length ε up to 1, the Zipf power lawasymptote ε > 0, follows directly from (4.13) [71], [72], [73]. A possible modeling function for such a bountiful probability distribution, forming a thin spike as x → 1, can be chosen in the form,with the probability density in the interval [0,1[,The straight line shown in Fig. 15 represents the hyperbolic decay of time intervals between the extreme events.

Fig. 15

The statistics of time intervals (in days) between the sequential extreme events for the fixed threshold values, and for positive and negative fluctuations of the log-return respectively. The solid line corresponding to the asymptotic quadratic decay (4.11) is given for reference.

Defining emergency scales by thresholds uncertainty.

Once we assume statistics for log- returns, or ranges for threshold values are defined using different methods shown in Section 4.1, we introduce uncertainty to the threshold value that depends on the method itself and the amount of data we use. We measure uncertainty to justify an emergency scale to represent the extreme events.

Using the rule of thumb (Section 4.2.3), we determined the ranges for the threshold values depending on the window of observation (Fig. 12) for both positive and negative log return values. For each admissible positive choice of a threshold u, we determine the probabilities η(u) that the threshold u can be changed on any given day (Fig. 14), degrees of uncertainty have been assessed by the means of the Shannon’s entropy for each threshold value and the window of observation [73] (5.1). The case with the negative values of the threshold is analogous.where η(u) is the probability that the threshold will be changed on any given day (Fig. 14).

We observed the Red Queen State and three emergency scales can be readily interpreted.

If amount of data is not sufficient (a solid line corresponding to the 18-day window of obervation in Fig. 16), then the uncertainty curve forms a skewed profile attaining a single maximum of 0.69295 for the threshold value 0.008272. In this situation, an observer’s perception of the events reminds Red Queen from “Through The Looking-Glass and What Alice Found There” by Lewis Carroll [21] who said “When you say hill, I could show you hills, in comparison with which you’d call that a valley”. Our understanding of the events whether they are extreme or not is very limited and uncertainty is blurry. As events become more severe, our uncertainty that the events are extreme decreases. The observer realizes that the events are extreme, but a precise point at which the events turn to be severe cannot be determined. This case is called the Red Queen State.

As the window of observation becomes larger, 25 days, for instance, the uncertainty curve exhibits two maxima indicating that the amount of data is sufficient. This happens because H(η) attains its maximum for and the η curve admits a value 1/2 twice (Fig. 14). As we further extend the window the curve torrents into sharp peaks (Fig. 16). The latter ones clearly separate the threshold values into three regions: three levels of emergency. The locations of peaks of the curves are summarized in Table 2 . The location of two peaks is not sensitive to the window of observation.

Table 2

Location of extrema points of the uncertainty curves from Fig. 16.

Window	First maximum		Second maximum
	Threshold	Uncertainty	Threshold	Uncertainty
18 days	0.00827	0.69296	–	–
50 days	0.00461	0.67411	0.01532	0.69290
2000 days	0.00534	0.69315	0.01672	0.69034
4000 days	0.00533	0.68217	0.01580	0.69314

Emergency scales for S&P 500.

Subcritical. In the region I, the threshold values are small to raise concern about extreme events. Then we can observe a spike with the degree of uncertainty attaining its first maximum. This extremum indicates a transition to the next kind of uncertainty. For the window of 2000 days, the region I lies in the interval [0, 0.00534) (see Fig. 16, Table 2).

Critical. In the region II, the interval [0.00534, 0.01672) for the 2000-day window, uncertainty is conceptualized with a question whether a magnitude of the event is already critical, extreme or not yet. Further, we see another jump of uncertainty reaching 0.69034. At this point we are certain that events are not regular anymore.

Extreme. In the case of the window of 2000 days, the interval [0.01672, ∞) constitutes the region III. We consider all events in this region extreme with our uncertainty decreasing as threshold values increase.

With the analysis presented above, we define an emergency scale of three levels based on the three regions of the threshold values corresponding to three peaks of the uncertainty curve. This emergency scale is not sensitive to the size of the window of sufficient amount of data considered.

Conclusion

The S&P500 times series in the period from January 2, 1980 till December 31, 2018 exhibits an asymmetrical skewness of the distribution with the right and left power law tails. Multifractal detrended fluctuation analysis of log return time series for S&P 500 index reveals a scale invariant structure for the fluctuations on both small- and large scale magnitudes, as well as its short- and long-range dependence on different time scales. Moreover, the segments with small fluctuations have a random walk like structure whereas segments with large fluctuations have a noise like structure.

We have reviewed different methods of threshold selection and studied the extreme events presented in the time series using different statistical approaches. We found that the distribution of the weekly-return data can be described by a combination of different distributions. Based on a graphical approach for threshold selection, we chose separate thresholds for the positive and negative values of the log return, 0.016 and respectively. With this choice, we registered 507 instances of extreme events corresponding to raise of market and 462 extreme events related to market declines. With a few exceptions, exceedances over (under for negative log return values) the threshold follow the GPD.

The rule of thumb showed that a threshold value depends on the width of observation window, and the threshold can change at any moment, once new data become available. Uncertainty of the threshold values can be determined by the probability of changing the threshold on any given day.

The moment we assume statistics of distributions or the dataset is fixed, it leads to uncertainty of the threshold value which can be resolved by the emergency scales rigid to variation on the size of the dataset.

We suggested a statistical model that describes registration frequency of extreme events under threshold uncertainty. Our model fits well the statistics of occurrence of the extreme values in the S&P 500 time series.

Declaration of Competing Interest

The authors declare no conflict of interests.