Universality of delay-time averages for financial time series: analytical results, computer simulations, and analysis of historical stock-market prices

We analyze historical data of stock-market prices for multiple financial indices using the concept of delay-time averaging for the financial time series (FTS). The region of validity of our recent theoretical predictions [Cherstvy A G et al 2017 New J. Phys. 19 063045] for the standard and delayed time-averaged mean-squared ‘displacements’ (TAMSDs) of the historical FTS is extended to all lag times. As the first novel element, we perform extensive computer simulations of the stochastic differential equation describing geometric Brownian motion (GBM) which demonstrate a quantitative agreement with the analytical long-term price-evolution predictions in terms of the delayed TAMSD (for all stock-market indices in crisis-free times). Secondly, we present a robust procedure of determination of the model parameters of GBM via fitting the features of the price-evolution dynamics in the FTS for stocks and cryptocurrencies. The employed concept of single-trajectory-based time averaging can serve as a predictive tool (proxy) for a mathematically based assessment and rationalization of probabilistic trends in the evolution of stock-market prices.


Evolution of stock-market prices and stochastic processes with multiplicative dynamics
Variations of stock-market prices are generally believed to be unpredictable, with completely random price fluctuations obeying memory-less random walks [1][2][3][4]. This was first argued in 1900 by Bachelier [1], who introduced (before Einstein and Langevin) the concept of arithmetic Brownian motion (BM), developed further in 1908 by Bronzin [2]. Bachelier and Bronzin laid the foundation of modern financial mathematics, the starting point for applying stochastic processes and statistical analysis to time-evolution of stock prices, S(t). Generations of brilliant scientists-economists, econometricians, econophysicists, specialists of financialtime-series (FTS) analysis, etc (including many Nobel-prize winners)-have been focusing their long-term efforts on unraveling the general underlying functioning principles of financial markets and predicting the governing laws of formation and evolution of asset prices  (see also the key books [46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63]). One of the main conjectures for financial-market models is that consecutive price changes and respective returns, r n = [S(t n ) − S(t n−1 )]/S(t n−1 ), are uncorrelated in time (Bachelier's first law [1]) and can be represented by independent identically distributed Gaussian random variables.
The ground-breaking discovery in financial mathematics was the development of the paradigmatic geometric BM (GBM) process (also called 'exponential' or 'economic' BM). It was 'revived' in 1965 by Samuelson [12] from Bachelier's works and the respective closed-form option-pricing formula utilizing the GBM process was invented by Black, Scholes, and Merton (BSM) in 1973 [21][22][23][24][25][26]. This famous BSM model was generalized i.a., by Cox, Ingersoll, and Ross [32], Hull and White [33], Wiggins [J B Wiggins, Option values under stochastic volatility: theory and empirical estimates, J. Financial Econ., 19, 351 (1987)], Stein and Stein [3]-even a tax-free trading would be hardly profitable. Real financial markets are extremely complex and outof-equilibrium systems [86], often operating at a limited 'liquidity'. Thus, some short-horizon trends of price evolution may potentially be predictable from the incoming information and based on statistical analyses of price-fluctuation patterns of market response to similar information in the past.
The traders' anticipation, numerous sources and types of noise [120], delays and memories are existential for trading and 'functioning' of the market. In traditional economic theories, the trading events act as an instrument to dynamically probe and approach the market equilibrium (Walras' tâtonnement process, see reference [86]), via regulating supply and demand [86]. This equilibrium can, however, be illusive or even artificially created, e.g., by 'big players' controlling the market [115]. These (usually much more informed investors) trade profitably knowing the 'true price' [35,115,118,120,127], often at the expense of less informed traders [128]. 6 The cost of information acquisition prevents the market from being 'informationally efficient' [118,129,130] impeding (or even destroying) the 'fair play' among the participants and often violating the market-efficiency hypothesis.
As an indication of private-information-driven changes of volatility, the fluctuations of FTS were demonstrated (over multiple time domains) to be far too large [131], and with price changes being only partially attributable to important financial news [127,132]. The variations of stock-market prices were also shown to fail in reflecting rationally changes in fundamental values [133]. This questions the relative role of exogenous (external) versus endogenous fluctuations (or 'self-generated' price changes) in such price-formation processes. The latter are often irrational and affected possibly more by behavioral-finance principles and human decision-making principles, rather than by real values.
Long-term speculative exponential GBM-like growth of stock-market prices is, thus, supplemented by disproportionably volatile and often irrational short-term price variations (with noise-versus information-based trends [86,120,134] being impossible to separate rationally). Moreover, market 'overreaction' to new financial information is well known [121,135]. This can be a useful instrument for the traders to artificially create bigger winning margins first in order later to 'capitalize' larger profits on a price 'bounce-back' 7 .
The exponential GBM-like price growth can, in addition, be interrupted by hardly predictable [85]-but possibly inevitable or even pre-programmed-price drops at times of market crashes and financial bubbles [80,81,85,[139][140][141][142]. The 'efficient market hypothesis' [12,19] implies that the variation of the unanticipated part of stock prices should be a martingale [25], with no correlations in price differences [4,12] (see also reference [136]). The speculative bubbles are based on overoptimistic expectations, the herding behavior of market participants [143][144][145][146]-particularly in periods of extreme volatilities (in high-risk-high-reward times), and pure greed-all accelerating the price dynamics and growth, with a possibly superexponential [147,148] price 'explosion' near the bubble. The burst of speculative bubbles and subsequent market crushes are inconsistent with the hypothesis of efficient markets [5,19,27] (the empirical evidence of the latter are, in fact, insufficient [133,136,138,149]).

Outline of the paper
The paper is organized as follows. In section 2 we recapitulate on the recent results of the time-averaged MSD (TAMSD) δ 2 i (Δ) analysis for GBM and historical stock-market prices [150]. In section 3 we present the essential details of the analytical model of GBM. In section 3.1 we introduce the formalism for GBM, in section 3.2 the key expressions for its delayed TAMSD δ 2 d,i (Δ) are derived. In section 4 we discuss the numerical solution of the stock-market price-evolution equation and present the details of computer simulations for GBM. The findings for δ 2 i (Δ) and δ 2 d,i (Δ) are described and compared to analytical results for GBM. In section 5 we present the main results of the FTS-analysis in terms of the delayed TAMSDs for a number of indices. In section 6 we summarize the results, list general conclusions, and discuss possible further developments.
Our main target here is to present the full analytical results for δ 2 d,i (Δ) (generalizing and extending those of reference [150]), to confirm the agreement of its behavior with the results of computer simulations and the analysis of evolution of stock-price data (at all lag times). We present some auxiliary figures in appendix A, some derivations in appendix B, and the data-driven procedure of parameters estimation in appendix C.

Previous results for the TAMSD of GBM: aged and delayed properties
In contrast to standard moving/'rolling' average-averaging prices over n preceding points along an FTS-the TAMSD operates with all points of an FTS S i (t) of length T yielding the trajectory-averaged displacements at varying values of the lag time, 0 < Δ < T, via integrating the squared stock-price increments as [150] Here, S i (t) is the price of the ith stock. In equation (1) the price increments shifted in time by Δ are squared and averaged over the FTS. The discrete interpretation of (1) for equidistantly sampled input data-such as in historical FTS-is straightforward. The TAMSD concept-ubiquitously used for the analysis of singleparticle-tracking data [154][155][156]-was employed to the FTS-analysis in reference [150], a foundation for the current study. The mean TAMSD for N independent trajectories constructing a 'satisfactory' statistical ensemble is We summarize below the basic concepts and main results of reference [150], focusing in particular on the concepts of the aged and delayed TAMSDs. First, we confirmed that the TAMSD behaved strictly linear with the lag time Δ for a large number of companies of various classes, with This behavior is consistent [150] with the predictions of the standard-GBM model highlighting, thus, weak ergodicity breaking [155] emerges for this process at short lag times, at Δ T. Specifically, the exponentially growing MSD for GBM (16) contrasts the linearly growing mean TAMSD in equation (3). We extend this analysis here for all lag times Δ and present nontrivial features of the TAMSD at later stages of the trajectories.
Second, we introduced the aged and the delayed TAMSDs and enumerated them as functions of the aging t a and delay t d time for a large number of FTS [150]. The modified TAMSDs (4) and (5) were also computed analytically for standard GBM with a constant volatility σ. The mean aged TAMSD was found (in the limit of short lag times, Δ T, and with no drift) to grow nearly exponentially with t a both for real FTS and for GBM predictions, namely The stock-specific or 'idiosyncratic' factor σ 2 i in equation (6) was, however, shown [150] to yield a spread of the distribution of log δ 2 a,i (Δ) δ 2 i (Δ) when this ratio is examined as a function of aging time t a for different stock-indices at a fixed lag time Δ. The most relevant difference between the aged and delayed TAMSDs is the fact that δ 2 d,i (Δ) assigns more 'weight' to the data points toward the end of FTS. This key difference yields (for often nearly exponentially growing FTS) some universal and parameter-free relations for the delayed TAMSD.
Third, the most interesting finding of reference [150] was the universal behavior predicted for the delayed TAMSDs versus t d for all stock indices examined, with in the limit of short delay times and short lag times, at The parameter-free master curve (7) obtained for many FTS is consistent [150] with the GBM-based solution in the same limit (8). The concept of the delayed TAMSD (5) introduced in reference [150] is, therefore, especially useful for the analysis of fast-growing FTSs and for assessing their universal parameter-free characteristic features, such as (7).
Our current analysis extends the consideration of δ 2 d,i (Δ) to the entire range of lag-and delay-times. We demonstrate, i.a., that the behavior of the later parts of the log δ 2 d,i (Δ) δ 2 i (Δ) versus t d curves with a fasterthan-linear growth are not inaccuracies of the basal law (7) observed at short lag times. This 'nonlinear' growth is rather a real feature of FTS, which agrees with both the GBM predictions and results of computer simulations, as shown below. We also clarify the implications of financial crashes on deviations from typical GBM-based evolution of FTS observed for the TAMSD in crisis-free times.

GBM: equation, general solution, model parameters, and moments
Since Bachelier [1], the evolution of a priori uncorrelated [3] stock-market prices is described using the concept of random walks. According to the BSM model, the asset price S(t) obeys a stochastic differential equation driven by multiplicative [152] noise, where W(t) is the standard Wiener process defined via zero-mean white Gaussian noise ξ(t) as Equation (9) is considered in the Itô formalism. The constant parameters μ and σ denote the drift and volatility of this process, respectively (see the discussion in section 6.2). After applying Itô's lemma [57,153], one can derive from (9) the partial differential equation yielding the Black-Scholes formula (the BSM model). The solution of equation (9) in Itô representation defines GBM as a Markovian process (exponential of BM X(t)) of the form where S 0 ≡ S(0) is the initial price, and The first and second moments of GBM, defined as are obtained via averaging (11) with the log-normal distribution 8 [22,23] The latter satisfies the Fokker-Planck-like partial-differential equation (see references [25,28,106,108]) with the δ-function-like initial condition, P(S(0), t = 0) = δ(S(0) − S 0 ). The averaging procedure (13) yields For brevity, we use hereafter the following definition The variance of GBM is given by 8 The log-normal distribution of stock prices after any time interval is another fundamental postulate of the BSM model. The other assumptions are, i.a., (on the 'model' level) that trading is continuous in time, the price variations are continuous and jump-free in time, the adjustment to a new information is instant and memory-less, and the actual price reflects all information available. Additionally (on the 'executional' level), no restrictions on short-term sells are imposed, no taxes or retail commissions are to be paid from possible profits, no general transaction/trading fees exist, etc (see reference [25] for option-pricing models with discontinuous price variations). Importantly, for GBM models the volatility is known parameter constant in time, that is a rather unrealistic assumption [120].
For GBM, the exponential growth of S 2 (t) in (16) stems from the BM of the log(price), with the variance of fluctuations growing linearly with time (as in the original Bachelier's study [1]).

TAMSD of GBM: nonaged, aged, and delayed cases
The relations (16) and (18) follow from the general 'stochastic' solution (11) and one-and twopoint probability-density functions of the Wiener process, (19) and (at t 1 > t 2 ) computing, respectively, and Similarly, the integrand of TAMSD (1) can be expressed via that yields Taking the final TAMSD integral over t from the last exponent in (24), one gets for the standard (nonaged), aged, and delayed TAMSDs of GBM, respectively, and From equation (25) in the absence of drift, at short lag times (Δ T) and for long traces (σ 2 T 1) one gets the fundamental relation (3) [150]. From equation (26) relation (6) follows at short lag and aging times, {t a , Δ} T. Finally, from (27) at short lag and delay times ({t d , Δ} T) one arrives at the fundamental law (7). The TAMSD analytical expressions (25)-(27) is our first key result. Although derived to get the short-time asymptotes before [150], the complete expressions are presented here for the first time. We refer to appendix B for an alternative TAMSD derivation.

Numerical integration scheme
As follows from (12), to simulate GBM the Wiener process W(t) needs to be exponentiated. The increments of W(t) are independent and normally distributed random variables with mean zero and variance (t 1 − t 2 ), so that As W(t 1 ) ∼ N (0, t 1 ) for t 2 = 0, using (12) for the mean and variance of X(t) we get, respectively, X(t) = μ − σ 2 /2 t and The subsequent values of the Wiener process are generated using where n = {0, 1, 2, . . . ,N − 1} is the discretization index and Z ∼ N (0, 1) obeys the normal distribution for the time points For a stochastic process X(t) defined by (12) using expression (30) we get On the equidistant time-grid with the step-size δt = t n+1 − t n the simple recursion formula [33] for GBM then becomes From equation (11) follows that ratio S(t)/S 0 is distributed log-normally, namely see equation (14). We checked that computer simulations (see figure 1) indeed produce the correct distribution (34), figure AA1. The time-step in the simulations δt was chosen to be 1 business day or 1/252 of the fiscal year (unless explicitly specified otherwise). This makes the results of our simulations straightforwardly comparable with the findings of the analysis of real FTS for classical stocks examined with the same time-step in section 5.1.

TAMSD
Extensive GBM-based computer simulations deliver novel features for the current study, as compared to the original data-focused study [150]. Some individual in-silico-generated GBM trajectories and the respective For longer lag times, a faster growth of the TAMSD and faster-than-linear scaling with the lag time are often observed. In this later region-due to the worsening statistics inherent for the TAMSD definition (1) [154,155]-considerably fewer increment values are available for time averaging and, therefore, the observed variations of δ 2 i (Δ) with Δ are much stronger. We mention that for a steadily increasing S i (t) realization of GBM in simulations the TAMSD magnitude increases for longer trace lengths T, while for stalling or dropping price realizations S i (t) the TAMSD loses this systematic trend, compare, e.g., figures AA2(B) and AA2(H).
In figure 1 we show N = 25 different TAMSD realizations for GBM-each set generated at the same magnitudes of drift and volatility-as functions of the lag time and for the trace length of T = 35 years. We recognize the linear scaling of the TAMSD in the region Δ T, equation (35). We also find that for relatively small σ-when the dynamics is mainly impacted by a nonzero value of the drift parameter μ (responsible for the exponential growth of the first moment of GBM, equation (16))-the TAMSDs reveal a faster-than-linear growth at later lag times. The spread of individual TAMSDs in a given set of generated curves decreases for progressively smaller σ values, as expected for the situation when the influence of 'randomness' or the underlying volatility of the process decreases, see figures 1(A) and (B).
For larger σ values the reproducibility of TAMSD realizations for GBM naturally decreases: each TAMSD becomes more volatile in magnitude and as a result their spread around the mean value increases dramatically, see figure 1(G). The mean TAMSDs are strongly shifted toward the top region of the distribution of individual TAMSDs, as illustrated in figure 1(G): we find that often a single trajectory with an extremely large magnitude dominates the mean TAMSD. With increasing σ the distributions of TAMSDs in figure 1 reveal larger spreads, typical for the GBM process with higher volatilities. T (shown for some GBM traces also in figures AA2 and 1) is clearly present after averaging. Note that a large set of N = 10 6 independent GBM trajectories was used for averaging here, and in most plots with the results of simulations. This very large statistical ensemble is necessary to compensate for outliers and large-deviation trajectories for this multiplicative and innately highly varying GBM process (such a trace is capable of biasing the mean TAMSD for the entire ensemble). The size of a proper averaging ensemble was argued to increase exponentially [197] with the number of points in the time series,N, for such processes.
We observe that with increasing volatility σ the magnitude of the TAMSD dramatically increases, see figure 2(A) for the δ 2 (Δ) variation for GBM. The behavior of the normalized TAMSD, namely of

Delayed TAMSD
The behavior of individual realizations δ 2 d (Δ)-computed for a single GBM trace for varying delay times t d -versus the lag time Δ is shown in the right column of figure 3. A key property of the delayed TAMSD [150] is that the later parts of the FTS contribute progressively stronger at longer delay times. We find that, similar to the trends of the standard TAMSD shown in figure 1(G), the variation of the magnitudes of δ 2 d (Δ) increases for larger volatilities, see figure 3(H).
The variation of the delayed TAMSD computed at the shortest lag time, δ 2 d (Δ 1 ), with delay time t d is shown in figure 4. We focus on the shortest lag time because the magnitudes of the standard and delayed TAMSDs are statistically most reliable at Δ/T 1 [150,155]. Here we find, as expected, also a wider spread of individual becomes negative-the fluctuations of the Wiener process cannot ensure any sustainable (exponential) price growth for the resulting stochastic process. For this situation, a finite fraction of GBM traces predict drastically falling prices S(t) at longer times, a feature of a collapse or bankruptcy, see figure 4(D). In economics, μ σ,eff is often called the expected growth rate [112]. The mean value δ 2 d (Δ) is, however, again often dominated by a single large-magnitude trajectory and, thus, only weakly affected by such price drops observed for some strongly volatile GBM trajectories.
At even larger values of σ used in simulations we find extremely volatile behavior of the resulting S(t) and of the respective δ 2 d (Δ) realizations. This fact, in turn, gives rise to large variations of δ 2 d (Δ 1 ) versus the delay time and, ultimately, to irregularities in the behavior of the mean δ 2 d (Δ 1 ) versus t d (not shown). In general, we find that δ 2 d (Δ 1 ) -dependence versus t d evaluated for different σ values mainly differs in magnitude, while keeping the overall functional form of variation with the delay time, figure AA3. Specifically, δ 2 d (Δ 1 ) reveals almost no variation for short-to-intermediate delay times and exhibits a significant increase of the magnitude at t d → T, see figure AA3. The agreement of analytical predictions of equation (27) with the results of computer simulations for δ 2 d (Δ 1 ) is excellent for small-to-moderate volatility values. For an ensemble of N = 10 6 GBM trajectories used for averaging, only the data for the largest μ and σ values are considerably lower than the analytical predictions. We 'scale' the results of simulations and shift them up in magnitude for such situations in order to reach the theoretical asymptote. Figure AA3 shows both the original and 'scaled' results of simulations for δ 2 d (Δ 1 ) . One can argue that for strongly volatile realizations of GBM at large σ values the magnitude of the mean delayed TAMSD can be influenced by a single extreme-magnitude trajectory (see figure 4(D) for an example) that dominates the mean δ 2 d (Δ 1 ) . Even larger ensembles of Figure 5. The same log-ratio log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) as in figure AA4 for the same μ and σ values, but shown for varying trajectory lengths T (see the legend). The curves are the analytical results (27). Parameters: μ = 0.1, σ = 0.2, N = 10 6 . The error bars are smaller than the symbol size for all data points.
GBM traces could be required to encounter such large-magnitude trajectories in simulations and to mitigate the observed disparity of theory-versus-simulations in the magnitude of δ 2 d (Δ 1 ) , as detected at larger σ in figure AA3.
The behavior of the log of the mean delayed TAMSD normalized to the standard TAMSD computed at Δ = Δ 1 , namely log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) , the universal quantity introduced in reference [150], is shown in figure AA4. The simulations reveal the linear regime of equation (7) for this log-ratio at t d T, as expected, and also confirm a faster-than-linear growth of log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) at intermediate-to-long delay times. Excellent agreement with the theoretical result (27) is observed in our computer simulations for the log-ratio of the normalized delayed TAMSD in the entire range of delay times, see figure 5. This is our third key result.
In figure 5 we present the detailed results for log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) for GBM evaluated as function of the delay time t d for varying trace lengths T. We find that for longer traces the magnitude of log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) drops and the region with a faster-than-linear scaling becomes more pronounced, extending into a larger domain of delay times, at t d T. These two trends for the behavior of the delayed TAMSD are fully supported by our analytical results, equation (27). We also mention that for longer trajectories the variation of μ and σ has a considerably smaller effect on the scaling behavior of the log-ratio log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) computed at short delay time t d , see figure AA5.

