Memory Effects, Multiple Time Scales and Local Stability in Langevin Models of the S&P500 Market Correlation

The analysis of market correlations is crucial for optimal portfolio selection of correlated assets, but their memory effects have often been neglected. In this work, we analyse the mean market correlation of the S&P500, which corresponds to the main market mode in principle component analysis. We fit a generalised Langevin equation (GLE) to the data whose memory kernel implies that there is a significant memory effect in the market correlation ranging back at least three trading weeks. The memory kernel improves the forecasting accuracy of the GLE compared to models without memory and hence, such a memory effect has to be taken into account for optimal portfolio selection to minimise risk or for predicting future correlations. Moreover, a Bayesian resilience estimation provides further evidence for non-Markovianity in the data and suggests the existence of a hidden slow time scale that operates on much slower times than the observed daily market data. Assuming that such a slow time scale exists, our work supports previous research on the existence of locally stable market states.


Introduction
The S&P500 is an aggregated index of the stocks of the 500 largest companies traded at US stock exchanges and therefore serves as an indicator of the overall US economic performance.Estimating and predicting the correlation between different assets is crucial for optimal portfolio selection and risk management and has been the focus of financial research since Markowitz's portfolio theory [27].
Because the economy is a system with a large number of interacting parts (traders and companies), it is amenable to the analysis tools from complex systems science [25].Due to the increasing availability of data for the economy, highly data-driven methods can be applied in this field of research, with many researchers choosing to focus on the correlations between the relative price changes of different stocks, i.e. the correlations of the stocks' returns [4,23,26,27,28,34].In [32,38,45] the authors identified states of the economy by analysing correlation matrices of daily S&P 500 data and found eight clusters.These can be interpreted as economic states reflecting e.g.times of economic growth or crises [32] and are found to be locally stable [38,45].Further analyses showed that exogenous events, precursors for crises and collective market behaviour can also be identified via the correlation matrix approach [14,13,15].Ever since Bachelier's seminal work introduced the random walk model to describe stock prices [2], the inherent stochasticity of financial data has been taken into account by researchers and practitioners alike.The Langevin equation is a model for stochastic differential equations that includes a deterministic drift and a random diffusion and can be used to describe such stochastic data.Much research has been devoted to estimating Langevin equations from data [7,43,35,6,36,20,54] and Langevin equations have found applications in various fields of research such as fluid dynamics [42], molecular dynamics [21] and meteorology [5] (cf.[8] for a review of applications).The generalised Langevin equation (GLE) expands the Langevin model by introducing a kernel function to include memory effects [31].Bayesian statistics includes a collection of methods that reverse the classical approach to statistics and focuses on calculating posterior distributions of model parameters as probabilities conditional on the observed data [44,49].This approach enables an efficient estimation of Non-Markovian Langevin equations [53,52].Also, other time series analysis methods can be implemented in the Bayesian frame work to e.g.detect change points in time series data [18].We show that an estimated GLE model manages to achieve a high goodness-of-fit for the correlation time series and implies a strong memory effect which has to be taken into account when predicting future market correlations.Furthermore, we perform a detailed comparison of a Markovian mono-time scale and a Non-Markovian two-time scale model with a hidden slow time scale which confirms the GLE analysis results regarding Non-Markovian memory effects.Additionally, the analysis supports the local quasi-stationary economic state theory discussed in Stepanov et al. [45] and provides some evidence that a hidden slow time scale might be involved in the mean market correlation dynamics.It is not far-fetched to assume that a complex system like the human economy involves a multitude of time scales and our findings coincide with such reasoning.The slow time scales might be connected to business and innovation cycles or similar economic dynamical mechanisms even though we could not extract a quantitatively reasonable magnitude of the hidden slow time scale.The remainder of this article is structured as follows: Section 2 explains the data gathering and preprocessing in 2.1, the Bayesian methodology in 2.2, the GLE estimation procedure in 2.3 and the resilience analysis in subsection 2.4.Section 3 describes the goodness-of-fit and the estimated memory parameters for the GLE model in 3.1 as well as the resilience results including the two-time scale discussion in section 3.2.Finally, the results of these analyses are interpreted and compared to the results from [45] in section 4.

