2.1. Distance matrix

To simplify, we consider time series updated in discrete time, in units of some elementary discretization step, taken as unity without loss of generality. Let us denote X(t₁):t ₁=0,…,N ₁−1 and Y(t₂):t ₂=0,…, N ₂−1 the two time series of which we would like to test their dependence. Note that the lengths N ₁ and N ₂ of the two series can in principle be different as our method generalizes straightforwardly to this case, but for the sake of pedagogy, we restrict ourselves here to the case N ₁=N ₂=N. These time series X(t ₁) and Y(t₂) can be very different in nature with largely differing units and meanings. To make them comparable, we normalize them by their respective standard deviations, so that both normalized time series have comparable typical values. From now on, the two time series X(t ₁) and Y(t₂) denote these normalized time series.

We introduce a distance matrix E _{X , Y} between X to Y with elements defined as The value | X(t ₁)−Y(t₂)| defines the distance between the realization of the first time series at time t ₁ and the realization of the second time series at time t ₂. Other distances can be considered and our method described below applies without modifications for any possible choice of distances. Depending on the nature of the time series, it may be interesting to use other distances, which for instance put more weight on large discrepancies | X(t ₁)−Y(t₂)| such as by using distances of the form |X(t₁)−Y(t₂)| ^q with q > 1. In the following, we do not explore this possibility and only use (1).

When Y(t) is the same time series as X(t), a matrix deduced from (1), by introducing a threshold so that entries of the matrix (1) smaller (respectively larger) than the threshold are set to 0 (respectively 1), has been introduced under the name of ‘recurrence plot’ to analyse complex chaotic time series (Eckmann et al. 1987). In the physical literature, the binary matrix deduced from (1) with the use of a threshold for two different time series is called a cross-recurrence plot. This matrix and several of its statistical properties have been used to characterize the cross-correlation structure between pairs of time series (Marwan and Kurths 2002, Marwan et al. 2002, Quiroga et al. 2002, Strozzia et al. 2002).

Consider the simple example in which Y(t)=X(t−k) with k > 0 is a positive constant. Then, ε (t ₁,t ₂) =0 for t ₂=t ₁+k and is typically non-zero otherwise. The detection of this dependence relationship then amounts in this case to finding the line with zero values which is parallel to the main diagonal of the distance matrix. This line defines the affine mapping t ₂=φ(t ₁) = t ₁+k, corresponding to a constant translation. More generally, we would like to determine a sequence of elements of this distance matrix along which the elements are the smallest, as we describe next.

2.2. Optimal path at ‘zero temperature’

When the relationship between X(t ₁) and Y(t ₂) is more complex than a simple constant lead-lag of the form Y(t)=X(t−k), the determination of the correspondence between the two time series is less obvious. An initial approach would correspond to an association of each entry X(t ₁) of the first time series to the value Y(t ₂) of the second time series which makes the distance (1) minimum over all possible t ₂ for a fixed t ₁. This defines the mapping t ₁→t ₂=φ(t ₁) from the t ₁-variable to the t ₂-variable as Note that this procedure analyses each time t ₁ independently of the others. The problem with this approach is that it produces mappings t₂ =φ (t₁ ) with two undesirable properties: (i) numerous large jumps and (ii) absence of one-to-one matching (φ is no longer a function since the curve can have overhangs and ‘cliffs’) which can also be viewed as a backward (non-causal) time propagation. Property (i) means that, in the presence of noise in two time series, with large probability, there will be quite a few values of t ₁ such that φ(t ₁+1)−φ(t ₁) is large and of the order of the total duration N of the time series. Most of the time, we can expect lags to be a slowly varying function of time and large jumps in the function φ are not reasonable. The second property means that, with large probability, a given t ₁ could be associated with several t ₂, and therefore there will be pairs of times t ₁<t ₁′ such that φ(t ₁) > φ(t ₁′): an occurrence in the future in the first time series is associated with an event in the past in the second time series. This is not excluded as lags between two time series can shift from positive to negative as a function of time, as in our example (12) below. But such occurrences should be relatively rare in real time series which are not dominated by noise. Obviously, these two properties disqualify method (2) as a suitable construction of a time correspondence between the two time series. This reflects the fact that the obtained description of the lag structure between the two time series is erratic, noisy and unreliable.