Ergodicity breaking
Note that the 'time-rearrangement trick' (B9) can also be useful [151] for evaluating higher moments of GBM, such as the fourth time-averaged moment that contributes to the ergodicity breaking parameter [154,155], Computer-generated GBM trajectories reveal a similar scaling for δ 2 (Δ) 2 and δ 2 (Δ) 2 with Δ and, thus, a roughly lag-time-independent value of EB. For trajectory length of T = 35 years we find EB GBM (Δ) ≈ 5T/2, as shown in figure AA6. A detailed analytical study of EB for arbitrary trace length and GBM-model parameters, the simulations-based EB analysis for GBM as well as its applicability to real FTS is the subject of a separate investigation [151]. From simulations we can conclude that GBM is a nonergodic process featuring a finite-rather than vanishing-value of the EB parameter in the limit of long trajectories and short lag times, see figure AA6 at Δ/T 1. Note that the nonergodicity of GBM studied in reference [113] differs in its definition from equation (37). The latter was extensively applied in recent years for assessing weak ergodicity breaking for the in-silico-generated trajectories of various stochastic anomalous-diffusion processes [157][158][159][160][161][162][163][164].

Analysis of historical FTS
With these theoretical concepts and results (in particular, regarding the δ 2 d,i (Δ)), we proceed now to the analysis of historical FTS. For real data we compute the same observables as for the simulated GBM trajectories and examine to what extent the GBM-model predictions are applicable.