Data Preparation
The S&P 500 is a stock index containing 500 large US companies traded at American stock exchanges.Daily stock data from these companies were downloaded via the Python package yfinance [1] for the time period between 1992 and 2012, which is the same time period as in [38].Only stock data of companies that were part of the S&P 500 index during at least 99.5% of the time were used for this analysis.If a company's price time series P t is not available for the full time period, the remaining 0.5% of the price time series data are interpolated linearly with the .interpolate()method in pandas [37,30].Overall, there are 249 companies' time series for 5291 trading days with one observation per date.Each company's returns R t = (P t+1 − P t )/P t are locally normalised to remedy the impact of sudden changes in the drift of the time series with the method introduced in [39] as Here, . . .n denotes a local mean across the n most recent data points, i.e. r t is subjected to a standard normalisation transformation with respect to the local mean and standard deviation (i.e.volatility σ).Following [32], n = 13 was chosen for the daily data.For each time step t and each pair of sectors i, j the local correlation coefficients are calculated over a time period of τ = 42 trading days like in [32] (the 42 working days correspond to 2 trading months) with the local standard deviations σ (i) τ .As shown via Principle Component Analysis in [45], the mean correlation already describes much of the variability in the data.Hence, it makes sense to restrict the analysis to this one-dimensional time series C(t) (shown in figure 1).The time t is here selected as the central value of the time window of length τ in order to have a symmetrical window.The preprocessed time series is available via [50].

Bayesian Statistics
Bayes' theorem for the conditional probability distributions f (•|•) of model parameters θ and observed data connects the standard statistical likelihood f (x|θ) and prior knowledge about the model parameters f prior to a posterior distribution f post of the unknown model parameters conditional on the observed data [24].A bundle of methods have been derived from this approach and are collectively referred to as Bayesian Statistics [44,49].Markov chain Monte Carlo algorithms (MCMC) can be used to infer the posterior distribution by simulating several random walkers that explore the (potentially high-dimensional) posterior distribution [11].It generates samples of the model parameters which are uncorrelated after cutting off the first n burn exploration steps as a burn-in period and thinning the remaining samples.The resulting MCMC samples can also be used to integrate the posterior distribution over all but one model parameter θ i to derive the marginal posterior distribution f (θ i |x) of this single parameter.Parameter estimation can then be done via the mean of the posterior distribution or via its maximum (abbreviated as MAP for maximum a posteriori) and credible intervals (CIs) of e.g.95% credibility can be directly derived from the quantiles of the samples for θ i .

Fitting a Generalised Langevin Equation with Memory Kernel
Sections 4 and 5 of [45] analyse the one-dimensional mean correlation time series by fitting a Langevin model d C with independent Gaussian noise Γ, a deterministic time-dependent drift function D (1) and a timeindependent diffusion function D (2) .Note that the Langevin equation is a Markovian model, i.e. it has no memory.As an alternative model, a Generalised Langevin Equation (GLE) includes previous values of the time series in a memory kernel K, but has only time-independent parameters with a functional equation A Bayesian estimation of the parameters of equation ( 6) is implemented in [53].It discretises the memory kernel K(k) to K k and discretises the observed data into n B different bins to significantly speed up the estimation procedure.Because an overlap of the windows used to calculate C can lead to artefacts in the memory effects (cf.section A in the appendix), we choose to calculate C with different parameters than τ = 42 and s = 1 from [45] (every s th is retained for the analysis and hence for s = 1, the full time series is used).Instead, we set τ = 5 and s = 5, meaning that the mean market correlation C is calculated over one trading week for the end of the trading week in disjoint windows (i.e.we create a time series whose length is 1/s = 20% of the original time series).Because using disjoint intervals automatically leads to a thinning of the time series, this seemed like a useful trade-off with a real-world interpretation as weekly correlation.Note that although only five trading days contribute to the calculation of C, there is still an order of magnitude of 10 4 correlations C i,j whose mean value is the desired C. I.e.although only a short time window is used to calculate each of the (i, j) correlation pairs, the sheer size of the ensemble of C i,j reaffirms our trust in the accuracy of their mean C. The resulting time series is shown in figure 1 together with the time series used in [45] and although it obviously becomes much noisier due to the shorter window size, it still retains some of the features of the less noisy time series like e.g. the position of spikes with high correlation.
As described in [53], the memory effects are aggregated to a quantity K k (named κ k in [53]) that measures the strength of all memory effects from 1 up to k time steps ago.If K k saturates towards a plateau at step k 0 , then k 0 is the maximum length of any reasonable memory kernel.This often also coincides with kernel values of K k0 ≈ 0 at k 0 .The estimated values for the model parameters are chosen by calculating the marginal posterior distribution and choosing the mean or MAP parameter estimation.The goodness-of-fit is also evaluated by using MCMC to simulate an artificial time series C(a) t with the estimated model and comparing the autocorrelation function and the distributions of C(a) t and the increments C(a) t − C(a) t−j to those of the original data.Finally, the inferred model structure can be trained on only a subset of the data (training data) and be used to predict the remaining test data to evaluate its predictive performance.