To address these two problems, we first search for a smooth mapping t ₁→t ₂=φ(t ₁): In the continuous time limit, this amounts to imposing that the mapping φ should be continuous. Then, the correspondence t ₁→t ₂=φ(t ₁) can be interpreted as a reasonable time-lag or time-lead structure of the two time series. For some applications, it may be desirable to impose even further constraints by ensuring the differentiability (and not only the continuity) of the mapping (in the continuous limit). This can be done by a generalization of the global optimization problem (4) defined below by adding a path ‘curvature’ energy term. Here, we do not pursue this idea further. Then, the dependence relationship between the two time series is searched for in the form of a mapping t₂ =φ (t₁ ) between the times t₁ of the first time series and the times t₂ of the second time series such that the two times series are the closest in a certain sense, i.e. X(t ₁) and Y(φ (t₁ )) match best, in the presence of these two constraints.

To implement these ideas, our first proposal is to replace the mapping (2) determined by a local minimization by a mapping obtained by the following global minimization: under constraint (3). Note that, without constraint (3), the solution for the mapping of minimization (4) would recover the mapping obtained from the local minimization (2), as the minimum of the unconstrained sum is equal to the sum of the minima. In contrast, the presence of the continuity constraint changes the problem into a global optimization problem.

This problem actually has a long history and has been extensively studied, in particular in statistical physics (see Halpin-Healy and Zhang (1995) for a review and references therein), under the name of the ‘random directed polymer at zero temperature’. Indeed, the distance matrix E _{X , Y} given by (1) can be interpreted as an energy landscape in the plane (t ₁ , t ₂ ) in which the local distance ε (t ₁,t ₂) is the energy associated with the node (t ₁ , t ₂ ). The continuity constraint means that the mapping defines a path or line or ‘polymer’ of equation (t ₁,t ₂=φ(t ₁)) with a ‘surface tension’ preventing discontinuities. The condition that φ(t ₁) is non-decreasing corresponds to the fact that the polymer should be directed (it does not turn backward and there are no overhangs). The global minimization problem (4) translates into searching for the polymer configuration with minimum energy. In the case where the two time series are random, the distance matrix (and thus energy landscape) is random, and the optimal path is then called a random directed polymer at zero temperature (this last term ‘at zero temperature’ will become clear in section 2.3). Of course, we are interested in non-random time series, or at least in time series with some non-random components: this amounts to having the distance matrix and the energy landscape to have hopefully coherent structures (i.e. non-white noise) that we can detect. Intuitively, the lag–lead structure of the two time series will reveal itself through the organization and structure of the optimal path.

It is important to stress the non-local nature of the optimization problem (4), as the best path from an origin to an end point requires knowledge of the distance matrix (energy landscape) E _{X , Y} both to the left as well as to the right of any point in the plane (t ₁ , t ₂ ). There is a general and powerful method to solve this problem in polynomial time using the transfer matrix method (Derrida et al. 1978, Derrida and Vannimenus 1983). Figure 1 shows a realization of the distance (or energy) landscape E _{X , Y} given by (1) and the corresponding optimal path.

The transfer matrix method can be formulated as follows. Figure 2 shows the (t ₁ , t ₂ ) plane and defines the notations. Note in particular that the optimal path for two identical time series is the main diagonal, so deviations from the diagonal quantify lag or lead times between the two time series. It is thus convenient to introduce a rotated frame (t,x) as shown in figure 2 such that the second coordinate x quantifies the deviation from the main diagonal, hence the lead or lag time between the two time series. In general, the optimal path is expected to wander around, above or below the main diagonal of equation x(t)=0. The correspondence between the initial frame (t ₁ , t ₂ ) and the rotated frame (t,x) is given in detail in the Appendix.

The optimal path (and thus mapping) is constructed such that it can either go horizontally by one step from (t ₁,t ₂) to (t ₁+1, t ₂), vertically by one step from (t ₁,t ₂) to (t ₁,t ₂+1) or along the diagonal from (t ₁,t ₂) to (t ₁+1,t ₂+1). The restriction to these three possibilities embodies the continuity condition (3) and the one-to-one mapping (for vertical segments the one-to-one correspondence is ensured by the convention to map t ₁ to the largest value t ₂ of the segment). A given node (t ₁ , t ₂ ) in a two-dimensional lattice carries the ‘potential energy’ or distance ε (t ₁,t ₂). Let us now denote E(t ₁,t ₂) as the energy (cumulative distance (4)) of the optimal path starting from some origin (t _1,0,t _2,0) and ending at (t ₁ , t ₂ ). The transfer matrix method is based on the following fundamental relation: The key insight captured by this equation is that the minimum energy path that reaches point (t ₁ , t ₂ ) can only come from one of the three points (t ₁−1,t ₂),(t ₁,t ₂−1) and (t ₁−1,t ₂−1) preceding it. Then, the minimum energy path reaching (t ₁ , t ₂ ) is nothing but an extension of the minimum energy path reaching one of these three preceding points, determined from the minimization condition (5). Then, the global optimal path is determined as follows. One needs to consider only the sub-lattice (t _1,0,t _2,0)×(t ₁,t ₂) as the path is directed. The determination of the optimal path now amounts to determining the forenode of each node in the sub-lattice (t _1,0,t _2,0)×(t ₁,t ₂). Without loss of generality, assume that (t _1,0,t _2,0) is the origin (0, 0). Firstly, one performs a left-to-right and bottom-to-up scanning. The forenode of the bottom nodes (τ₁,0) is (τ₁−1,0), where τ₁=1,…,t ₁. Then, one determines the forenodes of the nodes in the second-layer at t ₂ = 1, based on the results of the first (or bottom) layer. This procedure is performed for t ₂  = 2, then for t ₂  = 3, …, and so on.