Companies/indices
For the chosen FTS of stock-market prices, the log-ratios of the normalized delayed TAMSDs, log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) , demonstrate the universal behavior via collapsing onto a single master curve as functions of t d /T, see figure 6. The most pronounced feature is the fact that at short-to-intermediate delay times-in crisis-free times, when the GBM-model itself is applicable-the FTS yield the delayed TAMSD that follows the linear GBM-conform law (7) [150]. The deviation of log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) as function of t d from this law at later delay times FTS is also in excellent agreement with the GBM-based prediction (27) 10 . The latter is valid at all values of lag and delay times yielding our fourth key result. We stress here that all data points were used in the analysis, without any pre-selection, both for the crisis-containing and crisis-free times.
A severe drop of the delayed TAMSDs at the end of some trajectories in figure 6 is the result of the three latest economic crises; this region with rapidly falling prices is (naturally) in disagreement with GBM predictions. We quantify these deviations below via analyzing the data in log-linear scale, for only later segments of the FTS and for smaller delay-time increments, δt d . The pronounced deviations from the GBM predictions at later delay times are, naturally, due to the respective drastic crisis-induced price drops. In figure 7 we relate the delay  times t d to the actual calendar date of the respective FTS. This way, we relate the time of drops in δ 2 d,i (Δ) to the respective economic perturbations/crisis, compared to plots versus the t d itself, as in figure 6.
In figure 7 we demonstrate a quantitative agreement with the theoretical GBM-based predictions (27) for the actual historical FTS in crisis-free times. We also find that the drops of log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) versus t d indeed corresponds to the periods 1997-1999 and 2008-2009 (the Russian/Asian and financial crises, respectively). The ongoing 2020 world-wide economic decline triggered by the ongoing Covid-19 pandemic is also visible for the later parts of the analyzed FTS, as indicated in figure 8. The latest stock-price variations are particularly well visible and pronounced when the delayed TAMSD data are presented for the last prior-to-crisis months only and with the minimal step of δt d = 1 day in order to resolve daily fluctuations (results not shown). Naturally, rapid variations of the log-ratio of the delayed TAMSDs at later parts of FTS are then resolved much better, as compared to the case of δt d = 1 year shown in figure 7.
After the 2008-2009 crisis with a partial 'resetting' [165,166] of many stock-market indices to a certain extent we observe that some indices started to continue their (exponential) GBM-like price growth again, thereby reloading the price-explosion spring (results not shown). Note also that choosing the end point of FTS well after the peak of a financial crisis smears out the crisis-induced drop of prices when evaluating the TAMSD and reduces the deviations from the GBM log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) ∼ t d /T asymptote, see equation (7). This happens as typically higher prices occur at later stages in the FTS. This is why no decrease is visible, for instance, in the S & P 500 long-time data when plotted in terms of log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) versus t d in October 1987, at the time of S & P 500 crash, see figure AA7 (and also the results of reference [150]).