Resilience Estimation
By applying Bayes' theorem (4) we can deduce a quantitative measure of resilience under given model assumptions.Therefore, we infer the parameters of two models of stochastic differential equations in a rolling window approach that allows for resolving the time evolution of the resilience and noise level and accounts for the quasi-stationary nature of the time series that is observed in Stepanov et al. [45].Following the quasi-stationarity argument we assume to be in a fixed state per window, i.e. the fixed point C * which is approximated by averaging over the mean market correlation data per window.First, we estimate a Markovian Langevin equation ( 5) with the Taylor-expanded drift with the drift parameters θ 0,...,3 (t; C * ) ≡ θ 0,...,3 ( C * ) per rolling window and the constant diffusion We choose uninformed priors which are given by for the linear part of the drift function and the Jeffreys scale prior [49] p prior (θ 4 ) = 1 θ 4 (10) for the noise level θ 4 with suitable prior ranges.The priors of higher order parameters are chosen to be p prior (θ 2 ) = N (µ = 0, σ = 4) and ( 11) with Gaussian distributions N centred around the mean µ = 0 with a standard deviation σ = 4 and σ = 8.Since we consider economic systems to operate on multiple time scales which concurrently lead often to Non-Markovian time series, we additionally introduce a two-dimensional Non-Markovian model analogue to Willers and Kamps [54].The model takes the form with a hidden Ornstein-Uhlenbeck process (OU-process) λ, drift D (1) • λ and diffusion C ( C, t) and diffusion D C ( C, t) of the observed process C(t) remain unchanged.The Non-Markovian analogue to the constant noise level σ 2 of the Langevin equation is given by the composite noise level with the small discrete sampling time step h.
For the OU-process, an invariant prior of a straight line and a scale prior for the diffusion are multiplied: Furthermore, via the prior we introduce a pre-defined time scale separation of the time scales τ C and τ λ of the observed and unobserved process, respectively, i.e. we require either τ C > γ • τ λ or τ λ > γ •τ C with a scale separation coefficient γ.The characteristic time scales [46] are approximated by The priors for the model of the observed data C(t) remain unchanged as well apart from the term D C ( C, t) = θ 2 4 which corresponds to a coupling constant in the Non-Markovian model.We simply apply a Gaussian prior like the one in equation (11) to θ 4 in this case.Inspired by the formalism of linear stability analysis for both models we calculate the drift slope per window as a Bayesian parametric resilience measure.For the Markovian model the calculations are performed with the open-source Python toolkit antiCPy [16,17].In this modelling framework, a stable state corresponds to a negative drift slope ζ, whereas a changing sign indicates destabilisation via a bifurcation.More details on the procedure can be found in [18].3 Results

Estimated GLE Model
The data is split into n B = 10 equally wide bins and an initial modelling attempt with a rather long kernel k max = 10 is tested, i.e.K q = 0 for all q > k max .It shows a plateau emerging at around k ≥ 6 (cf.appendix B).Hence, a model with k max = 6 is used as a reasonable length for the memory kernel and its goodness-of-fit is evaluated.We then use a very conservative estimation of the memory up to k max = 3 to evaluate the GLE's predictive power.