The global minimization procedure is fully determined once the starting and ending points of the paths are defined. Since the lag–leads between two time series can be anything at any time, we allow the starting point to lie anywhere on the horizontal axis t ₂  = 0 or on the vertical axis t ₁  = 0. Similarly, we allow the ending point to lie anywhere on the horizontal axis t ₂=N−1 or on the vertical axis t ₁=N−1. This allows for the fact that one of the two time series may precede the other one. For each given pair of starting and ending points, we obtain a minimum path (the ‘optimal directed polymer’ with fixed end-points). The minimum energy path over all possible starting and ending points is then the solution of our global optimization problem (4) under constraint (3). This equation of this global optimal path defines the mapping t ₁→t ₂=φ(t ₁) defining the dependence relationship between the two time series.

2.3. Optimal path at finite temperature

While appealing, the optimization program (4) under constraint (3) has an important potential drawback: it assumes that the distance matrix E _{X , Y} between the time series X to Y defined by (1) is made only of useful information. But, in reality, the time series X(t ₁) and Y(t ₂) can be expected to contain a significant amount of noise or more generally irrelevant structures stemming from random realizations. Then, the distance matrix E _{X , Y} contains a possible significant amount of noise, or in other words irrelevant patterns. Therefore, the global optimal path obtained from the procedure of section 2.2 is bound to be delicately sensitive in its conformation to the specific realizations of the noises of the two time series. Other realizations of the noises decorating the two time series would lead to different distance matrices and thus different optimal paths. In the case where the noises dominates, this question amounts to investigating the sensitivity of the optimal path with respect to changes in the distance matrix. This problem has actually be studied extensively in the statistical physics literature (see Halpin-Healy and Zhang (1995) and references therein). It has been shown that small changes in the distance matrix may lead to very large jumps in the optimal path, when the distance matrix is dominated by noise. Clearly, these statistical properties would lead to spurious interpretation of any dependence relationship between the two time series. We thus need a method which is able to distinguish between truly informative structure and spurious patterns due to noise.

In a realistic situation, we can hope for the existence of coherent patterns in addition to noise, so that the optimal path can be ‘trapped’ by these coherent structures in the energy landscape. Nevertheless, the sensitivity to specific realizations of the noise of the two time series may lead to spurious wandering of the optimal path, which does not reflect any genuine lag–lead structure. We thus propose a modification of the previous global optimization problem to address this question and make the determination of the mapping more robust and less sensitive to the existence of noise decorating the two time series. Of course, it is in general very difficult to separate the noise from the genuine signal, in the absence of a parametric model. The advantage of the method that we now propose is that it does not require any a priori knowledge of the underlying dynamics.

The idea of the ‘optimal thermal causal path’ method is the following. Building on the picture of the optimal path as being the conformation of a polymer or of a line minimizing its energy E in a frozen energy landscape determined by the distance matrix, we now propose to allow ‘thermal’ excitations or fluctuations around this path, so that path configurations with slightly larger global energies are allowed with probabilities decreasing with their energy. We specify the probability of a given path configuration with energy ΔE above the absolute minimum energy path by a multivariate logit model or equivalently by a so-called Boltzmann weight proportional to exp[−ΔE/T], where the ‘temperature’ T quantifies how much deviation from the minimum energy is allowed. For T→0, the probability for selecting a path configuration of incremental energy ΔE above the absolute minimum energy path goes to zero, so that we recover the previous optimization problem ‘at zero temperature’. Increasing T allows us to sample more and more paths around the minimum energy path. Increasing T thus allows us to wash out possible idiosyncratic dependencies of the path conformation on the specific realizations of the noises decorating the two time series. Of course, for too large temperatures, the energy landscape or distance matrix becomes irrelevant and one loses all information in the lag–lead relationship between the two time series. There is thus a compromise as usual between not extracting too much from the spurious noise (not too small T) and washing out too much of the relevant signal (too high T). Increasing T allows one to obtain an average ‘optimal thermal path’ over a larger and larger number of path conformations, leading to more robust estimates of the lag–lead structure between the two time series. The optimal thermal path for a given T is determined by a compromise between low energy (associated with paths with high Boltzmann probability weight) and large density (large number of contributing paths of similar energies as larger energies are sampled). This density of paths contributing to the definition of the optimal thermal path can be interpreted as an entropic contribution added to the pure energy contribution of the optimization problem of section 2.2. In a sense, the averaging over the thermally selected path configurations provides an effective way of averaging over the noise realizations of the two time series, without actually having to resample the two times series. This intuition is confirmed by our tests below which show that the signal-over-noise ratio is indeed increased significantly by this ‘thermal’ procedure.