Cryptocurrencies
Some 'digital' assets, such as BitCoin and other cryptocurrencies, see figure AA9, are extremely speculative, volatile, and nontransparent [146] in their price-formation strategies. Their FTS can potentially feature superexponential [147] price evolution in crisis-free times and at 'bull-market'-conditions. This is a snow-balllike, heating-up, herd-driven [83] phase of price 'explosion' characteristic of financial pyramids (with no 'real value' supporting the growth). Our current preliminary analysis of some cryptocurrencies quite surprisingly, however, indicates a good agreement with GBM-like variation for δ 2 d,i (Δ), see figure 9, but with significantly elevated values of drift and volatility, see appendix C for the detailed description of the data-driven algorithms of determination/extraction of the GBM model parameters.
The dependencies of the trajectory-specific quantifier log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) for the FTS of the 'classical' stocks and of three cryptocurrencies versus the GBM predictions for log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) evaluated for the {drift, volatility} pairs determined from the data are shown in figure 9, both as functions of the delay time t d . The highly speculative nature of the 'contemporary' cryptocurrencies gets reflected in much larger respective drift and volatility values, as compared to the 'classical' stocks, see figures BB2 and BB3. The direct comparison in terms of the variation of the delayed TAMSD versus the delay time-demonstrates an expected much steeper growth of the log-ratio for the cryptocurrencies in the same range of delay times, as visible in figure 9. Despite this, the general trends of the growth of δ 2 d,i (Δ) with the delay time for the cryptocurrencies are well-described by the standard-GBM theoretical model.
We stress that such a comparison is performed in the same range of lag times implies similar time-scales for all the stocks and currencies examined here. The 'agility' and also the anticipated 'life-time' of a given stock or asset (or a class of those) are, however, inherently specific. Evidently, the market of cryptocurrencies is much more dynamic and speculative than that of 'classical' assets/commodities (some of them with a century-long history, such as that of trading and pricing gold). This determines, or at least affects, the characteristic 'internal' time-scale and also sets an 'effective temperature' that 'heats up' the growth and overall dynamics of a given stock or asset, an issue deserving future investigations.

