{'fraction_non_alphanumeric': 0.03338050059152381, 'fraction_numerical': 0.009812336171101162, 'mean_word_length': 4.757211538461538, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-based prediction and causality inference of nonlinear dynamics

Natural systems are typically nonlinear and complex, and it is of great interest to be able to reconstruct a system in order to understand its mechanism, which can not only recover nonlinear behaviors but also predict future dynamics. Due to the advances of modern technology, big data becomes increasingly accessible and consequently the problem of reconstructing systems from measured data or time series plays a central role in many scientific disciplines. In recent decades, nonlinear methods rooted in state space reconstruction have been developed, and they do not assume any model equations but can recover the dynamics purely from the measured time series data. In this review, the development of state space reconstruction techniques will be introduced and the recent advances in systems prediction and causality inference using state space reconstruction will be presented. Particularly, the cutting-edge method to deal with short-term time series data will be focused. Finally, the advantages as well as the remaining problems in this field are discussed.


Introduction
In the big data era, various types of time-series data are increasingly accumulated every day in various disciplines and in various scopes [14,18,55,75,87]. An outstanding problem is how to make systemlevel study based on the vast amount of data. Particularly, based on the measured/observed time series data, how to reconstruct the system and further study the mechanism behind the phenomenon has become an attracting problem and a hot topic in interdisciplinary science. Despite of the widely used approaches employing linear models [29], which are mainly based on correlation and regression, the natural systems usually show nonlinear nature: the components within one system intertwine with each other, and the dynamics of each component cannot be separated from one another. Such characteristic property distinguishes nonlinear systems from the linear stochastic systems, and nonlinear methods are thus more favorable in reconstructing complex systems.
Since the embedding theory proposed by Takens in 1980s [93], the State Space Reconstruction (SSR) technique has been developed and applied in nonlinear time-series analysis. The key idea of the SSR stems from the non-separability of the nonlinear system. The property of non-separability implies that the dynamics of one observable variable could contain all the dynamical information of the whole system. As a model-free method, the SSR based methods have been attracting persistent interest and under rapid development in the recent decades [30,42]. In this review, we will not cover the whole topics of system reconstructions using SSR methods, which can be referred to some excellent reviews [30,34,42,82,97], but focus on the forefront of the recent advances in predicting future dynamics and inferring causality using SSR framework. Particularly, the majority of SSR based methods usually require sufficiently long time-series data or a large number of samples, while the real data sets always have limited length due to the experimental restrictions or costs. We will introduce in detail the recently developed methods aiming to deal with short-term time series data.
This review is organized in the following manner. In Section 2, we first introduce the principle of various embedding theorems, which constitute the theoretical foundation of SSR. Based on the embedding theory, the framework of SSR as well as several related aspects in real applications are presented. In Section 3, the developments of predicting nonlinear system using SSR in the recent decades are reviewed. Particularly, how to make predictions from short-term high-dimensional time series data will be focused and introduced in details. Section 4 presents the directional interaction detection methods, i.e., the causality inference methods using SSR framework. Both the recently proposed Convergent Cross Mapping (CCM) method, using convergent of mutual predictions, and Cross Map Smoothness (CMS) method, particularly designed for short-term data, are introduced in details. We also speculate on several open problems, for instance, conditional causality and causal relations in delayed systems in this section. Conclusions and discussions can be found in Section 5.

Embedding theory
Consider a general dynamical system in the form of dx dt = F (x), (2.1) where x = [x 1 , x 2 , . . . , x n ] T is the system's state vector and F = [f 1 , f 2 , . . . , f n ] T with each f i (·) a nonlinear function. Starting from the initial value x(0) = x 0 , the trajectory x(t) describes the evolution of the system along time. After transient dynamics, the system's bounded trajectory will converge to an attractor, i.e., a bounded and invariant set A. Here the diverging situations are not under consideration. For a discrete system, the dynamical system is in the form of where F is a nonlinear map, and x(t), t = 0, 1, 2, . . . is the trajectory of the system, and the attractor can be defined analogously.
We consider the dynamics constrained in a compact manifold M containing the attractor A. For both continuous and discrete cases, the time evolution corresponding with an initial value x 0 ∈ M can be denoted uniformly by x(t) = F t (x 0 ): in the discrete case time t ∈ N and F t = (F ) t , while in the continuous case time t ∈ R and F t (x 0 ) is the integral curve through x 0 .
Intuitively, to study the property of the attractor in the state space, the values of all the system's state components x i (t), i = 1, 2, . . . , n along time should be known. However, the embedding theory implies that from a few or even one single time series observed from the system, it is possible to study the dynamical property of the whole system, i.e., the system can be reconstructed by observations. Actually, .. it was first discovered by Packard et al that using only one single coordinate of Rössler system one can reconstruct a faithful state-space representation of the dynamics in the original 3-dimensional space [68].
To be specific, the dynamical system of interest in [68] can be expressed in the form of Eq.(2.1) as and with time series of all the state components, the system's attractor can be obtained, as shown in Fig. 1(a) under x − y − z coordinates. As suggested by Packard et al, by sampling the time series of only single coordinate x(t), a variety of three independent coordinate sets can be obtained such as where τ is a sampling interval andẋ,ẍ are the first order and second order differences respectively. The attractors in such coordinates are shown in Fig. 1(b)(c) respectively. The geometry property of the attractors in Fig. 1(b)(c) are quite similar to that in Fig. 1(a). As a matter of fact, the dimensionality and the positive characteristic exponent can be further computed for each of these attractors [68] and the results confirm that the attractors obtained in Fig. 1(b)(c) can be regarded as reconstructions of the attractor obtained in Fig. 1(a). Such discovery was theoretically confirmed by Takens in his seminal work [93], in which the famous Takens embedding theory was proposed.
is an embedding where the integer L = 2d + 1.
Here an embedding of M means a smooth diffeomorphism from M onto its image ψ(M), i.e., a smooth one-to-one map which has a smooth inverse. The number L = 2d + 1 is the dimension of the space for the reconstructed attractor.
Actually, the embedding theory could be stated in other forms, such as derivative embedding and generalized embedding theorems [17,93,99] stated as follows.
is an embedding where the integer L = 2d + 1.
is an embedding.
It is clear that here the dimension d is important for the embedding theory. However, the dimension d is usually much higher than the essential dimension of the attractor A. Therefore, the embedding theory was later generalized into fractal dimensions by Sauer et al in [78]. Let d A be the box-counting dimension of the attractor A, which could be a non-integer and could be much smaller than the integer dimension d of the manifold M containing A, then for almost every smooth measurement function h, the maps defined in the Takens embedding theorems are still an embedding on each compact subset of a smooth manifold contained within A with L > 2d A .
Particularly, starting from one dimensional observed time series, a delay coordinate map could be introduced.
According to the fractal embedding theorem, the delay coordinate map is generically an embedding as long as L > 2d A . Thus it is straightforward to see that when the observation function h is simply chosen as coordinate function x i , the delay coordinate map with sampling interval τ is generically an embedding.
The one-to-one property of embedding is essential to the reconstruction of the original system because the state of a deterministic dynamical system, as well as its future evolution, are completely specified by a point in the state space. Therefore, the delay coordinate map actually sets up a bridge between the dynamics of the original system and the dynamics in the reconstructed space. Thus the delay coordinate is a state space reconstruction of the original system, which implies that the one dimensional observed time series can contain all the information of the dynamics of the whole original system. Moreover, the embedding theory implies that the dynamics in the reconstructed space is topologically conjugated with the dynamics in the original state space [42]. Thus, the nonlinear dynamic properties such as Lyapunov exponent, entropy, and even the unstable periodic orbit (UPO) are preserved from the original system to the reconstructed system, which forms the cornerstone of the state space reconstruction technique.

Issues related to state space reconstruction
When applying delay coordinate SSR to real data sets, the information of the original nonlinear system is generally unknown a priori, thus there are two key factors to be determined for successful SSR application: one is the embedding dimension L and the other is the delay τ . Dozens of methods for estimating optimal dimension L and delay τ have been developed in the past few decades based on time series data [8,21,25,43,44,46,47,49,50,69,77,94].
Theoretically, any τ > 0 could be used in the delay coordinate SSR process, but in practice, the delay coordinates will be strongly correlated with too small τ while the reconstructed attractor can fold over on itself with too large τ . Hence, one wants the smallest τ that maximizes the independence of the coordinates of the embedding vector. Therefore, several representative methods have been proposed to choose optimal τ by computing a statistic that measures the independence of τ −separated points in the time series. As τ increases, the first zero of the autocorrelation function of the time series, the first minima of the average mutual information, or that of the correlation sum [21,25,50] could be candidates of optimal τ . Actually, such τ maximizes general forms of independence, and consequently can be adopted in the delay coordinate.
Analogous to the case of choosing τ , the embedding dimension L should be neither too small nor too large. In theory, L should be sufficiently large so as to satisfy the condition of L > 2d A , however in practice, low dimensional reconstruction is more favorable. On the other hand, it is noticed that L > 2d A is a sufficient but not necessary condition. Thus the optimal L should be chosen as the smallest value that affords a topologically correct result. To this, the false near neighbor (FNN) algorithm by Kennel et al [44] is a feasible method generally accepted. In an FNN-based method, one tries to construct a delay coordinate embedding with embedding dimension k, computes each point's near neighbors, and repeats such process with increased embedding dimension k + 1. If any of the geometry relations changes, i.e., some neighbor in k dimension is no longer a neighbor in k + 1 dimension, then it is a false near neighbor and it indicates that the dynamics are not properly unfolded under embedding dimension L = k. Thus, the value appearing on the elbow of the decreasing FNN curve as k increases can be a candidate for optimal embedding dimension.
In the recent years, it has been noted that the product L * τ is the factor that really matters. Thus the both factors can be estimated at the same time and a unified approach has been proposed [69].
Furthermore, when applying SSR to practical measured time series, noise is unavoidable. The theory of state space reconstruction in the presence of noise was first considered in [11], which incorporated the observational noise and estimation error into the framework. Delay coordinate embedding theory for stochastic systems was later systematically studied by Stark et al in [85] and [86].

Prediction of nonlinear dynamics
Predicting nonlinear dynamics using state space reconstruction method has a long history, which can be dated back to 1969 when Lorenz proposed the method of analogues in the state space to forecast atmosphere [57]. Later the idea was shown to be applicable in the reconstructed state spaces by Pikovsky [71]. The nonlinear but smooth vector field or map in the state space implies that the evolution trajectories starting from two close states will still be close to each other in short time, and thus the information of unknown future value can be approximated by the future values of its near neighbors, as sketched in Fig. 2(a). The embedding theory guarantees that such property also holds for states in delay coordinate reconstructed state space. This has inspired the development of making predictions for nonlinear time series using nearest neighbors.

Nearest neighbors methods
The fact that close states will evolve to close future states implies that using near neighbors on the attractor can patch local information to approximate the true future values. Based on such local patching idea, a large number of creative strategies have been fruitfully proposed to predict future dynamics, including several influential works [10,19,91]. The research on this topic has even spawned a competition on nonlinear time series prediction [98].
To be specific, considering a state x(t) ∈ R d on the attractor A and its k nearest neighbors x(t i ), i = 1, 2, . . . , k where t i < t, the future value x(t + τ ) can be locally obtained by the values of x(t i + τ ). In [19], the zeroth-order approximation method was proposed where x(t + τ ) was simply approximated by where x(t ′ ) was the nearest neighbor of x(t) on the attractor A. In [10], the relation was assumed to be a polynomial and fitted by the k-nearest neighbors The simplex method was proposed in [91], where x(t + τ ) was approximated by the weighted average of its nearest neighbors' future values in the form of where ω i was set as an exponential weight to its original distance x(t)−x(t i ) and the number of nearest neighbors k was generally set to be larger than the reconstructed dimension d such that k > d. When the future value x(t + τ ) is obtained, the prediction of each component x i (t + τ ) is realized. Particularly, in the case of delay coordinate reconstruction, with the time series data of one dimensional observed variable x(t), the one step prediction of x(t + τ ) is realized, which is the first component of the reconstructed state vector These methods have inspired many data-driven analysis and have been fruitfully reported in the study of multi-disciplines from population study to ecological and biological data even to physiological data and climate data [5,15,20,26,28,30,37,40,42,45,56,65,81,83,89,95,98,101,104].
The research of using local information provided by nearest neighbors to make predictions has been active and attracting even to this day. In order to leverage information from multiple time series observations and overcome the curse of dimensionality, Ye et al proposed a new method in 2016 with name Multiple View Embedding (MVE) [104] which was claimed to be able to make predictions for noisy high-dimensional systems. According to the generalized embedding theorem stated in Theorem 2. L .
For a specific variable y, each reconstruction could be used to make a one-step prediction using nearest neighbors in the way y t+1 ≈ y nn i (t)+1 where nn i (t) is the time index of the nearest neighbor in the ith reconstructed attractor. Though each single prediction may have reduced precision due to the limited data length and noise, the number of predictions could be fairly large. The MVE method firstly ranks all these predictions according to in-sample forecast accuracy, and then the MVE forecast is defined as a simple average of the top k reconstructions in the waŷ where k is empirically set to k = √ m. The base of information leverage by MVE follows from two aspects. Large number of attractors could be reconstructed due to high-dimensional time series data; and each reconstructed attractor filters the dynamical information of the whole system in a different way. Therefore, combining the information from multiple reconstructed attractors could be advantageous in the manner of Eq.(3.4).
x(t+ ) x(t 1 + ) Despite of the rapid development of the nearest neighbors based methods in the recent years, one issue on whether nearest neighbors really reflect the local information near the true target is essential to the success of applying such methods. In Fig. 2(b)(c), two attractors are reconstructed for the same system by different lengths of observed time series. Obviously, when the observed time series is not long enough, the points on the reconstructed attractor are sparse and thus the nearest neighbors actually cannot provide local information of future dynamics. Thus, in order to get the true nearest neighbors, the attractor is necessary to be sufficiently sampled and the time series should be long enough. However, real data sets always have a limited length due to the experiments restriction or cost. Therefore, such classic methods are generally not suitable to be used in dealing with short-term data.
On the other hand, the short-term data are usually accompanied with high-dimensional measurements, i.e., the number of state variables n is much larger than the length of measurements m. Representative examples of such data are high throughput biological data, e.g. microarray data and gene sequencing data, where tens of thousands of gene probes can be measured simultaneously, but only some tens or even fewer successive time points can be measured for expression of each gene [72] due to various restrictions. Collective dynamics of different complex networks [88] such as social networks and biological networks also show high-dimensional but short time course features, where a large number of coupled dynamical systems evolve via various topological structures but only limited time points can be recorded. As a matter of fact, such kind of high-dimensional biological data have emerged in the big data era, from microscopic gene expression data to mesoscopic neural activity data, even to macroscopic ecological and atmosphere data [14,18,35,55,87].
In such high-dimensional nonlinear systems, the variables are not independent but intertwined with each other, and the interactions or correlations between them contain rich information contents, which may reflect the accumulated dynamical features of the target system, as illustrated in Fig. 3. Intuitively, such intertwined relationship, as a source of predictability, has been neglected by the nearest neighbors based methods and thus there is a gap between the classic methods and the large amount of big data these days [18].

Inverse embedding method
Prediction of short-term high-dimensional data was first addressed in [60], in which Ma et al proposed a new SSR based method aiming to transform the spacial information into a temporal domain and thus made it possible to predict short-term high-dimensional data in a nonlinear manner. As illustrated in Fig. 4(a), the classic SSR methods start with a single long-term time series and reconstruct a shorter but higher dimensional time series using delay coordinates. By contrast, the method developed in [60] uses the embedding theory in the inverse way, i.e., starting with short high-dimensional time series and recovering the corresponding longer low-dimensional time series. Theoretically, this new idea can be viewed as an inverse direction of the embedding theory and thus transfers the spacial/interaction information into temporal domain. The process of the inverse embedding method is sketched in Fig. 4(b), and the mapping  Figure 3 Illustrations of information in different types of data sets, from Ref. [60] with permission. (a) A short-term low-dimensional time series has insufficient information to reconstruct the attractor for making predictions. (b) A longterm low-dimensional time series has sufficient information to reconstruct the attractor based on the embedding theory. Many existing methods are based on such a scheme. (c) A short-term high-dimensional time series also contain sufficient information, similar to the long low-dimensional time series, which can be used to reconstruct dynamics or the attractor and thereby make predictions on the time series for some selected variables.
Φ between the original high dimensional attractor A and the reconstructed attractor M using delayed coordinate, as the key to the method, can be formulated in the form as shown in Fig. 5 where the shaded zones stand for the unknown/predicted values.
A simple way to construct the map Φ is solving it as a linear map, which can be obtained in a very efficient way and can be applied to the prediction of many real systems. However, generally the map Φ is an unknown nonlinear vector function and the constraints are limited according to Fig. 5. In order to derive the map in an accurate manner so that the transformation can be carried out, it is firstly expanded into a basis functional space, i.e., the power expansion, the RBF network, or the wavelet expansion. Considering the short-term property, the training problem will be generally categorized into an underdetermined problem where there are more coefficients to be determined than the number of independent equations. Note that in many situations, especially with a proper choice of basis functional space, the coefficients in the expanding series can be sparse, thus the recently developed idea of compressive sensing [9] can be used to solve this problem.
Using the compressive sensing technique, an algorithm to solve the map Φ as well as to make predictions was proposed in [60], which was shown to be particularly effective for short-term high-dimensional

Prediction of New Points
The Corresponding Attractor +τ +2τ +3τ +4τ ( ) +τ +2τ +3τ +4τ +5τ +6τ Time  biological data sets, for instance, microarray data or RNA-seq data, where time points as long as order O(10) were available while thousands of probes were measured and several successful predictions could be made using inverse embedding method. The method has been further considered in the process of big biological data [48] and predicting early-warning signals of critical transitions [52].
To compare the inverse embedding method and nearest neighbors based method, Fig. 6 shows a benchmark test with both methods. A coupled Lorenz system with 15 variables is considered and two criteria for prediction accuracy are depicted, i.e., the correlation ρ and the rooted mean squared error (RMSE) between the test series and the predicted series. In Fig. 6(a)-(b), it is clearly shown that the inverse embedding method can reach quite a good prediction accuracy with training data of less than 100 time points, while the nearest neighbor based method is still in the convergence process even with length of 400. It validates that the inverse embedding method avoids finding nearest neighbors, and thus can deal with short-term data. In Fig. 6(c)-(d), the robustness of the prediction methods against noise deterioration is also tested. It is shown that the nearest neighbors based method is quite robust against noise, which

Original Attractor
Reconstructed Attractor  is mainly attributed to the averaging process of several nearest neighbors. However, compared with the nearest neighbors based method, the inverse embedding method is sensitive to the noise. This is mainly caused by the compressive sensing algorithm adopted in [60], which may enlarge the noise in the process of making parameters sparse. As a matter of fact, though the inverse embedding method transforms the interactions or correlations between variables as spacial information into temporal domain to make predictions, such information actually has redundancy, which is also confirmed in [104]. Therefore, how to dig out such rich redundant information to combat noise could be future work for the inverse embedding method.

Causal Interactions
The interactions between variables within one complex system are essential to understand the mechanism behind the experimental phenomenon, and thus are the foundation of mathematical modelling in a system level. Particularly, data driven interaction inference is in an urgent need in order to model various levels of complex networks. The simplest measure of interaction is the correlation between time series data from two variables x and y. However, correlation only provides the concomitant covariation between x and y and it has been shown that even a perfect correlation between two observed variables is neither necessary nor sufficient to imply a causal relationship [90], especially when it is recognized that nonlinear dynamics are ubiquitous. Thus, compared to non-directional relationship as correlation, the directional interaction inference is more desired. Intuitively, the directional interaction from x to y could be understood as the causal effect from the dynamics of x to the dynamics of y, thus the concept of causality can be introduced. Though there still has been no universally accepted definition of causality, various measures of causality or causal relations have been reported and extensively studied. Among them the Granger Causality (GC), proposed by Clive W.J. Granger, the 2003 Nobel laureate in economy [23,24], is undoubtedly a widely accepted definition and it is also a framework to detect causal relationship from time series data. In the sense of Granger Causality, the definition of causality comes from the idea by Wiener in 1956 [100]: "For two simultaneously measured signals, if we can predict the first signal better by using the past information from the second one than by using the information without it, then we call the second signal causal to the first one". The standard test of GC developed by Granger is based on a linear regression model and later generalized into nonlinear regression [1,62]. Moreover, the original Granger Causality is defined in temporal domain, and the concept could be transformed into frequency domain [7,27,41]. The Granger Causality has been widely accepted in multiple disciplines, especially in economical fields and neural signal analysis, and fruitful applications have been reported in the past decades [4,22]. However, one key assumption in Granger Causality is the separation of effect from one variable to another, i.e., the effect of one variable could be separated and removed from the dynamics of the other, which is generally not true in real situations. Particularly, in complex systems when nonlinearity is dominated, every variable could carry all the system's information, as indicated by Takens' embedding theory, and thus the nonlinear systems are generally non-separative [90]. Moreover, the Granger Causality is defined for stochastic process, and could not be applied to deterministic systems directly.

Nearest neighbors based causality detection
From the viewpoint of information-flow, when x drives y, the dynamics information of x is injected into the dynamics of y, thus if one can discover such footprint, the driving effect from x to y could be detected. Based on this idea, the state space reconstruction framework provides another choice to consider causality in a purely nonlinear style. Following the embedding theorems, each attractor reconstructed by SSR from delayed coordinates of univariate gives one-to-one relation to the original system's attractor, and thus the two attractors reconstructed by x and y respectively are also one-to-one mapped with each other. Explicitly, such one-to-one map can be set up by the time index, i.e., where t is the time index as shown in Fig.7(a), and such map φ is named as a cross map. Moreover, for one particular point y(t 0 ) ∈ M y and the counterpart x(t 0 ) ∈ M x , the local neighborhood of y(t 0 ) can be defined on the attractor as N ǫ (y(t 0 )) = {y(t yi ) y(t yi ) − y(t 0 ) < ǫ} and the mutual neighbors of x(t 0 ) could be further defined as {x(t yi ) y(t yj ) ∈ N ǫ (y(t 0 ))}, as illustrated in Fig.7(a). The mutual neighbors of y(t 0 ) could be analogously defined as show in Fig.7(b). Intuitively, if x and y are causally linked, the mutual neighbors will have some specific properties due to the information injection. As a matter of fact, in the early development of SSR based causality detection, there are two controversial explanations on how to use such property to define the direction of causality.
In [74], Quyen et al claimed that if x drives y but y has no driving force to x, there will be a generalized synchronization relationship between the time series of x(t) and y(t) such that y = F (x) but not vice versa (e.g., no injectivity). Therefore, using the next step of mutual neighbors of y(t n ) ∈ M y , one can make one-step prediction for y(t n+1 ) in the waŷ where N is the number of nearest neighbors in the ǫ−neighborhood. In such a manner, the cross prediction error Pred x→y = ŷ(t n+1 )−y(t n+1 ) can be used to distinguish the causal relation between x and y: if the cross prediction error from x to y is much smaller than the inverse direction, i.e., Pred x→y ≪ Pred y→x , strong evidence of unidirectional causality X → Y is implied.
On the other hand, Schiff et al proposed another explanation in [79]. If x and y are unidirectionally coupled, i.e., information flows only from x to y, the reconstructed attractor M x will present the dynamics in x alone, while the attractor M y reconstructed by y will represent the combined dynamics of both x and y. Thus, x will be predictable by cross prediction using mutual neighbors. Obviously, this is completely contradictory in interpreting the differences in cross prediction errors. . . . . . . x(t-nτ) . . . x(t-nτ) x(t-τ) . . . Figure 7 Illustrations of mutual neighbors, cross map and smoothness, from Ref. [58] with permission. (a) For one point y(t 0 ) ∈ My and its counterpart x(t 0 ) ∈ Mx, one can find the nearest neighbors y(ty 1 ), y(ty 2 ), y(ty 3 ) for y(t 0 ) and define the mutual neighbors x(ty 1 ), x(ty 2 ), x(ty 3 ) for x(t 0 ). The map between the nearest neighbors and mutual neighbors is defined as cross map Φyx. In the case that x causally influences y, the cross map Φyx maps one neighborhood to one neighborhood. (b) In the case y does not causally influence x, the cross map Φxy does not necessarily map one neighborhood to one neighborhood. This issue has not been resolved until Sugihara et al proposed the Convergent Cross Mapping (CCM) method in [90]. In [90], the convergence of cross prediction error was firstly stressed as the length of time series data increases, and it was clearly shown that convergence was essential to distinguish the direction of causal link. If the time series data is not long enough, the variance of the cross prediction errors will be quite large and the confidence area of such errors for both directions would overlap with each other, so that it is not possible to distinguish the two directions, as illustrated in Fig.S5 of [90]. As the length of time series data increases, the reconstructed manifold becomes denser and the highly resolved attractor improves the accuracy of prediction based on nearest neighbors, and thus convergence is observed. With the concept of convergence, both theoretical analysis and benchmark test confirm that the direction indicated in [79] is correct, that is, if x drives y, then one can use the information of y to estimate x, and vice versa.

Cross Map Smoothness method
The convergence concept revealed in [90] actually coincides with the argument proposed in [58,60] as shown in Fig.2(b)-(c), that is, using nearest neighbors on the reconstructed attractor as local approximations essentially relies on the density of points on the reconstructed attractor. Only sufficiently long-term time series can bring faithfully reconstructed attractor with dense points on it and consequently the nearest neighbors can be used as local approximations. As shown in Fig. S5 of [90], to realize convergence, a typical length of order O(10 3 ) was required for a benchmark coupled nonlinear system. Such requirement is usually not practical for typical biological data, e.g., high throughput microarray or RNA-seq data for gene expressions of a biological process. Moreover, on some occasions even though long-term data can be measured, only short (recent) pieces can correctly reflect the causal relation between subsystems due to the nonstationary and fast switching property of the concerned systems. Therefore it is in urgent need of developing new methods to detect causality based on short-term data or a small number of samples.
In [58], Ma et al proposed a new method named Cross Map Smoothness (CMS) which was also based on the SSR technique but was able to infer nonlinear causality from short-term time series data. The base of CMS method lies on the observation that the nearest neighbor based methods actually measure the cross map's local property around each point on the reconstructed attractor, leaving the global information unused. To be specific, if x causally influences y, the mutual neighbors of x(t 0 ) can be used to predict the next value of x, implying that the nearest neighbors of y(t 0 ) are mapped to close states of x(t 0 ) , i.e., the cross map Φ yx : M y → M x is locally smooth around y(t 0 ), as shown in Fig.7(a). On the reverse direction, if y has no influence over x, the image of x(t 0 )'s neighborhood under the cross map Φ xy : M x → M y is not necessarily the neighborhood of y(t 0 ), thus the cross map Φ xy is not necessarily smooth around x(t 0 ), as shown in Fig.7(b). If the cross map is locally smooth around each point on the attractor, it is further globally smooth on the whole attractor, and vice versa, as illustrated in Figs.7(c) and (d). Thus, the global smoothness of Φ yx and Φ xy can be built from local properties. Moreover, as the coupling strength increases, information becomes more distinct between the causally influenced variables. As a result, the reconstructed attractors will contain stronger historical information from the causes. Thus, finding mutual nearest neighbors is equivalent to measuring the smoothness of the cross map, and the relative smoothness indicates the relative strength of causative effectiveness, i.e., the smoothness of Φ yx indicates the strength of causative effectiveness from x to y.
To quantify such smoothness without finding nearest neighbors, an efficient algorithm CMS was designed in [58]. The idea of CMS lies on the fact that any smooth map can be approximated by a neural network N while training a neural network to approximate an unsmooth map will fail with large training errors, as illustrated in Figs  The test on the benchmark system carried out in [58] implied that the CMS method could deal with time series data with length around O(10), which greatly reduced the length requirement of O(10 3 ) by nearest neighbor based methods. It is stressed that the SSR methods, including both the nearest neighbor based methods and the cross map smoothness method, are designed to detect causality with weak to moderate coupling in nonlinear systems. For too weak coupling, the footprints from the driving force are too low to be measured, while for strong coupling, various kinds of synchronization will appear [64], and causality cannot be defined between synchronized dynamics. Actually, the problem of synchronization, particularly transient synchronization, between two time series will cause problems in causality detection from short-term data, which should be avoided in practical applications.

Indirect causality
The causality detection from time series is mainly defined in a pair-wise style, which is the base of network topology inference, but there is still a gap between the existing methods and the network topology reconstruction. Though there are some discussions for the existing methods [58,90], the detection of complex causality patterns remains uncomplete in the literature. As indicated in [61], three types of motifs are generic in network inference: fan-out pattern, fan-in pattern, and cascading pattern, as shown in Fig. 9. Intuitively, in the fan-out pattern, multiple variables are driven by one common source and consequently they are correlated and have the same footprints from the common source, as shown in Fig. 9(a). Both the CCM method [90] as well as CMS method [58] have been shown to work on this pattern, and no significant causality has been detected between the driven nodes. This result confirms the difference between correlation and nonlinear causal relation. Moreover, Hirata et al designed a method specially to deal with common driving source in [32] based on the recurrence plot technique, a technique closely related to SSR technique in nonlinear time series analysis. This was demonstrated to work well in climate and brain network reconstruction [33]. While the fan-in pattern shows no possibility to generate spurious relation, the detected causal strengths are actually weaker as reported in [58]. The third motif is the cascading pattern in which x drives y and y drives z but x has no direct driving force over z, as illustrated in Fig. 9(c). In this pattern, since the information of x has also been injected into the dynamics of z via the intermediate variable y, the spurious indirect causality from x to z could be detected by using the SSR methods pair-wisely. Compared to such detected indirect causal relation, it is more desirable to be able to distinguish the direct causal relationship from the indirect causal relationship in the network reconstruction. To the best knowledge of the authors, there is no SSR methods specially designed to deal with this issue to date, which leaves space for future work in the field.

Causality with delay
In the concept of causality, it is generally accepted that "the cause occurs before the effect" [24,100]. This implies that the influence of one event on another cannot be instantaneous but delayed after some specific time. As a matter of fact, in the real world, time delays are intimately related to causation, which may result from the transmission time, switching speed, memory effect, or other effects. A realworld example is protein expression regulated by transcription factors, where the regulation process goes through a cascade of transcriptions and translations that take from minutes to hours to exert their influences. However, the SSR methods, for instance CCM, actually assumes that the time series are on the same footing through the implicit assumption that the directed influence occurs instantaneously, i.e., the detected causality is from x(t) to y(t). Thus, application of the original SSR methods to detect causality may lead to erroneous results. Recently, there are several efforts to incorporate time delay into the CCM scheme. In [103], Ye et al considered a single time shift in the underlying time series so as to distinguish time-delayed causal interactions by CCM method. However, it is found that the simple incorporation of time shift into CCM framework is problematic. Specifically, the detected time delay by CCM may drift under different embedding dimensions [59], which is definitely confusing. An analogous method was also independently proposed in [80], but it also suffered from the same problem. Quite recently, this problem has been fixed in [59], where a new score Cross Map Evaluation (CME) was designed to substitute the old score in the CCM scheme, and it was reported to be capable of detecting causality with single, multiple, and even distributed time delays. Despite of such progresses in incorporating time delay into SSR based methods, there are still theoretical gaps to be filled. When considering nonlinear systems with time delays, actually one is dealing with infinite-dimensional dynamical systems. However, the classic Takens' embedding theorem is proposed for finite-dimensional systems. Among the few studies on this topic, Robinson has set up an infinite-dimensional Takens theorem in [76], which partially gives theoretical foundation for applying SSR method to time delayed systems and sheds light on this problem.

Conclusion
The world is ruled by physical laws and each factor is evolving, together with other factors, under various governing laws. Thus the world is complex and no factor could be really independent from each other, but all the factors are intertwined by the dynamics of others. The state space reconstruction methods, stemming from the embedding theorem, has brought light to the the analysis of such complex and nonlinear world. Particularly, the SSR methods provide a bridge between the observed time series data and the dynamics of the whole hidden original system, making it possible to reconstruct the original system by simply measuring some observed data. Moreover, all these processes could be done in a model-free style and no equations governing the system are required a priori. The SSR methods could work as a framework and be effective in various ways of reconstructing the system such as determining dimensionality, characteristic exponent, steady states, unstable periodic orbits, etc. In this review, we focus on two data based applications of the SSR methods: prediction and causality inference for nonlinear systems. Despite the long developing history and fruitful results, these two fields have been still active and attracting even to these days, and new methods keep on getting developed in the recent years. Particularly, in this review, we present in detail the forefront methods to deal with short-term time series data, and several perspectives for future work have been provided. Such developments would create an indispensable tool to big data analysis and systems reconstruction. Also, these types of methods can be applied to the problems of inferring dynamical network structures [54,105] and detecting dynamical criticality of nonlinear systems before the abrupt changes of system states [12,52,53].