Let us now describe how we implement this idea. It is convenient to use the rotated frame ( t, x ) as defined in figure 2, in which t gives the coordinate along the main diagonal of the (t ₁ , t ₂ ) lattice and x gives the coordinate in the transverse direction from the main diagonal. Of course, the origin (t ₁=0,t ₂=0) corresponds to (x=0,t=0). Note that the constraint that the path is directed allows us to interpret t as an effective time and x as the position of a path at that ‘time’ t. Then, the optimal thermal path trajectory ⟨x(t)⟩ is obtained by the following formula In this expression, G _◃(x,t) is the sum of Boltzmann factors over all paths emanating from (0, 0) and ending at ( x, t ) and G _◃(t) = ∑ _x G _◃(x,t). In statistical physics, G _◃(x,t) is called the partition function constrained to x while G _◃(t) is the total partition function at t. Then, G _◃(x,t)/G _◃(t) is nothing but the probability for a path to be at x at ‘time’ t. Thus, expression (6) indeed defines ⟨ x⟩ as the (thermal) average position at time t. It is standard to call it ‘thermal average’ because G is made up of the Boltzmann factors that weight each path configuration. The intuition is to imagine the polymer/path as fluctuating randomly due to random ‘thermal kicks’ in the quenched random energy landscape. In the limit where the temperature T goes to zero, G _◃(x,t)/G _◃(t) becomes the Dirac function δ[x-x_DP(t)] where x_DP(t) is the position of the global optimal path determined previously in section 2.2. Thus, for T→0, expression (6) leads to ⟨x⟩ = x_DP(t), showing that this thermal procedure generalizes the previous global optimization method. For non-vanishing T, the optimal thermal average ⟨x(t)⟩ given by (6) takes into account the set of neighbouring (in energy) paths, which allows one to average out the noise contribution to the distance matrix. The Appendix gives the recursion relation that allows us to determine G _◃(x,t). This recursion relation uses the same principle and thus has the same structure as expression (5) (Wang et al. 2000).

Similarly to expression (6), the variance of the trajectory of the optimal thermal path reads The variance gives a measure of the uncertainty in the determination of the thermal optimal path and thus an estimate of the error in the lag–lead structure of the two time series as seen from this method.

3.1. Construction of the numerical example

We consider two stationary time series X(t ₁) and Y(t ₂), and construct Y(t ₂) from X(t ₁) as follows: where a is a constant, τ is the time lag and the noise η∼N(0,σ_η) is serially uncorrelated.

The time series X(t ₁) itself is generated from an AR process: where b<1 and the noise ξ∼N(0,σ_ξ) is serially uncorrelated. The factor f=σ_η/σ_ξ quantifies the amount of noise degrading the dependence relationship between X(t ₁) and Y(t ₂). A small f corresponds to a strong dependence relationship. A large f implies that Y(t ₂) is mostly noise and becomes unrelated to X(t ₁) in the limit f→∞. Specifically, and

In our simulations, we take τ = 5, a = 0.8, b = 0.7 and σ_ξ =1 and consider time series of duration N = 100.

For a given f, we obtain the optimal zero-temperature path by using the transfer-matrix method (5) explained in section 2.2 for 19 different starting positions around the origin and similarly 19 different ending positions around the upper-right corner at coordinate (99,99). This corresponds to solving 19×19 transfer matrix optimization problems. The absolute optimal path is then determined as the path which has the smallest energy over all these possible starting and ending points. We also determine the optimal thermal paths ⟨x(t)⟩, for different temperatures, typically from T=1/5 to 10, using relation (A3a) for the partition function and the definition (B1a) for the average transverse path trajectory (given in the Appendix).

Figure 3(a) shows that transverse trajectory x(t) as a function of the coordinate t along the main diagonal for f=1/10 and for temperatures T = 0, 1/5, 1 and 10. This graph corresponds to the case where we restrict our attention to paths with fixed imposed starting (origin) and ending (coordinates (99,99) on the main diagonal) points. This restriction is relaxed as we explain above and is applied below to prevent the boundary effects clearly visible in figure 3(a). Figure 3(b) shows the corresponding standard deviation defined by (7) of the thermal average paths.