Summary of the main results
We presented the time-averaging-based analysis of the historical FTS of stock-price evolution. We first compared the analytical GBM-based predictions and reported the results of extensive GBM-based computer simulations. Our main focus was on the behavior of the novel observable, the delayed TAMSD [150]. This singletrajectory quantifier-the term 'inherited' from the data analysis in single-particle-tracking experiments-is demonstrated here to reveal a universal behavior for the log-ratio log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) at the shortest lag time Δ = Δ 1 as a function of the delay time t d . This universality in the entire range of delay times (in crisis-free times) complements the initial analysis at t d T in reference [150].
Additionally, our FTS-analysis revealed a transition from a linear to a faster-than-linear growth regime for log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) as function of the delay time t d , in full agreement with GBM-based predictions (in crisis-free times). Certain deviations from this 'ideal' GBM-like behavior found for the delayed TAMSD in real FTS were shown to be clearly attributable to market collapse associated with the dramatic price drop during respective financial crises.
One more novel element of the current analysis compared to reference [150] is the procedure of determination and optimization of the GBM-model parameters and direct fit of δ 2 d,i (Δ) versus delay time, as detailed in sections C2 and C3 of appendix C. We compared several methods to extract σ and μ and employed finally the most rational and consistent approach, see the results of figure 9. We found that the historical FTS for the highly speculative BitCoin demonstrated the extreme growth in magnitude of the log δ 2 d (Δ 1 ) δ 2 (Δ 1 ) within a few years of delay time t d . The similar-in-magnitude growth is achieved by the classical stocks only within ≈30 . . . 40 years. This once again indicates a bubble-like nature and lack of sustainability of BitCoin and other cryptocurrencies. This parameter-determination procedure enabled us to fit and quantitatively compare the TAMSDs for individual historical FTS, that is the fifth key result of our study.