Goodness-of-fit
A time series of length 10 5 is simulated via Euler-Maruyama integration of the GLE model to compute the autocorrelation function (ACF) and to compare it to the ACF of the original time series.The best model is chosen via MAP estimation in the Bayesian framework of [53], but selecting the mean estimation yields almost identical results.Both ACFs are computed via the function statsmodels.tsa.stattools.acffrom the Python package statsmodels [41].Figure 2 shows that the two ACFs show very good agreement up to lags of 10 trading weeks and decent agreement up to lags of 20 weeks.The distributions of the time series and the first two increments for the original and the simulated data are shown in figure 3 for the MAP estimation and also for the mean estimation.Figure 3 shows an almost perfect overlap between the two estimation procedures and a good agreement between the estimated time series and the original time series, especially for the increment distributions.Overall, these diagnostics indicate that the estimated model with memory kernel length 6 manages to reproduce these important statistical properties of the original time series.Because the realised values of the time series C are not distributed uniformly, we also use a modelling procedure with unequal bin widths so that each of the 10 bins has the same amount of data.However, this barely changes the model diagnostics shown in figures 2 and 3.The only noticeable change is that for long kernels with length ≥ 10, the ACF seems to be captured a little worse in the shown range of figure 2, but manages to fit more closely to the empirical ACF for large lags at around r ≈ 100.While the regular LE without memory has a very similar increment distribution, the ACF is much worse than for the GLE.

Estimated Memory Kernel
To make further inference on the memory kernel and to estimate its quantitative effect, the posterior distribution of the model with memory length 6 is sampled via MCMC.100 walkers are simulated for 10 5 time steps and after a burn-in period of 450 initial time steps is discarded, the chains are thinned by only keeping every 450 th step to obtain uncorrelated samples.Bayesian CIs can now be calculated at the 95% level for each parameter in the kernel function and the results are shown in figure 4. The CIs of K 4 are already very close to the zero line and those of K 5 include zero, meaning that there is no evidence for a memory term at five weeks distance.Interestingly, the CIs for K 6 clearly exclude the value zero.Because it is difficult to exactly identify the beginning of the plateau in the memory aggregation in figure 9, it may be possible that the plateau already emerges at k = 5 and that the nonzero memory kernel K 6 is therefore unreliable.However, the results in figure 4 clearly imply a nonzero memory effect for four time steps with a clearly nonzero effect strength for memories up to k = 3 trading weeks.combined with the uncertainty about the beginning of the plateau in figure 9 imply that the nonzero memory effect for K 6 may be misleading.All previous memory kernels K 1,2,3,4 have nonzero values at the 95% credibility level.

Prediction via the GLE with Kernel Length 3
With a conservative interpretation of the results in section 3.1.2,a model with memory length k = 3 is used to evaluate the GLE's power to predict a future value y t+1 by forecasting an accurate prediction ŷt+1 .It is tested against a regular Langevin equation (LE) model without any memory effects (corresponding to k = 0 and also estimated with the code in [53]) and against the naive benchmark of predicting the next time step y t+1 by simply setting it to the last previously known value: ŷt+1 = y t .To test these three methods, the Langevin models are trained on the first α% of the time series (the training data) and the predictions of the GLE, the LE and the naive forecast are evaluated on both the training data (as in-sample predictions) and on the remaining 1−α% test data (as out-of-sample predictions).The coefficient of prediction ρ 2 is used to evaluate their predictive accuracy.If n observations y 1,...,n are forecasted as ŷ1,...,n , then the coefficient of prediction is given by with the ȳ denoting the mean value.It takes the value of ρ 2 = 1, if the prediction is always exactly true, and ρ 2 = 0, if the prediction is only as accurate as always using the mean ȳ, and ρ 2 < 0, if it is less accurate than using the mean and is computed via the function sklearn.metrics.r2score from [33].The results in table 1 show that the GLE model consistently achieves the highest accuracy on in-sample and out-of-sample predictions for the three chosen test data sizes.Notably, the negative ρ 2 of the naive method for the test data indicates that the out-of-sample prediction is by no means trivial, meaning that the low, but positive ρ 2 of the GLE on the test data is nevertheless a good performance.The GLE achieves slightly better results than the LE, indicating that the memory effect should be taken into account for prediction tasks.Figure 5 shows the predictions of the LE and GLE for the in-sample and out-of-sample predictions with α = 90%.The same visual comparison between naive forecast and the GLE can be found in the appendix D.

