Introduction

Throughout human history, several pandemics have significantly reduced world population (^{Hays, 2005}; ^{Huremović, 2019}). For instance, a century ago the world was afflicted by the Spanish flu epidemic (H1N1), which began and peaked in 1918 and persisted until 1920 (^{Hays, 2005}; ^{Huremović, 2019}). This pandemic was one of the first in which data and statistics (from 48 countries) allowed health authorities to estimate that 40.1 million people died, i.e., 2.1% of the world's population at that time (^{Barro, et al., 2020}). On the other hand, the COVID-19 coronavirus disease (SARS-Cov-2) pandemic is the first with data available in real-time (^{Dong, et al., 2020}) on the time series of the total confirmed cases, confirmed cases per day, total deaths, deaths per day, active cases, total test results, etc., for each country (^{Worldometers, 2021}). The availability of time series for confirmed cases and deaths (^{GitHub, 2021}) from January 1, 2020, to September 15, 2021, allowed us to detect some properties also found in a wide variety of complex systems.

Recently, different aspects of the spread of the COVID-19 pandemic have been studied from the point of view of complex systems (^{Holovatch, et al., 2017}). Just to mention a few of these studies, researchers analyzed the spatial evolution of the Covid-19 pandemic around the world based on the coronavirus cases time series per country using a particular type of neural network called self-organizing maps to show that it has behaved similarly in certain groups of countries (^{Melin, et al., 2020}). In another study, time series of daily infections per country were analyzed using a methodology developed to study network problems showing that there are groups of countries with differentiated contagion dynamics (^{Álvarez, et al., 2020}). Moreover, by using a toy model for social networks and the so-called small-world network, researchers explained why most COVID-19 infection curves are linear (^{Thurner, et al., 2020}). On the other hand, three clusters of countries and territories with similar profiles of standardized COVID-19 time dynamics were identified by using a machine learning technique (^{António, et al., 2021}). Besides, times series of active and cumulative cases from three countries were analyzed in the context of a continuous effective model based on a modified diffusion equation under Lifshitz scaling to describe some aspects of the SARS-Cov-2 virus diffusion and to reproduce the spreading of the pandemic between two cities of Brazil (^{Anacleto, et al., 2020}). Additionally, based on confirmed cases in the US distributed across more than 2,000 counties, researchers found that the right tail of this distribution exhibited a power law with Pareto exponent close to one (^{Beare & Toda, 2020}). Furthermore, the formalism of the propagation of a wave packet in a dispersive medium was used to describe the spreading of the COVID-19 pandemic in Italy (^{Feoli, et al., 2020}). Finally, a mathematical framework was developed to study the economic impact of pandemics like the COVID-19 by integrating epidemiological dynamics with a kinetic model of wealth exchange (^{Dimarco, et al., 2020}).