Extensions and future developments of GBM-based models
A more detailed analysis, that would potentially include some generalized GBM-based models featuring superexponential MSD growth, with β > 1, a 'scaled-GBM' process recently introduced in reference [167]. The nonlinear nature of another type was included in the models with the return rate and volatility growing nonlinearly with the price, such that ∼ σS(t) m dW(t) (39) enters the right-hand side of equation (9) (see equation (24) in reference [139]). This simple nonlinear generalization of GBM accounts for a positive feedback and herding behavior (a price-driven model of speculative bubbles). Certain bursts of volatility accompany the approach to and the transition across the 'bubble point' in this kind of models. Anomalous behavior of the markets, particularly with a superexponential price growth [147,148] at and near the bubble corresponds to m > 1 choice in equation (39).
For the standard GBM model, the parameters μ and σ are assumed to have constant values, while in reality they can be complicated functions of time, price, and numerous other factors. Often, the real markets reveal, e.g., an increased volatility when the prices drop [168,169], nonlinear price-volatility correlations [93,170], anticorrelations of returns and future volatility values [171], the periods of persistent volatility [172,173], correlations of the trading volumes and intraday volatility [175], and rich dynamics of (nonstationary) intraday price increments [174,176]. The effects of 'volatility clustering'-with large price changes being followed by large ones (and vice versa)-and possible long-term memory in volatility are the general 'stylized facts' [177] for some FTS. We refer to references [32,33,36,178] for the models of stochastic volatility and to reference [179] for the analysis of finely and coarsely defined volatilities. This feature is treated, i.a., in the models of generalized autoregressive conditional heteroscedasticity (GARCH) [44,119,149,[180][181][182][183][184] and autoregressive moving average (ARMA) (see also references [185,186] for the integrated, actual, and realized volatility as well as higher moments of volatility [185]; the multifractal analysis of financial markets is presented recently in reference [187]). Large volatilities-peaking at the time of big financial crashes [169]-imply rapid price chances and risky trading [173,188]. In particular, the future values of volatility of an asset are hard to assess: this leads to often imprecise predictions and necessitates modified BSM models [38,39,88,189] via including, e.g., 'implied volatility' concepts [190].
Specifically, the description of the 'Joseph', 'Noah', and 'Moses' effects-accounting, respectively, for longtime correlations (or nonindependence), nonfinite variance (or 'fat tails'), and time-dependence (or nonstationarity) in the distributions of increments of a stochastic process-for the FTS-analysis was presented in reference [176]. A pronounced nonstationarity of the variance of increments for the daily stock-market prices was also demonstrated [174,175]. The distinct, U-shape-like, stark volatility variations obtained from the high-frequency-data analysis-shown to be universal for a number of indices and currencies [174,175], being also positively correlated with the respective trading volumes-are inherent to the intraday volatility and stocks-trading activity. This time-dependent diffusion [174] is vital for understanding the rich and complicated dynamics of daily evolution of stock prices and trading volumes. Particularly high volatilities were also observed at the start of the morning trading sessions [175]. The related effects on longer time scales, with higher returns on Mondays and during the first weeks of January (connected to larger price variations), are also documented [138].
The 'diffusion coefficient' or volatility was also shown to depend on time and currency exchange rate, with several intervals during the day when the standard deviation of increments revealed distinct, multiple powerlaw-like scaling behaviors with time, also featuring drastic variations of the intraday-volatility magnitude [174,191,192]. This can gives rise to 'spurious' [174] non-Gaussianity in the distributions of price increments, including Laplace-like forms [193]. Such non-Gaussian features were ubiquitously observed [194,195] when averaging over longer time-scales was employed and simple/standard time-independent model parameters were used in the fitting analysis. Note that this complicated but repeatable intraday dynamics of the GBMmodel parameters favors ensemble-averaging-based analysis [192,196]-with a large set of daily trading data composing a statistical ensemble of (independent) realizations,-as compared to the methods of time averaging, such as those employed in the current study. The standard GBM model with constant drift and volatility should, thus, be essentially modified to account, e.g., for a power-law variation of volatility [174,191,192], σ(t) ∼ t γ , advocating for a new multiplicative stochastic process of 'scaled' GBM [167].

Author contributions
All the authors have made a substantial intellectual contribution to the work and approved it for publication.

Competing interests
The authors have no competing interests to declare.

Acknowledgments
A G C is deeply grateful to Humboldt University of Berlin and to Professor I M Sokolov personally for hospitality and support. A G C thanks V G Cherstvy and S Thapa for their help with some graphs and A Serafimovich for providing reference [198]. A G C thanks H Safdari, J Shin, S Thapa, and W Wang for their long-time help with providing scientific articles inaccessible from University of Potsdam. RM acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG Grants ME 1535/7-1 and ME 1535/12-1). RM also thanks the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej) for support within an Alexander von Humboldt Polish Honorary Research Scholarship.