Hidden Slow Time Scale and Non-Markovianity
Applying the resilience analysis method, described in subsection 2.4, we can deduce some interesting evidence for multiple time scales and further confirmation of Non-Markovianity present in the considered economic time series.The results for both the simple Markovian and multi-scale Non-Markovian model (cf.( 5) and ( 13), respectively) are compared in figure 6 (a-c).For better readability the parameters of the calculations are listed in table 2 of appendix E.
Stepanov et al. [45] argue for quasi-stationary economic states which are occupied over finite time periods, before they transition into another quasi-stationary economic state.The approximated state potentials from the mean market correlation in the article [45] suggest that there might be shifts in the fixed point positions over time, but no bifurcation-induced tipping (B-tipping) is involved, i.e. no qualitative change of stable and unstable fixed points or attractors is observed.Instead of B-tipping mechanisms, the intrinsic economic stochasticity drives the jumps between alternative quasi-stationary states, which is a mechanism basically related to noise-induced tipping (N-tipping).
If we choose a model that captures the key features of the data generating process, we should thus be able to uncover generally negative drift slope estimates ζ, corresponding to data of a locally quasi-stationary state in each rolling window of the mean market correlation C(t) the data of which is shown again in figure 6 (a).In this spirit, we find the Non-Markovian model with an unobserved slow time scale, i.e. τ λ > γ • τ C with γ = 2, to yield the expected result of negative drift slope estimates ζ as presented in figure 6 (b) and indicated by the red solid line with green credibility bands (CBs).This result is supported by rather intuitive qualitative considerations: Economic processes are wellknown to operate on various fast and slow time scales, on the one hand e.g.high frequency trading, trading psychology and sudden external events, whether that might be political, sociological, severe weather or other impacts.These fast-scale processes could be termed as "economic weather" in a climatologist's metaphor.On the other hand, long-term "economic climate evolution" takes place on much slower time-scales.Examples could be innovation processes and technical revolutions like the invention of the steam engine or the internet, economic cycle theories like that of Marx (τ ∼ 2 ato10 a) [29,10], Keynes, Schumpeter or Kondratjew [3,22,40], cycles of fiscal and demographic developments [48], cultural evolution [51] and generation changes influencing economic reasoning, adaptions to climate change and the scarcity of resources and much more.Keeping that in mind, the trading day resolution of the data is rather fine-grained and it might be reasonable to assume that a hidden slow time scale is present in the data.The slow time scale of the presented Non-Markov model is determined by the upper boundary of the prior range which is chosen to correspond roughly to 6 a which is the time of an average business cycle and coincides with the first zero-crossing of the autocorrelation function of the mean market correlation C(t).If the prior is chosen broader the model converges to roughly 700 ato4000 a which is not a plausible magnitude of the time scale of economic evolution (cf.Appendix E).However, the main results of local quasistationary economic states is not affected by the prior choice.The estimation of a more reasonable magnitude of the slow time scale without prior restriction might be prohibited by the very limited amount of data per window, since the window range does not even include one complete business cycle which would be the smallest proposed slow time scale candidate.Further note that the twotime scale Non-Markovian model in eq. ( 13) incorporates noise in the slow non-observed variable λ.In that way it could formally reflect intermediate dynamics on a time scale τ N with τ C < τ N < τ λ .But since it is coupled to the mean market correlation C(t) on trading day resolution, we consider  6 (c) following the same color-coding.The noise level seems to increase over the years with a clear increase in the periods of the global financial crisis and the Euro crisis which accounts for a higher N-tipping probability in this highly turbulent economic period.The noise plateau of roughly one window length around the end of the Asian and the beginning of the Russian financial crises around 1998 is probably the results of the outliers that are incorporated into the windows.
Since the discussed observations alone only allow for relatively weak qualitative deduction of the discussed features of multi-scaling and Non-Markovianity, we additionally perform an analogous analysis on a synthetic time series x that shares the key features of a hidden slow time scale and Non-Markovianity with the original one.In figure 7 we provide a comparison of the per definition stationary first differences, (a) of the original mean market correlation C(t) and (b), of the synthetic time series x.The noise level in x is assumed to increase over time to mirror the noise level evolution of the mean market correlation C(t), suggested by the estimates in figure 6 (c), and is adjusted to cover almost the range of the first differences in C(t).Only the positive trend of the mean market correlation visible in figure 6 (a) is not included in the simulations of x.In that way the PDFs of the two time series' first differences are shaped similar apart from the fact that the highly centered probability mass of the mean market correlation C(t) with steep tails due to rare outliers is a bit more smeared out into flatter tails of the synthetic time series x.More simulation details can be found in the Appendix E.
The resilience analyses on the synthetic time series x are shown in figure 6 (d-f) and are in very good agreement to the original analyses results in figure 6 (a-c).That appears at least as an independent qualitative confirmation of our findings' interpretation.Moreover, this interpretation is supported by two additional facts: First, the estimation of a Markovian model on the Gaussian kernel detrended version of the mean market correlation C(t) results in similar results to the multi-scale Non-Markovian model.This strengthens our previous findings, because the detrending subtracts a non-stationary slow process suspected to be present in the data.We notice that weaker detrending leads to positive trends in the drift slope estimates ζ which could be due to increasing distortions due to incomplete detrending of the non-stationarity.Second, we fit a multi-scale Non-Markovian model with inverse time scale separation τ C > γ • τ λ with γ = 2 (i.e. the observed trading day time scale of the mean market correlation C(t) is considered to be at least two times slower than the hidden time scale).This leads to results similar to the Markovian model with and without detrending.This is an expected result, since the prior restriction of the time scale separation basically restricts the model to the Markovian case in which only the trading day time scale can be resolved apart from an even faster stochastic contribution.In other words, the prior assumption τ C > γ • τ λ with γ = 2 prohibits the incorporation of an unobserved slower time scale even if it is present in the data.The discussed results are presented in more detail in Appendix E. The FD capture essentially the same range as the FDs of C(t), but are less centered arount zero.The time series is modeled with increasing coupling strength to imitate the increase of the noise level found for the mean market correlation C(t) in figure 6 (c).(c) The FD PDFs of both time series are comparable.The PDFs only deviate to some extent in the flatter tails of the synthetic orange histogram with less dense probability density around zero.The positive trend of C(t) is not modelled in the simulated time series x.