As a result of the first stage into a research program in progress about the spreading of the COVID-19 pandemic based on the analysis of empiric data from January 2020 to September 2021, it is evident that the COVID-19 pandemic has spread around the world exhibiting complex system properties such as lognormal laws, temporal fluctuation scaling, and time correlation. In contrast with ^{Beare & Toda's results (2020)}, where the cumulative confirmed cases distributed among the US counties exhibited a power law, in the present study, the daily cumulative numbers of confirmed cases and deaths distributed among the countries of the world exhibited lognormal laws. In this respect, it should be mentioned that lognormal laws emerge in a wide variety of complex systems in different science fields (Limpert, et al., 1991) because of the presence of stochastic multiplicative processes in them (^{Mitzenmacher, 2004}). It is also clear in this paper that temporal fluctuation scaling (TFS) behaves as a power-law relation between the standard deviation and the mean when these are calculated for the daily cumulative number of confirmed cases and deaths distributed among countries in the world. With respect to this property, from a general perspective, fluctuation scaling affirms the existence of a power-law relation between the standard deviation and the mean of a time series (^{Eisler & Kertész, 2006}; Eisler, et al., 2008; ^{Fronczac & Fronczac, 2010}). Fluctuation scaling can be seen in two different forms (Eisler, et al., 2008): The ensemble fluctuation scaling and the TFS (^{Eisler & Kertész, 2006}; ^{Fronczac & Fronczac, 2010)}. Since the discovery of the fluctuation scaling by Smith (1938) and Taylor (1961), this property has been found in complex systems in areas such as social science, life science, ecology, economics, financial markets, climatology, city traffic, complex networks, physics, etc. (Eisler, et al., 2008; ^{Fronczac & Fronczac, 2010}). On the other hand, temporal fluctuation scaling has also been observed in the daily return time series of cases and deaths per day in the world, but this time following the Lévy alpha-stable distribution. The time correlation for these daily return time series has been detected by calculating the Hurst exponents. Regarding this property, it is important to mention that time correlation means that the autocovariance function decays exponentially by a spectral that tends to infinity (^{Stanley, 1971}; ^{Mantegna & Stanley, 2000}; ^{Grau-Carles, 2000}) in such a way that, at the critical point, the said exponential decay turns into a power-law decay (^{Varotsos, et al., 2002}). It should be emphasized that time series associated with different types of complex systems exhibit time correlations as evidenced in nucleotide sequences, neuron activity, cardiac interbeat intervals, earthquake processes, variation of daily maximum temperatures, literary texts, stock index variations, electric signals, etc. (^{Martín-Montoya, et al., 2015}). It is important to insist that time correlation is detected by calculating the Hurst exponent (Martín-Montoya, et al., 2015), a property also detectable if the time series is not large enough (^{Sánchez-Granero, et al., 2008}; ^{Trinidad-Segovia, et al., 2012}). Power-law correlations in time can be measured by a scaling parameter or fractal dimension, but if the time series is self-similar, this fractal dimension is related to the Hurst exponent (Feder, 1994; ^{Chen, et al., 2008}).

The present study showed that the daily cumulative number of confirmed cases and deaths in world countries from March 31, 2020 to September 15, 2021 followed lognormal distributions with two parameters using a similar methodology to that developed to fit the gamma probability density function to wealth distributions obtained in the context of a non-conservative kinetic model of wealth exchange with production saving (^{Quevedo & Quimbay, 2020}). As for the daily cumulative numbers of confirmed cases and deaths among world countries for a given period, it was evident that the standard deviations and the means were related as power laws, which implies the existence of TFS in the system using an analogous procedure to that presented by ^{Abril & Quimbay (2021)} to show that the TFS and the temporal evolution of its exponent in some financial time series can be described by using the path-integral formalism.

In the second part of the paper, the same procedure was used to show that the TFS is also observed in the daily return time series of cases and deaths per day, but this time with a Lévy stable distribution. For these daily return time series, which are not large enough, the existence of time correlation was detected by calculating the Hurts exponents using a modification of the classical R/S method as in ^{Martín-Montoya, et al. (2015)} in their calculation of the Hurst exponents for seismic time series in Colombia given that they were not large enough.

It is important to emphasize that at this stage, the present study does not aspire to explain why the complex system properties identified emerged, as that is the objective of later phases of the research program. In the future, then, the hope is to formulate stochastic models aimed at explaining from a microscopic perspective the emergence of those properties. It is also relevant to mention that the second stage of the aforementioned research program will aim at showing how the properties identified in the spread of COVID-19 in the world also express themselves at the country level, or that such properties exist inside countries when the spread is considered at this level.

Data related to the COVID-19 pandemic spreading in the world

The starting point of this work is the data from confirmed cases and deaths by COVID-19 available in ^{GitHub (2021)}. Thus, the cumulative number of confirmed cases per country has been organized daily starting on January 1, 2020, with 27 confirmed cases in China and ending on September 15, 2021, with 226.477.919 confirmed cases in 193 countries and 721 confirmed cases in an international zone. The cumulative number of deaths per country were also organized daily, starting on January 11, 2020, with one death in China and ending on September 15, 2021, with 4.666.083 deaths in 186 countries and 15 deaths in the international zone. On the other hand, the number of new confirmed cases per day was also daily organized per country starting on January 2, 2020, with no confirmed new cases in China and ending on September 15, 2021, with 563.272 confirmed cases in 160 countries. Finally, the number of new deaths was daily organized per country starting on January 12, 2020, with no confirmed new deaths in China and ending on March 15, 2021, with 10.358 new deaths in 132 countries.

Figure 1 shows that the increase in the number of countries with confirmed cases of COVID-19 (curve in blue color) was large for the period between January 1, 2020, and March 31, 2020, increasing from 1 to 180 (179 countries and 1 international zone) while it was much smaller between April 1, 2020, and September 15, 2021, from 180 to 194. Likewise, the number of countries with deaths by COVID-19 (curve in orange color) passed from one on January 11, 2020, to 127 on March 30, 2020 (126 countries and 1 international zone) while it increased from 130 on April 1, 2020, to 187 on September 15, 2021. Figure 1 also shows a similar behavior for both curves except for a threshold of an approximate average of 18 days: the blue curve grew from 32 on February 22, 2020, to 148 on March 17, 2020, while the orange curve grew from 27 on March 10, 2020, to 147 on April 7, 2020.

Figure 1 Evolution of the number of countries with confirmed cases and deaths by COVID-19 

As shown in figure 1, the spread of the COVID-19 pandemic had a vertiginous progress in the world as the number of countries with confirmed cases and deaths grew strongly from January 1, 2020 to March 30, 2020. During this period, the daily evolution of the number of countries had an analogous behavior to the time evolution of the entropy in the case of a close economical system when states are out of equilibrium (^{Drăgulescu & Yakovenko, 2000}). The economic system described with the kinetic model of money exchange by Drăgulescu & Yakovenko (1990) reaches the Boltzmann-Gibbs distribution of money only when the entropy is at a maximum and the system is in equilibrium. It has been proven that from January 1, 2020 to March 30, 2020 the daily distributions of cumulative numbers of confirmed cases and deaths among world countries do not fit by using any of the probability density functions. Moreover, in the next section it will be evident that only after March 30, 2020 these daily distributions were very well fitted to lognormals.

Temporal fluctuation scaling

The daily cumulative numbers of confirmed cases and deaths in the world distributed among countries fitted by lognormal distributions with two parameters of the form

where o is the shape parameter and |r the location parameter. The goodness of fit was checked in all the cases by the following significance tests: Kolmogorov-Smirnov, Anderson-Darling, and Chi-square. In all cases, the goodness of fit was very good from March 31, 2020 on. According to the significance tests, it should be noted that these fits to lognormal distributions were not completely valid for the period before March 31, 2020 possibly because before this day the system was in a non-equilibrium state as suggested by the results in figure 1.

Table 1 presents the results of the fit parameters to lognormal distributions of two parameters for daily cumulative numbers of confirmed cases and deaths in the world by taking their distributions among the countries and shows the values of the fit parameters to lognormal distributions o and |r for 36 days. As already mentioned, for all cases, the goodness of the fit was very good. The daily values of the mean (M) and the standard deviation (SD) are also shown in the Table.

This system presents TFS because SD and M are related as a power law of the form

where C is a factor and D is the TFS exponent. Using the SD and M values in table 1, for the 15 data from the period between March 31, 2020 and October 31, 2020, the fit parameters to the power law for cumulative numbers of confirmed cases were C=5.2073 and D=0.9839, with a coefficient of determination R²=0.9981, while for the cumulative number of deaths they were C=7.0154 and D=0.9281 with R²=0.9984. Additionally, for the 9 data from the period between November 15, 2020 and March 15, 2021, the fit parameter to power law for the cumulative numbers of confirmed cases were C=2.4953 and D=1.0369 with ii!²=0.9973, while these parameters for the cumulative numbers of deaths were C=3.0786 and D=1.0107 with R²=0.9984. Finally, for the 12 data from the period between March 31, 2021 to September 15, 2021, the fit parameter to power law for the cumulative numbers of confirmed cases were C=21.556 and D=0.8732 with R²=0.9905, while for the cumulative numbers of deaths they were C=11.667 and D=0.8733 with R²=0.9951. It can be seen that for the first period (March 31, 2020-October 31, 2020), the TFS exponents were less than one for both cases, for the second period (November 15, 2020-March 15, 2021) they were greater than one, and for the third period (March 31, 2021-September 15, 2021) they were less than one.

Astonishingly, this system presents another temporal property, but in this case, the location parameter |i of the fit to a lognormal distribution and the mean (M) are related also as a power law of the form

where E is a factor and F is the exponent. For the first period (March 31, 2020- October 31, 2020), the fit parameters to the power law for the cumulative numbers of confirmed cases were E=1.4442 and F=0.1561 with R²=0,9782 while for the cumulative numbers of deaths they were E=0.4382 and F=0.2979 with R²=0,9982. For the second period (November 15, 2020-March 15, 2021), the fit parameters to the power law for the cumulative numbers of confirmed cases were E=2.7213 and F=0.1017 with R²=0.9917, while for the cumulative numbers of deaths they were E=1.7543 and F=0.1419 with R²=0.9755. For the third period (March 31, 2021-September 15, 2021), the fit parameter to the power law for the cumulative numbers of confirmed cases were E=2.2306 and F=0.1167 with R²=0.9823, while for the cumulative numbers of deaths they were E=0.8404 and F=0.2172 with R²=0.9713.

For periods of half a month (T), the standard deviations (SD) and the means (M) calculated daily both for the cumulative numbers of confirmed cases and deaths in the world (cumulative numbers that are distributed among the countries) were related as power laws of the form

where C(T) is the factor depending on the period T and D(T) is the TFS exponent depending on T. This power law relation means that for periods of half a month this system also presented TFS. The fit parameters to the power law are presented in table 2; it is evident that the R ² in all the cases is remarkably close to 1 indicating an excellent fit to a power law in all the cases. For the case of cumulative numbers of confirmed cases, the smallest TFS exponent (D-cases) was 0.5344 for the period July-I-21 while the largest was 1.1704 for the period May-I-21. On the other hand, in the case of cumulative numbers of deaths, the smallest TFS exponent (D-deaths) was 0.6192 for the period July-II-21 while the larger was 1.1943 for the period January-I-21.

Figure 2 shows the time evolution of the TFS exponents for cumulative numbers of confirmed cases (D-cases) and those of deaths (D-deaths) from the values in table 2: TFS exponents for cumulative numbers of confirmed cases (curve in blue color) had five peaks in the periods April-I-20, July-I-20, January-I-21, May-I-21, and September-I-21 and those for cumulative numbers of deaths (curve in orange color) also had five peaks but in the periods April-II-20, August-II-20, January-I-21, June-I-21, and September-I. It should be noted that these peaks were associated with the COVID-19 pandemic peaks registered until now considering the time series of confirmed cases per day and deaths per day in the world (GitHub, 2021). Based on these time series, the following procedure was applied to each of the time series: (i) The average of the data for 20 consecutive days served to construct the average time series; (ii) then, these were normalized with respect to their maximum value (normalized average time series), and finally, (iii) attention was focused on these normalized average times series for cases per day and deaths per day. For the cases per day, it was evident that the peaks were on April 13, 2020, July 30, 2020, January 8, 2021, May 6, 2021, and September 1, 2021, while for deaths per day, the peaks were on April 24, 2020, August 16, 2020, January 29, 2021, May 9, 2021, and September 3, 2021.

Figure 2 Time evolution of TFS exponents for the cumulative numbers of confirmed cases and deaths in periods of a-half month 

The correlation coefficient between the TFS exponents represented in the two curves in figure 2 is 0.421. If it is assumed that the orange curve has a one-period lag with respect to the blue curve, then the correlation coefficient is 0.727, which suggests that the behavior of confirmed cases per day would be reflected approximately 15 days after in the behavior of the deaths per day and, therefore, the existence of a 14-day lag between the times series of confirmed cases per day and the deaths per day may be confirmed. To do so, the correlation coefficients are calculated as a function of the lag in the number of days for the original time series of confirmed cases per day and deaths per day (curve in blue color in figure 3), and for the normalized average time series of confirmed cases per day and deaths per day (curve in orange color in figure 3). In figure 3 the correlation coefficient shows a maximum value in both cases with a 14-day lag so the maximum value for the time series of confirmed cases per day is 0.871 while that for the normalized average time series of confirmed cases per day is 0.924. It is clear that the correlation coefficients calculated for the time series of confirmed cases per day and deaths per day exhibit an oscillating behavior with a 7-day period.

Figure 3 Correlation coefficients as a function of the lag in the number of days both for the original time series of confirmed cases per day and deaths per day and for the normalized average time series of confirmed cases per day and deaths per day 

Time correlation

Considering the time series of confirmed cases per day and deaths per day in the world (^{GitHub, 2021}), the daily return times series of confirmed cases per day CD(t) and deaths per day DD(t) were constructed by defining the daily return of confirmed cases per day (DRCD) in the following form: DRCD(t+1) = [CD(t+1)-CD(t)]/CD(t); similarly, the daily return of deaths per day (DRDD) was expressed as DRDD(t+1)=[DD(t+1)-DD(t)]/DD(t). These definitions implied that the daily return represented the percentage change of the respective quantity between two successive days: t + 1 and t. In the next step, January 22, 2020, was defined for both daily return time series as the first date and September 15, 2021, as the last date, each with a total of 603 data. Starting from these daily return time series, time subseries having different numbers of data were constructed for each case so that for all the time subseries, the first data was on January 22, 2020. The number of data of the first daily return time subseries was 60, corresponding to the period from January 22, 2020, to March 21, 2020. The second daily return time subseries had 75 data, corresponding to the period from January 22, 2020, to April 5, 2020, which meant that the increase in the number of data between two consecutive time subseries was 15. Thus, the third time subseries had 90 data, corresponding to the period between January 20, 2020 and April 20, 2020. Taking the possible increments of 15 data, the total number of time subseries that could be constructed was 37. For instance, the time subseries number 36 had 585 data corresponding to the period from January 22, 2020 to October 28, 2021, and the last one had 603 data corresponding to the period from January 16, 2020 to September 15, 2021, which matched the original daily return time series.

For each one of the time subseries of confirmed cases per day and deaths per day, the following quantities were calculated: (i) mean; (ii) standard deviation; (iii) the four fit parameters to a Lévy distribution notated as α, β, μ, and σ, and (iv) the Hurst exponent (H) calculated by using a modification of the classical R/S method as in ^{Martín-Montoya, et al. (2015)}. The results of these calculations are presented in table 3.

The existence of another type of TFS is apparent given that SD(t) and M(t) are related as a power law of the form

where K is a factor and L is the TFS exponent. Using the values of SD(t) and M(t) shown in table 3, the fit parameters to the power law for daily returns of confirmed cases per day were K=5.2595 and L=0.5372, with R²=0.9997, while these parameters for daily returns of deaths per day were K=6.0527 and L=0.5159, with R²=0.9997.

Astonishingly, this system had another TFS because the standard deviation of the daily returns of deaths per day SDDRD(t) and the mean of the daily returns of confirmed cases per day DMRCC(t) were related as a power law of the form

where P is a factor and Q is the TFS exponent. Using the values of SDDRD(t) and MDRCC(t) in table 3, the following fit parameters to this power law were obtained: P=10.905 and Q=0.5225, with R²=0.9994. The reason for this result is that SDDRD(t) and the standard deviation of daily returns of confirmed cases per day SDDRCC(t) had a linear relation of the form SDDRED(t) = 2.062* SDDRCC(t) + 0.0581, with R²=0.9999. This implies that the correlation coefficient between the standard deviations SDDRED(t) and SDDRCC(t) in table 3 should be 1 and, effectively, this correlation coefficient was 0.9999718. This result implies that the behavior of the daily returns of deaths per day was highly correlated to the behavior of the daily returns of confirmed cases per day as validated by the fact that the mean of the daily returns of deaths per day MDRD(t) was related with the mean of the daily returns of confirmed cases MDRCC(t) through the linear relation MDRD(t) =3.074* MDRCC(t) - 0.0226, with R² =0.9999.

On the other hand, the Hurts exponents of the daily returns of confirmed cases per day and those of the daily returns of deaths per day in table 3 are represented in figure 4 where the evolution on time for the Hurts exponents of the daily returns of confirmed cases per day (blue color curve) and for those of the daily returns of deaths per day (orange color curve) are shown. The time correlation behavior as anti-persistent (H<0,5) was detected for the daily returns of confirmed cases per day in the following periods: January 22, 2020, to March 21, 2020; January 22, 2020, to April 5, 2020; January 22, 2020, to June 4, 2020; January 22, 2020, to June 19, 2020; January 22, 2020, to July 19, 2020, and January 22, 2020, to August 3, 2020. The behavior of the time correlation was persistent (H>0,5) for the other 31 periods. This persistent behavior grew continuously on time as the Hurst exponent for the period between January 22, 2020, and August 18, 2020, was 0.5204 while between January 22, 2020, and September 15, 2021, it was 0.6665. Additionally, the behavior of the short-long correlation in time was persistent (H>0,5) only for the daily returns of deaths per day for the 37 periods, a behavior that grew continuously with time as the Hurst exponent for the period from January 22, 2020, to March 31, 2020 was 0.5286 and from January 22, 2020 to September 15, 2021 it was 0.8348.

Figure 4 Evolution on time for Hurts exponents of the daily returns of confirmed cases per day (blue color curve) and the daily returns of deaths per day (orange color curve). 

Conclusions

The availability of the time series of confirmed cases and deaths (^{GitHub, 2021}) from January 1, 2020, to September 15, 2021, allowed for the detection of complex system properties such as lognormal laws, temporal fluctuation scaling, and time correlation in the spreading of the COVID-19 pandemic around the world.

During the period from January 1, 2020, to March 31, 2020, the spread of the COVID-19 pandemic had vertiginous progress in the world as the number of countries with confirmed cases and deaths increased. Besides, during the period from January 1, 2020, to March 30, 2020, the daily distributions of cumulative numbers of confirmed cases and deaths among countries did not fit any probability density function because the system behaved as if it were in a non-equilibrium state. Consequently, only after March 30, 2020, these daily distributions perfectly fitted lognormals of two parameters as the system behaved in an equilibrium state.

For periods of half a month, the standard deviations and the means calculated daily both for the cumulative numbers of confirmed cases and for those of deaths in the world were related as power laws as a manifestation of the TFS. Specifically, the values of the TFS exponents depended on the period considered and their values were between 0.53 and 1.19. Moreover, the TFS exponents for cumulative numbers of confirmed cases showed five peaks in the periods April-I-20, July-I-20, January-I-21, May-I-21, and September-I-21 and the TFS exponents for cumulative numbers of deaths also had five peaks but in the periods April-II-20, August-II-20, January-I-21, June-I-21, and September-I and they were associated with the five peaks that the COVID-19 pandemic around the world had had until September 15, 2021.

Starting from the daily return times series of confirmed cases and deaths per day and taking all the possible increments of 15 data between two successive time subseries for each case, the total number of daily returns time subseries constructed was 37. For all these daily return time subseries, the mean, the standard deviation, the fit parameters to a Lévy distribution and the Hurst exponent were calculated. It should be mentioned that the Hurst exponents were calculated by using a modification of the classical R/S method used when the time series is not large enough. The standard deviations and the means of all the daily returns time subseries for confirmed cases and deaths per day were related as power laws indicating the existence of another type of TFS in the system. Also, the standard deviations of the daily returns of deaths and the means of those of confirmed cases per day were related as a power law meaning that the daily returns of deaths per day were highly correlated to those of confirmed cases per day.

The presence of complex system properties in the spreading of the COVID-19 pandemic at world level was determined by considering the world as constituted by countries, although they were also found when smaller regions were considered. In another study in progress, it has been shown that these properties are also present at country level considering countries as constituted by states and, similarly, they are present as well at state level when it is assumed these are constituted by counties indicating that the properties identified here are invariant when the spatial scale changes.

These results are consistent with the stochastic nature of the spread of the COVID-19 pandemic around the world. Consequently, a stochastic model aimed at explaining it should lead to the emergence of the complex system properties identified in this paper.