The impact of the temperature is nicely illustrated by plotting how the energy of an optimal thermal path depends on its initial starting point x(0)=x ₀ (and ending point taken with the same value x(99)=x(0)). For a given x ₀ and temperature T, we determine the thermal optimal path and then calculate its energy e_T (x ₀) by the formula By construction, the time lag between the two time series is τ = 5 so that we should expect e_T  (x ₀) to be minimum for x ₀=τ =5. Figure 4 plots e_T (x ₀) as a function of the average of the path ⟨x(x₀)⟩ with different starting points x ₀ for different temperatures T respectively equal to 1/50, 1/5, 1/2, 1, 2, 5 and 10 and for f=1/2. One can observe a large quasi-degeneracy for small temperatures, so that it is difficult to identify what is the value of the lag between the two time series. The narrow trough at ⟨x(x ₀)⟩ = 5 for the smallest temperatures, while at the correct value, is not clearly better than negative values of ⟨x(x₀)⟩. In contrast, increasing the temperature produces a well-defined quadratic minimum bottoming at the correct value ⟨x(x ₀)⟩ =τ =5 and removes the degeneracies observed for the smallest temperatures. This numerical experiment illustrates the key idea underlying the introduction of the thermal averaging in section 2.3: too small temperatures lead to optimal paths which are exceedingly sensitive to details of the distance matrix, these details being controlled by the specific irrelevant realizations of the noise η in expression (8). The theoretical underpinning of the transformation from many small competing minima to well-defined large scale minima as the temperature increases, as observed in figure 4, is well understood from studies using renormalization group methods (Bouchaud et al. 1991).

Figure 5 further demonstrates the role of the temperature for different amplitudes of the noise η. It shows the position as a function of T for different relative noise levels f. Recall that ⟨x(t)⟩ is the optimal thermal position of the path for a fixed coordinate t along the main diagonal, as defined in (6). The symbol expresses an additional average of ⟨x⟩ over all the possible values of the coordinate t: in other words, is the average elevation (or translation) of the optimal thermal path above (or below) the diagonal. This average position is an average measure (along the time series) of the lag/lead time between the two time series, assuming that this lag/lead time is the same for all times. In our numerical example, we should obtain close to or equal to τ = 5. Figure 5 shows the dependence of as a function of T for different values of f.

Obviously, with the increase of the signal-to-noise ratio of the realizations which is proportional to 1/f, the accuracy of the determination of τ improves. For a noise level f, approaches the correct value τ = 5 with increasing T. The beneficial impact of the temperature is clearer for more noisy signals (larger f). It is interesting to notice that an ‘optimal range’ of temperature appears for large noise level.

3.2. Test on the detection of jumps or change-of-regime in time lag

We now present synthetic tests of the efficiency of the optimal thermal causal path method to detect multiple changes of regime and compare the results with a standard correlation analysis performed in moving windows of different sizes. Consider the following model In the sense of definition (8), the time series Y is lagging behind X with τ = 0, 10, 5, −5 (this negative lag time corresponds to X(t) lagging behind Y(t)) and 0 in five successive time periods of 50 time steps each. The time series X is assumed to be the first-order AR process (9) and η is a Gaussian white noise. Our results are essentially the same when X is itself a white Gaussian random variable. We use f=1/5 in the simulations presented below.

Figure 6 shows the standard cross-correlation function calculated over the whole time interval 1≤i≤250 of the two time series X and Y given by (12), so as to compare with our method. Without further information, it would be difficult to conclude more than to say that the two time series are rather strongly correlated at zero time lag. It would be farfetched to associate the tiny secondary peaks of the correlation function at τ =±5 and 10 to genuine lags or lead times between the two time series. And since, the correlation function is estimated over the whole time interval, the time localization of possible shifts of lag/leads is impossible.

Before presenting the results of our method, it is instructive to consider a natural extension of the correlation analysis, which consists of estimating the correlation function in a moving window [i+1−D,i] of length D, where i runs from D to 250. We then estimate the lag–lead time τ _D (i) as the value that maximizes the correlation function in each window [i+1−D,i]. We have used D = 10, 20, 50 and 100 to investigate different compromises (D = 10 is reactive but does not give statistically robust estimates while D = 100 gives statistically more robust estimates but is less reactive to abrupt changes of lag). The local lags τ _D (i) thus obtained are shown in figure 7 as a function of the running time i. For D = 10, this method identifies successfully the correct time lags in the first, third, fourth and fifth time periods, while τ _D (i) in the second time period is very noisy and fails to unveil the correct value τ = 10. For D = 20, the correct time lags in the five time periods are identified with large fluctuations at the boundaries between two successive time periods. For D = 50, five successive time lags are detected but with significant delays compared to their actual inception times and with, in addition, high interspersed fluctuations. For D = 100, the delays of the detected inception times of each period reach about 50 time units, that is, comparable to the width of each period, and the method fails completely for this case.