Data availability statement
The data used in this study were accessed from https://finance.yahoo.com.

Appendix A. Supplementary figures
Here, we present some auxiliary figures supporting the claims in the main text.
Clearly, equation (B4) also follows from This yields the second moment of GBM as    straightforward two-parameter fit of the variation of versus t d for a given index or stock to find the values of drift and volatility values often yields unsatisfactory results as to the accuracy and consistency of the fit (results not shown). Such a fit is particularly problematic for the FTS encompassing a period of a financial crisis or severe price drops. A more rational way to examine the discrepancies between the theoretical GBM-based predictions and real-data results is to compute the sum

C1.2. 'Classical stocks'
The procedure (C2) gives one point in the {μ, σ} plane, see the example of such contour plot for the stockprices of BA in figure BB1. We find that for all the indices studied the model-versus-data discrepancy quantified by (C2) increases for small σ and small μ values (not shown). In the opposite domain, for large μ and large σ values, the discrepancy increases as well, although often less steeply. At intermediate values of μ and σ an 'optimal valley' featuring minimal model-versus-data deviations is often realized, as illustrated in figure BB1. Figure BB2. The same as in figure BB1 but for the three main cryptocurrencies. The maximal t d listed in the plots is the value used to compute the deviations in equation (C2). The log-ratios are always calculated from the full FTS available, while the maximal delay times or the corresponding maximal years in the analyzed data are fixed (the respective values are listed in the plots). The data are taken till the last day of the 'max years' denoted in the graphs.
We find that for the 'classical stocks' the position of this optimal valley in the {μ, σ}-plane often does not change strongly upon varying the terminal point of a given FTS and the delay time t d (a company-universal 'signature'). These plots are the 'fingerprints' reflecting certain functioning principles employed by a given company (in terms of the deterministic μ-based and stochastic σ-based features of resulting variations of its stocks price). The depth of the minimum in figure BB1 for a given valley is an index-specific feature; the depth along the valley often is only weakly sensitive to the variations of μ and σ along the 'valley' (results not shown).

C1.3. Cryptocurrencies
In figure BB2 the optimal {μ, σ}-valleys for the three main cryptocurrencies are shown. Their S(t)-dynamics is much more rapid and volatile as compared to the classical stocks: the valleys are significantly shifted toward larger values of drift and volatility. In this plot, the data up to the end of 2017 are only considered in the analysis because after this data a number of crashes on the cryptocurrency markets happened and, as a consequence, the results cannot be expected to follow the GBM model. The latter data were, thus, excluded from the analysis of optimal μ and σ values for the cryptocurrencies. This speculation-driven shift of the valley for the cryptocurrencies toward larger μ and σ is also clearly visible in figure BB3. Note that, as cryptocurrencies are traded 24/7, the delay times given in figures BB2-BB4 were recalculated (via multiplying by ≈252/365) to  make them comparable to those for classical stocks (traded for ≈252 days per year). One fiscal year and one t d -year contain, therefore, different numbers of calendar days for classical stocks and cryptocurrencies.
We find that, contrary to the classical and well-established stocks, as shown in figure BB1, for a highly speculative dynamics behind the price evolution of the cryptocurrencies the terminal date of their FTS entering the analysis of equation (C2) indeed changes the location of the optimal {μ, σ}-valley. Specifically, when the latest data for BitCoin or BTC are included in the analysis, the position of the optimal valley shifts toward lower drift and volatility values, see figure BB4. This likely is a manifestation of a 'cooling down' effect after the crash in late December 2017 (as well as several later drops of the BitCoin price in the period 2018-2020) that made cryptocurrencies generally less promising for speculative short-term 'in-and-out' investments during those periods. In contrast, the optimal valleys for the 'heating-up' phase in the dynamics of the cryptocurrencies (with the terminal years being 2015 and 2016 in figure BB4) are located in the range of considerably larger values of drift and volatility (artificial price 'heating-up' as a preparation for a later crash of the pyramid and to extraction of profits).
These features of the {μ, σ}-plots in figures BB2 and BB3 may be viewed as some speculative instruments enabling 'big players' to control, manipulate, and direct the market, while 'small players' can only react to trends/news. Such small investors act delayed, the disadvantage vital for highly volatile markets of cryptocurrencies. Additional (natural or artificial) delays in the BitCoin selling scheme-like those occurred during the super-hype with sky-rocketing BitCoin prices in late December 2017, with lag times of up to a week(!) [198] required to complete the transactions-further impede small investors from timely monetizing their profits.
The minimum in the optimal {μ, σ}-valley is often attained for an infinite set of μ-and-σ pairs, with a certain functional connectivity. This 'over-determination' or 'degeneracy' of parameters, not surprising for a two-parameter fit, suggests a better strategy, where one parameter is determined directly from the data (μ), while another one (σ) is then found from the computed 'optimal valleys', see appendix C2.

