Rare and extreme events: the case of COVID-19 pandemic

Complex systems have characteristics that give rise to the emergence of rare and extreme events. This paper addresses an example of such type of crisis, namely the spread of the new Coronavirus disease 2019 (COVID-19). The study deals with the statistical comparison and visualization of country-based real-data for the period December 31, 2019, up to April 12, 2020, and does not intend to address the medical treatment of the disease. Two distinct approaches are considered, the description of the number of infected people across time by means of heuristic models fitting the real-world data, and the comparison of countries based on hierarchical clustering and multidimensional scaling. The computational and mathematical modeling lead to the emergence of patterns, highlighting similarities and differences between the countries, pointing toward the main characteristics of the complex dynamics.


Introduction
Many complex systems generate outputs that are characterized by a frequency-size power law behavior over several orders of magnitude [1,2]. The power laws have been associated with scale-invariance, self-similarity, and fractality and are consistent with self-organized criticallity, that is a process in which a system, by itself, converges to a state characterized by a coherent global pattern, created by local interactions between low-level elements [3,4].
The power laws are characterized by heavy-tails, giving non-negligible probability to large events. However, some extreme events, labeled 'dragon kings', while predictable, cannot be foreseen by the extrapolation of power law distributions [5,6]. 'Dragon kings' may be associated with positive feedback, bifurcations, and regime changes in out-of-equilibrium complex systems. 'Dragon kings' are often discussed in contrast with 'black swans', which denote unpredictable catastrophic rare events [7]. These outliers are pervasive in many areas, namely economy, finance, earth sciences, and biology.
The recent Coronavirus disease 2019 (COVID-19) outbreak is an example of an extreme event. We must highlight that the occurrence of a rare event and the actual description of its evolution are, however, distinct matters. This paper attempts to understand the dynamics of the spreading across different countries of COVID-19, but not the prediction of its outbreak or conclusion.
The first case of COVID-19 [8][9][10] was officially reported in China on December 31, 2019, in Wuhan of Hubei province. At an early stage, the Chinese authorities seemed not to give importance to the problem. However, with the rapid emergence of new cases, the attitude changed dramatically. The Chinese government took a series of strong measures to contain the disease and gave the world an example of commitment and effectiveness. The growth rate of new cases of COVID-19 in China has slowed significantly, and the situation appears to be under control at the moment of writing this paper [11,12].
In the meantime, new cases have been emerging in many countries. In particular, the rapid evolution in Iran, South Korea, Italy and Spain became the most dramatic cases. COVID-19 was gradually reaching all continents, with cases confirmed all over the world, while having 'alarming levels of inaction' by some countries, in the words of the director of the World Health Organization (WHO).
More recently, by March 11, 2020, WHO officially declared the COVID-19 a global pandemic, just when the number of known cases reached approximately 121,000 and caused 4300 deaths, and after the cases outside of China spread by a factor of 13 and the number of countries affected tripled in just two weeks.
Panic begins to spread in some populations [13], fueled by the massive and speculative news broadcast by the media and social networks. World governments have taken drastic measures, such as closing schools, entertainment venues and restaurants, and the movement of people has slowed dramatically, in parallel with thousands of canceled flights [14]. The tourism sector and commercial flights are the most affected and the world seems to be heading towards an economic recession, with the GDP of some countries being able to drop double-digit figures.
No one can say whether the measures being taken are sufficient [15], or what the evolution of the pandemic will be, but this appears to be the public health crisis of a generation. However, we cannot forget that, for example, the H1N1 flu of 2009 caused between 151,000 and 575,000 deaths worldwide [16]. The COVID-19 has still a long way to go to reach the H1N1 levels. The world faced other flu pandemic crises in the past [17] and the scientific knowledge has never been so well prepared as today to give appropriate answers to health crises.
The analysis of the evolution of the confirmed cases versus time has considerable interest from the point of view of delivering good information to health organizations and to the general public. Several statistics have been presented, adopting different forms for organizing and visualizing the data. However, a comprehensive representation of the COVID-19 spreading dynamics across different countries is still missing.
In epidemiology, mathematical modeling plays an important role in understanding the mechanisms that govern the transmission of contagious diseases. The work by Kermack and McKendrick [18] formulated the general theory of susceptible-infected-recovered (SIR). A SIR model involves a system of coupled equations relating the numbers of susceptible, infected and recovered people over time, and computes the theoretical number of infected people in a closed population. Many variations of the original SIR model were proposed, such as the susceptible-infectious-susceptible (SIS) and the susceptible-exposed-infectious-recovered (SEIR), based on ordinary [19], stochastic [20] and fractional order [21,22] differential equations. These recent versions were adopted for studying the spread of distinct infectious diseases [23,24], including the COVID-19 [25]. However, the main concern of such approach is the model validation, since it requires to compare the results with real data. Contrary to modeldriven, data-driven approaches rely on data series for deriving adequate fitting functions. These heuristic models describe well one stage of the epidemic, but fail when the disease evolves toward a different phase. We must also note that the heuristic model is useless in the initial epidemic phase, due to insufficient data [26,27].
The paper addresses the statistical comparison and visualization of COVID-19 country reported cases in the period December 31, 2019, up to April 12, 2020. The study does not aim to be a contribution tailored for medical treatment or prevention of the disease. Therefore, in a first phase, we adopt a nonlinear least-squares technique to determine possible candidate heuristic models for describing the data regarding COVID-19 infections. In a second phase, we use distinct metrics for processing the data both in the time and frequency domains. The information is visualized using hierarchical clustering (HC) and multidimensional scaling (MDS) for comparing the COVID-19 evolution in the different countries. The HC and MDS generate loci of Fig. 1 Geographic map of the COVID-19 spread for 165 countries. The color map is proportional to the number of days elapsed since the occurrence of the first case for each country for the period of time τ points in 2-and 3-dimensional spaces representing the number of infections for each country. The positioning and the patterns formed by the points lead to direct interpretations of the results. The study is data-driven, and the models are applicable only at some stages of the outbreak and when enough data points are available.
In this line of thought, the paper is organized as follows. In Sect. 2 we introduce the dataset adopted in the follow-up. In Sect. 3 we analyze the data by means of regression modeling. In Sect. 4 we compare and visualize the COVID-19 spreading data in various countries. In Sect. 5, we discuss the possibility of foreseeing the future evolution. Finally, in Sect. 6 we discuss the results and summarize the main conclusions.