Let us now turn to our optimal thermal causal path method. We determine the average thermal path (transverse trajectory x(i) as a function of the coordinate i along the main diagonal) starting at the origin, for four different temperatures T = 2, 1, 1/2 and 1/5. Figure 8 plots x(i) as a function of i. The time lags in the five time periods are recovered clearly. At the joint points between the successive time periods, there are short transient crossovers from one time lag to the next. Our new method clearly outperforms the above cross-correlation analysis.

The advantage of our new method compared with the moving cross-correlation method for two time series with varying time lags can be further illustrated by a test of predictability. It is convenient to use an example with unidirectional dependence lags (only positive lags) and not with bidirectional jumps as exemplified by (12). We thus consider a case in which X leads Y in general and use the following model At each instant i considered to be the ‘present’, we perform a prediction of Y(i+1) for ‘tomorrow’ at i + 1 as follows. We first estimate the instantaneous lag–lead time τ( i ). The first estimation uses the running-time cross-correlation method which delivers τ(i)=τ _D (i). The second estimation is the average thermal position τ(i)=max[x(i)],0 using the optimal thermal causal path method where the operator [·] takes the integral part of a number. We construct the prediction for Y(i+1) as In this prediction set-up, we assume that we have full knowledge of the model and the challenge is only to calibrate the lag. The standard deviations of the prediction errors are found, respectively, for the cross-correlation method equal to 2.04 for D = 10, 0.41 for D = 20 and 1.00 for D = 50. Using the optimal thermal path, we find a standard deviation of the prediction errors of 0.45 for T = 2, 0.39 for T = 1, 0.33 for T=1/2 and 0.49 for T=1/5. Our optimal causal thermal path method thus outperforms and is much more stable than the classic cross-correlation approach.

4.1. Revisiting the dependence between the US treasury bond yield and the stock market antibubble since August 2000

In a recent paper (Zhou and Sornette 2004), we have found evidence for the following dependence in the time period from October 2000 to September 2003: stock market → Fed Reserve (Federal funds rate) → short-term yields → long-term yields (as well as a direct and instantaneous influence of the stock market on the long-term yields). These conclusions were based on (1) lagged cross-correlation analysis in running windows and (2) the dependence of the parameters of a ‘log-periodic power law’ calibration to the yield time series at different maturities (see Sornette and Johansen (2001), Sornette and Zhou (2002) and Sornette (2003) for a recent exposition of the method and synthesis of the main results on a variety of financial markets).

Let us now revisit this question by using the optimal thermal causal path method. The data consist of the S & P 500 index, the Federal Funds rate (FFR) and ten treasury bond yields spanning three years from 9 September 2000 to 9 September 2003. The optimal thermal paths, x(i)'s, of the distance matrix between the monthly returns of the S & P 500 index with each of the monthly relative variations of the eleven yields are determined for a given temperature T, giving the corresponding lag–lead times τ(i)=x(i)'s as a function of present time i. Figure 9 shows these τ( i )'s for T = 1, where positive values correspond to the yields lagging behind or being caused by the S & P 500 index returns. The same analysis was performed also for T = 10, 5, 2, 1, 1/2 and 1/5, yielding a very consistent picture, confirming indeed that τ is positive for short-term yields and not significantly different from zero for long-term yields, as shown in figure 9. One can also note that the lag τ( i ) seems to have increased with time from September 2000 to a peak in the last quarter of 2003.

We also performed the same analysis with weakly and quarterly data of the returns and yield changes. The results (not shown) confirm the results obtained at the monthly time scale. This analysis seems to confirm the existence of a change of regime in the sign of the lag of the dependence between the S & P 500 index and the Federal Funds rate: it looks as if the Fed (as well as the short term yields) started to be influenced by the stock market after a delay following the crash in 2000, waiting until mid-2001 for the causality to be revealed. The positivity of the time lag suggests that the yields follow (are caused by ?) the stock index. This phenomenon is consistent with the evidence previously presented by Zhou and Sornette (2004) and thus provides further evidence on the dependence flowing from the stock market to the treasury yields. The instantaneous lag–lead functions τ( t ) actually provide much clearer signatures of the dependence structure than our previous analysis: compare for instance with the cross-correlation coefficient shown in figure 10 of Zhou and Sornette (2004). From an economic view point, we interpret this lag relationship as meaning that the FRB is causally influenced by the stock market (at least for the studied period); it should be cautioned however that we have not proved causality stricto sensu, only the existence of a clear lagged relationship between FRB and the stock market. Persuing our tentative interpretation, we take the established relationship as an indication that the stock markets are considered as proxies of the present and are conditioning the future health of the economy, according to the FRB model of the US economy. In a related study, causality tests performed by Lamdin (2004) also confirm that stock market movements precede changes in yield spread between corporate bonds and government bonds. Abdulnasser and Manuchehr (2002) have also found that Granger causality is unidirectionally running from stock prices to effective exchange rates in Sweden.