Discussion and Conclusion
The estimated GLE model manages to reproduce the statistical properties of the original data for the end-of-week correlations as shown in figs. 2 and 3.The estimated memory kernel parameters show that even with a highly conservative interpretation of the 95% credible level, there are clearly nonzero memory effects for memory terms for all lags as far back as a lag of 3 weeks.Therefore, it is advised to use a model with memory to describe the correlation of the S&P500 market, which is an improvement to the Markovian Langevin model estimated in [45].Moreover, the GLE estimation presented in this article achieves a high goodness-of-fit for the entire time series, whereas Stepanov et al. used a time-dependent Langevin model by splitting the time series into different intervals and estimating Markovian Langevin equations for each of them.Our work shows that this procedure can be circumvented by using a model with memory of at least 3 trading weeks.The major advantage of our method is its possible application in predicting future market correlation: The time-dependent drift estimation in [45] has no clear or smooth functional dependence on time and therefore, little information can be inferred about future values of the correlation time series.Our method needs no time dependency, generalises over the entire time series and can be used to predict future correlation values which can be used for portfolio risk assessment.As shown in section 3.1.3,the memory kernel helps the GLE to achieve better prediction accuracy than the regular Langevin equation and much better results than the naive forecasting method of using the last observation as the predicted value.Notably, the existence of memory effects in the market's correlation structure can be interpreted in the context of volatility clustering.It is a well-known stylised fact from empirical research on financial markets that the volatility of a stock's returns tends to cluster: periods of high volatility are often followed by periods of high volatility and vice versa for low volatility [9,12].The correlation ρ X,Y between two asset returns r X and r Y with expectation values µ X , µ Y and volatilities σ X , σ Y with is defined as and therefore directly includes the volatility values σ X and σ Y .Because the time series of volatility estimators σ X (t) shows a well-known memory effect, it is not far-fetched to assume a similar memory effect in the correlations ρ X,Y between two assets or, as we have discussed in this article, in the mean correlation of the market as a whole.
These considerations are complemented by a resilience analysis that involves the estimation of a mono-time scale Markovian model and a two-time scale Non-Markovian model.In contrast to the memoryless Markovian model, only the Non-Markovian model exhibits the negative drift slopes which are in line with the hypothesis of locally quasi-stationary economic states postulated and observed in Stepanov et al. [45].An independent change point analysis approach also supports this view [19].Overall, these findings provide new evidence for the existence of such locally quasistationary economic states and for the presence of a significant non-Markovian memory effect.Interestingly, the resilience analysis yields some evidence that a second time scale which is slower than the trading day time scale of the mean market correlation data, is involved in the underlying economic dynamics.Economic processes operate on various fast and slow time scales, e.g.day trading, trading psychology and sudden external events -whether that might be political, sociologic, severe weather or other impacts -may be incorporated in the fast trading day resolution of the mean market correlation data.We refer to these fast-scale processes in terms of "economic weather" to employ a metaphor from climatology (cf.also the distinction in climate-like and weather-like tasks in [47]).In contrast, the long-term evolution of the "economic climate" might involve innovation processes and technical revolutions like the invention of the steam engine or the internet, economic cycle theories like that of Marx (τ ∼ 2 ato10 a), Keynes, Schumpeter or Kondratjew [29,10,3,40,22,48], cultural evolution [51] and generational changes influencing economic reasoning, adaptions to climate change and the scarcity of resources and much more.However, we were not able to derive an economically reasonable magnitude of the slow time scale which we would expect to lie in the range of decades up to hundred years corresponding to well-known economic cycle theories or cultural evolution processes.Instead, our applied MCMC model estimation without prior range restriction converges to a hidden slow time scale of roughly 700 ato4000 a.Nevertheless, our results suggest that there should be involved at least two time scales in the data-generating process which is an interesting starting point for future research.Notably, it is not particularly surprising that we could not quantify the hidden time scale, since we employ a very simple model parametrisation, have only access to one variable of the high-dimensional economic state space and perform our estimation on small windows that not even include the smallest economic cycle time scale of roughly 2 ato10 a that typically correspond to business cycles.Against this background it might be a very interesting challenge of future research to develop more realistic models and estimation procedures that perform reliably under the circumstances of limited data per window and incomplete variable sets to uncover the manifold of hidden time scales in the complex system of human economy.The distributions look very similar to the ones in figure 3.