The dataset
The COVID-19 data are made available by the European Centre for Disease Prevention and Control (https:// www.ecdc.europa.eu/en). The dataset is provided in Excel format, containing the number of infected and the number of deaths for each country, on a daily basis. Data for the period from December 31, 2019, up to April 12, 2020, were collected for analysis. This period of time will be denoted as τ henceforth. Figure 1 depicts a geographic map where the colormap is proportional to the number of days elapsed since the occurrence of the first case in each country. We verify that the COVID-19 is particularly severe in the northern hemisphere and, thus, seems not to follow the same pattern of other serious diseases that affected mainly the underdeveloped countries. Therefore, some possible synchronization between countries, or, even, the emergence of new waves of spread in the future are still unclear and techniques such as the Kuramoto model [28] for assessing that hypothesis should be considered.
Let x i (t) denote the time series of confirmed COVID-19 daily infections for the ith country, i = 1, . . . , M, where t = 1, . . . , T represents time with one day resolution, within the time interval τ . Therefore, the signals x i (t) evolve in discrete times, t, and can be interpreted as one manifestation of a complex system.
For the sake of statistical significance and accuracy of the mathematical tools used for processing the data, we just consider the countries with time and we approximate the histogram of r i (t) by a symmetric α-stable distribution [29]. We recall that a probability distribution (PD) (and the corresponding random variable ξ ) is said to be 'stable' if a linear combination of 2 independent random variables with such PD has also an identi-cal PD, up to the scale and location parameters, c and μ, respectively [29,30]. A given family of stable distributions is often called Lévy alpha-stable distribution, after Paul Lévy [31]. The Lévy, Gaussian and Cauchy PD of a random variable ξ are particular cases of the α-stable distribution family with the parameter value α = 1 2 , 1 and 2, respectively [32]. The α-stable distribution is a four parameter family of distributions and is (usually) denoted by S(α, β, c, μ). The first parameter α is of particular relevance and describes the tail of the distribution. We have α ∈ (0, 2] for the stability (or characteristic exponent), β ∈ [−1, 1] representing the skewness, c ∈ (0, ∞) standing for the scale, and μ ∈ (−∞, +∞) for location parameters. With exception  Fig. 2 The histogram of the the log returns r i , i = 14, and αstable approximation, with tail characteristic exponent α = 0.26, for China during the period τ of the cases when α ≤ 1 and β = ±1, we have that for α < 2 the asymptotic behavior is described by [33,34]: where Γ denotes the Gamma function. We verify the presence of the so-called heavy or fat tails that cause the variance to be infinite for α < 2. Figure 2 depicts the histogram and the α-stable approximation for China, r 14 , in the period τ . Twelve bins were considered for having statistical significance. In this case, we have approximately f (r ) ∼ 1/|r | 1.26 , that is, α = 0.26 corresponding to an extremely small value, which entails a huge probability for extreme values of the return. The alternative of an asymmetrical PD was tested, but the resulting improvement was minor and by consequence merely the symmetrical version is depicted for the sake of simplifying.

Regression models for describing the spread of COVID-19
Let y i (t) represent the time series of cumulative number of infections of x i (t), that is y i (t) = t n=1 x i (n). We adopt the nonlinear least-squares [35,36] to examine the behavior of y i (t) for a variety of functions. We selected the 'Logistic' and 'Richards' models: for approximating the data of China (i = 14) and Italy (i = 36), respectively, where a, b, c, d ∈ R are parameters adjusted by means of a nonlinear least square fit numerical algorithm. These models were selected from a large number of heuristic functions simply because they (i) adjust adequately to real data, and (ii) involve a limited number of parameters. Therefore, no special biological meaning was intended when using such functions. Figure 3 illustrates the data time series, y 14 (t) and y 36 (t) and the corresponding approximationsŷ 14 We verify a good fit in both cases, with coefficient of determination R 2 = 0.99, but a single model with limited number of parameters is not able to fit well the time seriesŷ i (t) for all countries. Obviously, we can adopt other models involving a larger number of parameters for achieving a better fitting to a givenŷ i (t). Nonetheless, only analytical expressions requiring a limited set of parameters are of relevance [37]; otherwise, their comparison and interpretation becomes unclear. On the other hand, the use of distinct models for different countries lack generality when comparing results.

Global comparison of the COVID-19 spreading
We now analyze the COVID-19 spreading data of M = 79 countries both in the time and frequency domains. In the time domain, we compare the pair (i, j) of countries by the corresponding time series of the cumulative In the frequency domain, the pairs of countries are compared by the daily number of infec- . . , ω K , F{·} denotes the Fourier transform, ω represents the angular frequency and ı = √ −1. We adopt the Canberra distance to measure the dissimilarity between pairs (i, j) for the time and frequency domains: Real data Fitting model that is, distances based on the variables y i (t) = t n=1 x i (n) and X i (ıω), respectively, where Re{·} and Im{·} denote the real and imaginary parts. The Canberra distance has the relevant property of being relatively insensitive to the simultaneous presence of large and small values.
Obviously, other distances are possible [38] and several of them were also tested. However, further distances are not included herein for sake of parsimony, since d t,y C i j and d f,X C i j illustrate adequately the proposed concepts.

Hierarchichal clustering and visualization of COVID-19
For visualizing the relationships between the 79 countries, we first adopt the HC computational approach. The HC is a technique that groups similar objects [39]. Given M objects in a q-dimensional real-valued space and a dissimilarity metric, . . , M, of object to object dissimilarities is determined [40]. The HC generates a structure of object clusters, using Δ as input, that is represented graphically either by a hierarchical tree or a dendrogram. We have two alternatives to generate a hierarchy of clusters, namely the agglomerative and divisive clustering iterative techniques. In the first, each object starts in its own cluster and the successive iterations merge the pair of most similar clusters until there is a single cluster. In the second technique all objects start in one cluster and, during the iterations, the 'outsiders' are removed from the least cohesive cluster, until each object is in a separate cluster. In both cases the HC requires a linkage criterion, that is a function of the distances between pairs of items, for quantifying the dissimilarity between clusters. Metrics such as the maximum, minimum and average linkages are often used. The distance d (x R , x S ) between two objects x R ∈ R and x S ∈ S, in the clusters R and S, respectively, can be assessed by means of several metrics such as the average-linkage given by [41]: For assessing the clustering quality, the index cc is mostly used [42]. Let us consider that the original object X i is described by a HC representation T i . Additionally, let x(i, j) and t (i, j) stand for the the distances between the X i and X j original observations and the HC points T i and T j , respectively. If we havex = av (x (i, j)) andt = av (t (i, j)), where av(·) denotes average, then cc is given by: The closer the value of cc is to 1, the better the clustering reflects the original data. The results are represented in a Shepard chart that compares the original and the cophenetic distances. The closer to the 45 degree line the points, the better the obtained clustering. For example, in MATLAB, the cophenetic correlation coefficient can be obtained by means of the command cophenet.
Herein, the agglomerative clustering and averagelinkage method are adopted for visualizing the two resulting matrices of item-to-item distances based on Eqs. (5) and (6), respectively [43].
In Fig. 4, we see a clear position of China followed by the sub-cluster formed by Iran, Italy, Spain, Korea, United States, France and Germany. On the other hand, the tree based on the frequency response gives more importance to a sub-cluster formed by the United States and China, followed by a second group including United Kingdom, Italy, Germany, France, Spain, and Iran. Moreover, the second tree separates better those countries with a smaller impact from the virus spread. In Fig. 4, the countries with a smaller (larger) number of occurrences and a smaller (larger) number of days since the first case are to the left (right). In Fig. 5, the distribution for smaller (medium/larger) values is located on the right (left/middle bottom) sides. Figure 6 represents the Shepard plot for assessing the HC tree for the 74 countries and the item-to-item dissimilarities d t,y C i j . The chart reflects an accurate clustering of the original data. For the index d f,X C i j , the Shepard plot is identical to the one in Fig. 6 and, therefore, is not presented.

Multidimensional scaling and visualization of the COVID-19 dataset
The MDS is a computational technique for clustering and visualizing multidimensional data [44]. As for the  We have several variants of the MDS, such as the metric, non-metric and generalized MDS. In the case of the metric MDS, the iterative algorithm minimizes the stress cost function S. We can have for example the residual sum of squares: The Sammon criterion can be also adopted The MATLAB command cmdscale and stress criterion Sammon were adopted. The interpretation of the MDS locus is based on the patterns of points. Similar (dissimilar) objects are represented by points that are close to (far from) each other. Therefore, the information retrieval is not based on the point coordinates, nor the shape of the clusters. This means that it is possible to magnify, translate and rotate the MDS locus. The axes of the MDS plot have neither units nor a special physical meaning. The quality of the MDS can be assessed through the stress and Shepard diagrams. The stress chart represents S versus d. The plot is monotonically decreasing and choosing a given value of d is a compromise between obtaining low values of S and d. The values d = 2 or d = 3 are usually adopted, because they allow a direct representation. The Shepard diagram compares φ i j and δ i j for a given value of d. A narrow scatter represents a good fit between φ i j and δ i j . Figures 7 and 8 depict the 3D MDS maps for d t,y C i j and d f,X C i j , respectively, for the period τ . As before, the size of the dots is proportional to the logarithm of the number of infections (i.e., ln[y i (T )]) and the color is proportional to the time between the first reported case in each country and T . The clusters are the same obtained with the HC, however, for their clear visualization we need to rotate the 3D maps.
As for the HC we verify that the time domain approach makes a better distinction of the countries with a larger number of infections, while the MDS based on the frequency domain has a more eclectic distribution, that is with a larger dispersion, and leaves some room for distinguishing the countries with a smaller number of infections. In both figures, the countries to the left (right) have a smaller (larger) number of infections and more (less) time elapsed since their first case in each country.   when varying time [40,[45][46][47][48]. Procrustes is a statistical method that takes a collection of shapes and transforms them (using translation, rotation, and amplification/reduction of size) for maximum superposition. The comparison is performed for a shorter period of time τ from t = 40 up to t = T , so that there is significant non null data for the analysis.
Let y c (t), c = 1, . . . , 16, denote the data time series representative of the countries in C for periods of time starting at t = 40 and increasing up to k ≤ T . The individual 3D MDS maps (one for each value of k) are generated using the item-to-item dissimilarities d t,y C i j and d f,X C i j , and processed with Procrustes. We must note that we are now using T − 39 matrices Δ of dimension 16 × 16. Therefore, the MDS locus produced for each k is not a magnification of the previous charts. For each index, the collection of MDS maps yields one global chart (Figs. 10, 11, respectively) that represents world recent evolution of the COVID-19. We verify different behaviors of the time and frequency domains MDS loci being, apparently, slightly more representative the first. In Figs. 12 and 13, we redraw the MDS charts of Figs. 7 and 8 (for the initial period of time τ ) with a connection between those countries that lie close to each other in the MDS locus [49]. Therefore, the lines do not represent clusters, and, instead they indicate that in the near future the evolution of a given country will probably be similar (in the sense of the adopted distances) to the neighboring countries. We observe also some discontinuities in the lines. Again they do not represent different clusters. The discontinuities simply mean that some neighbor is closer than the other and, therefore, it is likely that its evolution is closer.

Is it possible to foresee?
It is written 'The future belongs to God, and it is only he who reveals it, under extraordinary circumstances' [50]. To the authors best knowledge, the HC and MDS techniques are not designed to make predictions. Indeed, they allow a better interpretation of the past and present. Nonetheless, we can take advantage of the MDS computational visualization to trace some similarities between the items represented in the loci. From Figs. 10 and 11, where time is a parametric variable, we verify that we do not obtain our intuitive feeling of time as a smooth and continuous variable embedded and synchronizing all events. In fact, we see that the time instant for the beginning and end of the trajectories vary considerably from country to country. Therefore, if we adopt a critical view we can interpret that the 1-dimensional time continuum (if exists) is not adequately represented by the two technique combination (i.e., MDS and Procrustes). However, it was already noticed in previous studies [51,52], addressing a distinct phenomenon and applying only MDS, that datasets under the influence of social and human factors phenomena exhibit a relativistic behavior and eventually different velocities. Let us name the 'relativistic time' to make it distinct from the standard notion of constant speed physical time.
Let us consider the Canberra and Lorentzian distances [38] to measure the dissimilarity between the pairsx i (t) andx j (t): wherex represents the ith vector of N consecutive values of x(t) obtained from non-overlapping time windows. Therefore, the time series is subdivided into identical periods of time giving rise to R = T /N windows, where · denotes the integer part of the argument. Obviously, the larger the width of the window the better the filtering, but the weaker the notion of time 'instant'. Now, the MDS has for input a matrix Δ of dimension R × R and produces a single locus that compares the set of N -dimensional R vectors, corresponding to the different time windows. The MDS with MATLAB command cmdscale was adopted. We verify again the emergence of areas involving a large variability, coincident with large transients, and zones with a smoother evolution, corresponding to a continuous dynamics. In both cases we observe clearly four phases: (i) an initial transient, (ii) a fast progress up to a peak, (iii) a return back (but not exactly to the initial state), and (iv) the emergence of a new, unclear, second wave. As usual with the MDS technique the distinct measures highlight different aspects but the overall conclusions are similar.
For estimating future outcomes we now propose a technique embedding the trendline and the MDS datadriven techniques for estimating the evolution. In what follows we adopt the case of Italy (i = 36) as our test bench and we shall tackle the values of the number of daily infections x i (t). The main idea is to fit a set of trendlines to the available data and, based on them, to extrapolate the future behavior. In a second phase we adopt MDS to compare the real-world data and the estimations provided by the trendlines.
We consider four models, namely the 'Hoerl', 'Reciprocal quadratic', 'Gaussian' and 'Vapor' given by: respectively, where a, b, c ∈ R are parameters and t ∈ τ . We emphasize again that these models have no specific meaning and are just some functions that fit adequately the available data. In fact, the proposed heuristic models follows the common sense that number of infections will diminish in the future. However, these trendlines are just for estimating the near future and we shall consider M 1 and M 3 as representing a 'optimistic' scenarios, while M 2 and M 4 stand for 'pessimistic' future outcomes. Moreover, models M 1 , M 2 and M 4 have an asymmetric evolution about the peak, while M 3 con-  Figure 15 represents the curve fitting for the period τ and an estimation up to an extended period of time τ = 80 days.
For assessing the quality of the estimations we compare the real data and the results provided by M j , j = 1, . . . , 4, for Italy using the Canberra and the Lorentzian distances defined previously in (11)- (12). As before, the heuristic models cover an extended period of τ = 80 days for providing an estimation of the near future. We adopt the stress S for assessing the conformity between the data and the trendlines including their estimation. The MATLAB command cmdscale and stress criterion Sammon were adopted. Figure 16 depicts the stress yielded by the MDS loci generated by the Canberra and the Lorentzian and the models M j , j = 1, . . . , 4 for the case of Italy.
We verify that for d = 3 the models M 2 and M 3 are superior when considering the Canberra distance. However, the models M 1 and M 3 are superior in the perspective of the Lorenzian distance. Figures 17 and 18 show the resulting MDS loci for the Canberra and the Lorenzian distances for the two best models in each case.
The Lorentzian distance is more sensitive and gives more space to the artifact represented by the value x 36 (t = 23) = 90 (the neighbor values are x 36 (t = 22) = 2547 and x 36 (t = 24) = 6230), but, on the other hand, shows more clearly the recent final phase. The Gaussian model, M 3 seems a good compromise between the two distances and produces a trajectory consistent with the data series without exhibiting acceleration/deceleration time periods distinct from those available at the time of writing the paper. Nonetheless, the authors highlight that the COVID-19 evolution has an underlying plethora of phenomena going from cultural and economical up to political and geographical issues. As someone said, 'Prediction is very difficult, especially if it's about the future' [53].

Conclusions
This paper investigated an example of an extreme event, namely the dynamics of the COVID-19 spreading. Two approaches were considered for the period from December 31, 2019 up to April 12, 2020. In a first phase, heuristic models were used to fit the time series of the number of infections verified in a set of 79 countries. In a second phase, two metrics were used for comparing the countries data both in the time and frequency domains, and the HC and MDS techniques were adopted for clustering and visualization. The time evolution was also considered for a group of countries exhibiting a more dramatic spread of the COVID-19. The combination of Procrustes and MDS showed that besides the number of infections, the dynamic characteristics play an important role that is not evident in standard representations. In fact, the computational and mathematical modeling lead to the emergence of patterns both highlighting the main clusters and the similarities or dissimilarities between them. Given the potential of the techniques discussed here we can think of their application in several research directions such as the subdivision of data sets according with distinct criteria such as the age, or the geographical origin of the patients. Additionally, analysis considering a large set of influential factors besides merely the casualties can be tried, such as infected people staying in general or in intensive care in hospital, or recovered based on some type of treatment. Therefore, if sufficient and assertive information is collected, then research can follow the aforementioned computational methods to unravel space and time nonlinear dynamics embedded the data.