4.2. Are there any dependence relationships between inflation and gross domestic product (GDP) and inflation and unemployment in the USA?

The relationship between inflation and real economic output quantified by GDP has been discussed many times in the last several decades. Different theories have suggested that the impact of inflation on the real economy activity could be either neutral, negative or positive. Based on Mundell's (1963) story that higher inflation would lower real interest rates, Tobin (1965) argued that higher inflation causes a shift from money to capital investment and raises output per capita. Conversely, Fischer (1974) suggested a negative effect, stating that higher inflation resulted in a shift from money to other assets and reduced the efficiency of transactions in the economy due to higher search costs and lower productivity. In the middle, Sidrauski (1967) proposed a neutral effect where exogenous time preference fixed the long-run real interest rate and capital intensity. These arguments are based on the rather restrictive assumption that the Phillips curve (inverse relationship between inflation and unemployment), taken in addition to be linear, is valid.

To evaluate which model characterizes better real economic systems, numerous empirical efforts have been performed. Fama (1982) applied the money demand theory and the rational expectation quality theory of money to the study of inflation in the USA and observed a negative relation during the post-1953 period. Barro (1995) used data for around 100 countries from 1960 to 1990 to assess the effects of inflation on economic output and found that an increase in average inflation led to a reduction of the growth rate of real per capita GDP, conditioned by the fact that the inflation was high. Fountas et al. (2002) used a bivariate GARCH model of inflation and output growth and found evidence that higher inflation and more inflation uncertainty led to lower output growth in the Japanese economy. Apergis (2004) found that inflation causally affected output growth using univariate GARCH models with a panel set for the G7 countries.

Although cross-country regressions explain that output growth often obtains a negative effect from inflation, Ericsson et al. (2001) argued that these results are not robust and demonstrated that annual time series of inflation and the log-level of output for most G7 countries are co-integrated, thus rejecting the existence of a long-run relation between output growth and inflation. A causality analysis using annual data from 1944 to 1991 in Mexico performed by Shelley and Wallace (2004) showed that it is important to separate the changes in inflation into predictable and unpredictable components whose differences respectively had a significant negative and positive effect on real GDP growth. Huh (2002) and Huh and Lee (2002) utilized a vector autoregression (VAR) model to accommodate the potentially important departure from linearity of the Phillips curve motivated by a strand of theoretical and empirical evidence in the literature suggesting nonlinearity in the output–inflation relationship. The empirical results indicated that their model captured the nonlinear features present in the data in Australia and Canada. This study implies that there might exist a nonlinear causality from inflation to economic output. It is therefore natural to use our novel method to detect possible local nonlinear lagged-dependence relationships.

Our optimal thermal causal path method is applied to the GDP quarterly growth rates paired with the inflation rate updated every quarter on the one hand and with the quarterly changes of the inflation rates on the other hand, for the period from 1947 to 2003 in the USA. The GDP growth rate, the inflation rate and and the inflation rate changes have been normalized by their respective standard deviations. The inflation and inflation changes are calculated from the monthly customer price index (CPI) obtained from the Fed II database (federal reserve bank). Eight different temperatures T=50, 20, 10, 5, 2, 1, 1/2 and 1/5 have been investigated.

Figure 10 shows the data used for the analysis, that is, the normalized inflation rate, its normalized quarterly change and the normalized GDP growth rate from 1947 to 2003.

Figure 11 shows the lag–lead times τ(t)=x(t)'s (units in year) for the pair (inflation, GDP growth) as a function of present time t for T = 2 and for 19 different starting positions (and their ending counterparts) in the (t ₁ , t ₂ ) plane, where positive values of τ(t)=x(t) correspond to the GDP lagging behind or being caused by inflation. This figure is representative of the information at all the investigated temperatures. Overall, we find that τ is negative in the range −2 years ≤τ≤0 year, indicating that it is more the GPD which leads inflation than the reverse. However, this broad-brush conclusion must be toned down somewhat at a finer time resolution as two time periods can be identified in figure 11.