D GLE: Prediction
The results of the naive forecasting method are compared to the GLE in this figure, showcasing the superior accuracy of the GLE, which can also be seen by comparing the ρ 2 values in table 1.

E.3 Parameter Documentation
In this subsection we summarise the parameter sets involved in the resilience calculations of the main article and the appendices in order to guarantee reproducibility without affecting the readability of the main article.The parameters are listed in table 2 and prior configurations are given in table 3.

Figure 1 :
Figure 1: The mean correlation of the S&P500 market calculated.The time series are the mean values of correlation matrices which were calculated based on moving windows of length τ days plotted against the centre of the τ -days-interval.This figure depicts τ = 5 and τ = 42.

Figure 2 :
Figure 2: Comparison between the autocorrelation functions of the original time series (solid line) and the simulated data from the model with kernel length 6 (dashed line) simulated via MAP estimation of the MCMC-integrated marginal densities.Up to lag 30, the ACFs show good agreement.The alternative simulation via mean estimation of the density yields an almost identical ACF.However, the ACF based on the model with no memory kernel (dotted line) fails to capture the empirical ACF.

Figure 3 :
Figure3: KDE plots to compare the distributions between the real data and the MCMC samples for the two estimated GLE models with kernel length 6 for the time series data Ct (left), the increments Ct − Ct−1 (centre) and Ct − Ct−2 (right).Differences between the two simulated models are hardly visible and overall, there is a good overlap with the real data.For the LE without memory, the respective distributions are depicted in appendix C.