C2. Fitting δ 2 d,i (Δ) data and optimal {μ, σ}-valley
The value of μ can be found via direct linear fit of log[S i (t)] versus t for a given stock in the same time domain as used previously to calculate δ 2 d,i (Δ). In virtue of a roughly exponential price evolution of S(t) (if described by GBM and after neglecting the σ-containing term in (11)), one can estimate the effective drift μ σ,eff (36) and using the landscape of the 'optimal valley' find the company-specific value of volatility. In figure BB5 the linear fit of log[S(t)] for CAT is presented. Note here a considerable arbitrariness in the length of the FTS used  for the fit and certain uncertainties with regard, e.g., to the value of log[S(t)] rapidly changing in magnitude. This procedure gives reasonable estimates and errors for the optimal μ and σ for some of the FTS examined.
However, for the value of μ estimated this way no uniquely corresponding value of σ can sometimes be found from the {μ, σ}-contour plots. Moreover, the contour plots of GBM-versus-data residuals (C2) for all the indices studied (results not shown) indicate that estimating σ from the data first and finding the respective value of μ afterward can be a more appropriate procedure, see appendix C3. Note that some more 'local' estimates for μ of GBM can be obtained via fitting log[S(t)] not for the whole FTS-data available, but rather, e.g., only for a crash-free time-domain (provided the GBM model is applicable there). Lastly, one can in principle set μ = 0 and (for most indices) find the corresponding volatility via a one-parameter fit, see figure BB6 for CAT.

C3. Assessment of parameters using log-returns
A more strict and also previously used [199] procedure is to determine σ first using the concept of time-local log-returns. The latter, defined at discrete times t n along a given FTS S(t) as r(t n ) = log[S(t n )/S(t n−1 )], presents a standard quantifier of price fluctuations of a stocks index. The log in (C3) 'removes' the expected exponential growth of price in time (supposing S(t) develops according to GBM). For small price variations the log-returns (C3) approach the standard returns, r n . The cross-sectional standard deviation of returns defines the dispersion of individual returns with respect to the mean for a given period of time, the averaging window. During periods of market stress [145], the dispersion should increase due to higher volatilities of rapidly fluctuating prices, but the herding behavior has, in contrast, an opposite impact: people tend to follow the trends thus reducing the dispersion [145,200]. The average volatility σ for allN points used in the analysis of a given FTS is given through the trajectoryaverage return defined as by [199] We calculate σ from the daily FTS with Δt = 1/252 (for ≈252 trading days per year [for the classical stocks]). The annual volatility is then Clearly, via using (C5) one can also compute the volatility for each year of the FTS-data separately if, e.g., time-dependent annual features of the price-evolution dynamics are to be considered. Plotting the variation of the squared residuals (C2) versus μ for the volatility value found from (C5) for each of the FTS we easily find a unique value of the drift parameter μ, as illustrated in figure BB7 (due to the existence of a clear minimum for the curve). This method provides a direct, rational, and unique way to find the optimal pair of drift-and-volatility for a given FTS and in an arbitrary time-domain making the current method superior to the approaches discussed in appendices C1 and C2).
The parabolic-looking curves in figure BB7 represent the cross-sections of the two-dimensional plots of figure BB1 for the same index/company after one has set the volatility to σ = σ given by equation (C5) for each of the stocks. The parameters μ and σ are bound to the time unit used in data-based computation and simulations 11 . We observe, e.g., that a sharper minimum for IBM in figure BB7 symbolizes a smaller 'tolerance' 11 If time is in years (as in most of our simulations in section 4, or t = 1/252 being one day), then μ is the yearly ('per annum') interest rate and σ is the annual volatility (for classical stocks). Thus, if not stated otherwise, we always refer to the annual values of μ and σ to compare the outcomes to the results of GBM simulations. The volatility-defined as the standard deviation of log returns, distributed normally with N ((μ − σ 2 /2)t; σ 2 t) for GBM, equation (34)-scales with √ time window. To compute 'annual' volatility, we can calculate the standard deviation from the yearly log-returns for sample size N = 60 (for a 60 years-long FTS), or from the daily log-returns for the same data with the sample size N = 15 120 days and then scale the value with √ 252 (for 'classical' stocks), as in equation (C6). of this index (in a given time domain) to deviations of parameters. Namely, upon variation of the drift parameter the GBM-based predictions for the delayed TAMSD quickly deviate from the determined minimum. This sharp minimum corroborates the behavior of the squared residuals (C2) for IBM when visualized as the density plots in the {μ, σ}-plane (not shown).