1. From 1947 (and possibly earlier) to the early 1980s, one can observe two clusters, one with negative −2 years ≤τ = x(t)≤0 years implying that the GDP has a positive lagged effect on future inflation, and another with positive 0 years ≤τ = x(t)≤4 years implying that inflation has a lagged effect on GDP with a longer lag.

2. From the mid-1980s to the present, there is no doubt that it is GDP which has had the dominating impact on future inflation lagging behind by about 1–2 years.

In summary, our analysis suggests that the interaction between GDP and inflation is more subtle than previously discussed. Perhaps past controversies in which one causes the other may be due to the fact that, to a certain degree, these different time series have entangled dynamics, with different time lags. Any measure of a dependence relationship allowing for only one lag is bound to miss such subtle interplay. It is interesting to find that GDP impacts on future inflation with a relatively small delay of about one year while inflation has in the past influenced future GDP with a longer delay of several years.

Figure 12 shows the lag–lead times τ(t)=x(t)'s (units in year) for the pair (inflation change, GDP) as a function of present time t for T = 2 and for 19 different starting positions (and their ending counterparts) in the (t ₁ , t ₂ ) plane, where positive values of τ(t)=x(t) correspond to the GDP lagging behind or being caused by inflation change. Due to the statistical fluctuations, we cannot conclude on the existence of a significant dependence relationship between inflation change and GDP, except in the decade of the 1980s for which there is a strong causal effect of a change of inflation on GDP. The beginning of this decade was characterized by a strong decrease of the inflation rate from a two-digit value in 1980, following a vigorous monetary policy implemented under the Fed's chairman Paul Volker. The end of the 1970s and the better half of the 1980s were characterized by an almost stagnant GDP. In the mid-1980s, the GDP started to grow again at a strong pace. It is probably this lag between the significant reduction of inflation in the first half of the 1980s and the rise of the GDP growth that we detect here. Our analysis may help in improving our understanding in the intricate relationship between different economic variables and their impact on growth and on stability and in addressing the difficult problem of model errors, which Cogley and Sargent (2005) have argued to be the cause for the lack of significant action from the Fed in the 1970s.

Figure 13 shows the lag–lead times τ(t)=x(t)'s (units in year) for the pair (inflation, unemployment rate) as a function of present time t for T = 2 and for 19 different starting positions (and their ending counterparts) in the (t ₁ , t ₂ ) plane, where positive values of τ(t)=x(t) correspond to the unemployment rate lagging behind or being caused by inflation. We use quarterly data from 1948 to 2004 obtained from the Fed II database (federal reserve bank). This figure is representative of the information at all the investigated temperatures.

1. From 1947 (and possibly earlier) to 1970, one can observe large fluctuations with two clusters, suggesting a complex dependence relationship between the two time series, similarly to the situation discussed above for the (inflation, GDP) pair.

2. From 1970 to the present, there is no doubt that inflation has predated unemployment in the sense of the optimal thermal causal path method. It is also noteworthy that the lag between unemployment and inflation has disappeared in recent years. From a visual examination of figure 13, we surmise that what is detected is probably related to the systematic lags between inflation and employment in the four large peak pairs: (1970 for inflation; 1972 for employment), (1975 for inflation; 1976 for unemployment), (1980 for inflation; 1983 for unemployment) and (1991 for inflation; 1993 for unemployment).

One standard explanation for a lagged impact of inflation on unemployment is through real wage: if inflation goes faster than the adjustment of salaries, this implies that real wages are decreasing, which favours employment according to standard economic theory, thus decreasing unemployment. Here, we find that surges of inflation precede the increase and not the decrease of unemployment. Rather than an inverse relationship between synchronous inflation and unemployment (Phillips curve), it seems that a better description of the data is a direct lagged relationship, at least in the last thirty years. The combination of increased inflation and unemployment has been known as ‘stagflation’ and caused policymakers to abandon the notion of an exploitable Phillips curve trade-off (see for instance Lansing (2000)). Our analysis suggests a more complex multivariate description which requires taking into account inflation, inflation change, GDP, unemployment and their expectations by the agents, coupled all together through a rather complex network of lagged relationships. This retrieves within a non-parametric method the conclusion obtained by macroeconomists who nowadays use various classes of endogenous growth models with different usage of physical and human capital and with different exchange technologies (Gillman and Kejak 2005). But, the fact that we are retrieving such conclusions within a non-parametric method is interesting, we think, since the issue is still actively investigated empirically by modern authors (see for instance Barr (1995), Fountas et al. (2002) and Apergis (2004)). Our contribution is to add to this evidence by using a new original method.