Figure 4 :
Figure4: Values for the memory kernel K k in the model with memory length 6. Credible intervals were estimated via MCMC at the 95% level.The inclusion of 0 in the credible interval of K 5 combined with the uncertainty about the beginning of the plateau in figure9imply that the nonzero memory effect for K 6 may be misleading.All previous memory kernels K 1,2,3,4 have nonzero values at the 95% credibility level.

Figure 5 :
Figure 5: Comparison between one-step-ahead forecasts of the GLE model with Kernel against a regular Langevin estimation without memory effects on the training data (left) and the test data (right).The identity f (x) = x is given as a benchmark for perfect predictive accuracy.Here, α = 90% of the time series were used as training data.

Figure 6 :
Figure 6: Results of the resilience analyses based on the Markovian Langevin equation 5 and the Non-Markovian model 13 with an unobserved slow time scale τ λ , i.e. under the prior assumption τ λ > 2 • τ C : (a) For the mean market correlation C(t).(d) For a synthetic dataset x that shares the multi-scaling features and Non-Markovianity with the economic time series C(t).(b) The Markovian model suggests persistent latent instability, i.e. ζ ≈ 0, whereas the Non-Markovian slopes agree with locally quasi-stationary economic states as observed in Stepanov et al. [45].(c)The noise level estimates almost identical and tend to increase over time which might reflect more turbulent economic times towards the the financial crisis 2008.The plateau of ca.one rolling window length around 1998 might be due to the incorporation of some outliers in the time of the ending Asian crisis and the beginning of the Russian financial crisis.(e) Drift slope results of the synthetic dataset x analogue to (b).The qualitative results are in good agreement to (b) which might be a hint that the mean market correlation C(t) exhibits Non-Markovian features and is additionally governed by processes on a slower time scale τ λ .(f) The estimated noise levels of the synthetic dataset x are almost identical which is also observed for the noise levels of the mean market correlation C(t) in (c).For more details see the running text.

Figure 7 :
Figure 7: First differences (FD) and corresponding PDFs of the raw mean market correlation data C(t) and the synthetic dataset x which shares the key features of Non-Markovianity and an unobserved slow time scale with the mean market correlation data C(t).(a) The FDs of the mean market correlation are highly centered around zero with some outliers.(b)The FD capture essentially the same range as the FDs of C(t), but are less centered arount zero.The time series is modeled with increasing coupling strength to imitate the increase of the noise level found for the mean market correlation C(t) in figure6 (c).(c) The FD PDFs of both time series are comparable.The PDFs only deviate to some extent in the flatter tails of the synthetic orange histogram with less dense probability density around zero.The positive trend of C(t) is not modelled in the simulated time series x.

Figure 10 :
Figure 10: The same plots as in figure 3 for the Langevin Equation without any memory kernel.The distributions look very similar to the ones in figure 3.

Table 1 :
Comparing the predictions of the naive forecast ŷt+1 = y t , the Langevin equation (LE) and the GLE with memory kernel length k = 3 via their ρ 2 score on in-sample training data and out-of-sample test data (α% training data and the last 1 − α% of the time series as test data).Note that the partially negative ρ 2 for the test data indicates that the test data is quite difficult to predict, which corresponds to the high fluctuations in the final part of the time series in figure1.

Table 2 :
Rolling window and MCMC parameters of the resilience calculations.

Table 3 :
Prior configurations of the resilience calculations.TSS is an abbreviation for Time Scale Separation.