Common dynamical properties of business cycle fluctuations are studied in a sample of more than 100 countries that represent economic regions from all around the world. We apply the methodology of multivariate singular spectrum analysis (M-SSA) to identify oscillatory modes and to detect whether these modes are shared by clusters of phase- and frequency-locked oscillators. An extension of the M-SSA approach is introduced to help analyze structural changes in the cluster configuration of synchronization. With this novel technique, we are able to identify a common mode of business cycle activity across our sample, and thus point to the existence of a world business cycle. Superimposed on this mode, we further identify several major events that have markedly influenced the landscape of world economic activity in the postwar era.

Despite a long tradition of systematically analyzing cyclic behavior in economic data,1–3 the nature of aggregate fluctuations is still one of the most controversial topics in macroeconomics.4,5 Although the emergence of business cycle synchronization across countries has been widely acknowledged6–9—especially in view of the ongoing globalization of economic activity—there is still no agreement on basic issues like the quantification of comovements. Over the years, economic developments, changes in resource availability, and changes in political systems can result in more or less drastic changes in the patterns of economic activity, and therefore change many aspects of synchronization. In the present work, we apply multivariate singular spectrum analysis (M-SSA) to identify common spectral properties in a sample of macroeconomic time series from over 100 countries that represent economic regions from all around the world. An M-SSA extension introduced herein helps us explore the cluster configuration of synchronization in this sample, as well as identify several major events that have markedly influenced world economic activity in the postwar era. A common mode of business cycle activity is found and it points to the existence of a world business cycle.

Over the last quarter-century, multivariate singular spectrum analysis (M-SSA) has proven its efficiency and reliability in the spatio-temporal analysis of large datasets in several fields of the geosciences and of other disciplines.10 M-SSA provides insight into the unknown or partially known dynamics of the underlying phenomena by decomposing the delay-coordinate phase space of a given multivariate time series into a set of data-adaptive orthogonal components. These components can be classified essentially into long-term trends, oscillatory patterns, and residual noise, and they allow one to reconstruct a robust “skeleton” of the dynamical system's structure.10–12 While this skeleton does not yield, in general, the dimension of the system's attractor,11,13,14 it can greatly help phase synchronization analysis and provides considerable insight into the mechanisms of rhythm adjustment.15 

Phase synchronization refers, in general, to an adjustment of rhythms of coupled oscillators that is reflected in a locking of both their frequencies and phases.16–18 In the presence of spiral behavior, the phase is typically defined as an angle of rotation with respect to an origin in phase space.19–21 For low-dimensional systems, this phase definition can be based on visual inspection; for high-dimensional systems and in the presence of noise, it may become more difficult to formulate a phase definition that is both useful and robust. In practice, approaches relying on the definition of an angle usually depend on a priori knowledge about the analyzed system.

The data-adaptive M-SSA approach, on the other hand, is able to automatically identify oscillatory modes and detect cluster synchronization in large systems of coupled oscillators, while no detailed knowledge of individual subsystems nor a suitable phase definition for each of them is required.15,22–24

Despite the algorithm's proven skill in identifying coupled oscillatory modes, the standard M-SSA approach assumes stationarity in the cluster configuration of the underlying system, a stationarity that does not appear to be very realistic in the context of economic activity. Over the years, economic developments, changes in resource availability, or changes in political systems can result in more or less drastic changes in the patterns of economic activity. These changes in the way a market or an economy operate can result in changes of the synchronized portion of activity and in the cluster configuration of synchronized subsystems.

In the present paper, we introduce an extension of the M-SSA approach to the analysis of structural changes in the patterns of synchronization. The two-stage approach implemented herein combines the M-SSA analysis with a subsequent complex empirical orthogonal function (EOF) analysis. Our approach helps systematically decompose the possibly intricate spatio-temporal structure of the identified oscillatory modes into a set of distinct clusters in time.

The interested reader might note the strongly interdisciplinary character of the paper—which combines (i) an advanced method of time series analysis, namely M-SSA; (ii) applied dynamical systems theory, in the shape of synchronization analysis; and (iii) a state-of-the-art aspect of macroeconomic theory, namely synchronization of business cycles—and the challenges this interdisciplinarity may pose upon a first reading. Every effort has been made to present the methodology, as well as the results, as simply and self-consistently as possible. Remaining difficulties are entirely the authors' fault and may be overcome upon a second reading.

In Sec. II, we give a brief overview of the economic background and the synchronization of economic activity. In Sec. III, we review the main properties of the standard M-SSA methodology and introduce a two-stage approach that combines a complementary M-SSA algorithm with a subsequent complex EOF analysis. In Sec. IV, we show how this two-stage approach can help detect changes in the pattern of synchronization in a chain of coupled oscillators, while in Sec. V we study the synchronization of economic activity in the large sample of macroeconomic indicators selected herein. A summary and concluding remarks make up Sec. VI, and two  Appendixes provide further technical details.

Macroeconomic time series are dominated, over many decades, by a long-term upward trend; they also exhibit smaller but still significant shorter-term fluctuations that are often associated with business cycles.25 The causes and characteristics of these cycles have been extensively studied in modern economic theory, while the debate on their nature—such as the endogenous vs. exogenous nature of business cycles and their propagation mechanisms—is still very much alive.5,25–27

A number of approaches have been proposed to separate the shorter-term fluctuations from the long-term trend.28–30 Since there is, however, no fundamental theory—and hence no generally accepted definition—of the trend, the resulting residuals have to be analyzed very critically.

On the other hand, it is widely acknowledged that business cycles are multi-national phenomena, showing common characteristics across countries.31,32 Still, there is no agreement on basic issues like the quantification of comovements, the existence of supranational cycles6–9—for instance at the European Union or G7 level—and the determinants of economic synchronization. For this reason, many theoretical and empirical studies disagree in their results, due to different datasets as well as different methodologies. Hence, the lack of agreement on the extent and nature of macroeconomic synchronization persists.

Univariate singular spectrum analysis has already been applied to study business cycles of individual macroeconomic indicators in a single country,33,34 while a first application of M-SSA across indicators has been shown to give deeper insights into the dynamical properties of US business cycles.35 Only quite recently, though, has business cycle synchronization been studied in a joint M-SSA analysis of a small set of three European economies and the US.36 In the latter study, common cyclical characteristics have been identified, both across countries and across indicators.

The present paper extends the study of business cycles and of the determinants of economic synchronization to the global scale. In an M-SSA analysis of more than 100 countries, covering economic regions from all around the world, we investigate here, in Sec. V, common dynamical properties of business cycle fluctuations and their spectral properties across countries, regions, and the world. The comprehensive dataset is analyzed in a unified M-SSA analysis, which provides us with a consistent separation into a permanent trend component and transient fluctuations that are orthogonal to it.

In this subsection, we summarize the main aspects of the standard textbook-version of the M-SSA algorithm.10,37 The algorithm involves four main steps, (1) embedding, (2) decomposition, (3) rotation, and (4) reconstruction; these steps are outlined in the following.

1. Embedding

We rely here on the trajectory-matrix approach,38,39 which starts by embedding each channel d of the multivariate time series x={xd(n):d=1,,D;n=1,,N}—with D channels of length N—into an M-dimensional trajectory matrix by using lagged copies

(1)

Each matrix Xd has M columns of reduced length N=NM+1, and we form the augmented trajectory matrix by concatenating all the D channels

(2)

The full trajectory matrix X has, therefore, DM columns of reduced length N.

2. Decomposition

The M-SSA algorithm proceeds by performing a singular value decomposition (SVD)

(3)

where the prime E denotes the transpose of E. The normalization factor η equals max{N,DM}. The matrix Σ has κ non-vanishing diagonal elements, which are its singular values {sk:k=1,,κ}, where κ=min{N,DM}. The two matrices of singular vectors P and E have both rank κ and they provide a set of empirical orthogonal functions (EOFs). A schematic diagram of the SVD decomposition in Eq. (3) is provided in Fig. 1(a).

FIG. 1.

Schematic diagram of the time series analysis methodology. (a) The standard M-SSA algorithm in Eq. (3); (b) the complementary M-SSA algorithm in Eq. (8), followed by (c) the subsequent complex EOF analysis in Eq. (13). The diagram illustrates the simple case of two variables, D = 2.

FIG. 1.

Schematic diagram of the time series analysis methodology. (a) The standard M-SSA algorithm in Eq. (3); (b) the complementary M-SSA algorithm in Eq. (8), followed by (c) the subsequent complex EOF analysis in Eq. (13). The diagram illustrates the simple case of two variables, D = 2.

Close modal

It becomes clear that the matrix E of right-singular vectors is composed of D consecutive segments Ed of size M×κ

(4)

each of which is associated with a channel Xd in X. These vectors are referred to as space-time EOFs (ST-EOFs) and represent a set of multivariate data-adaptive filters of length M.

The window length M of the filter is equivalent to the embedding dimension in Eq. (2) and it is typically chosen to cover more than one oscillation period. This way, oscillatory behavior is captured via the so-called oscillatory pairs of two ST-EOFs in phase quadrature.11,40 These pairs are the analogs of sine-and-cosine pairs in Fourier analysis, and M here needs to be no less than a small multiple of its period for any such pair to be reliably identified by the M-SSA analysis.

The matrix P of left-singular vectors has size N×κ. These vectors are referred to as temporal EOFs (T-EOFs) and reflect the corresponding temporal behavior of an oscillation, as captured through the filter lens of the ST-EOFs. This filtering step becomes more apparent via the equivalent projection of the trajectory matrix onto the ST-EOFs

(5)

from which we obtain the principal components (PCs); i.e., the PCs are proportional to the T-EOFs in matrix P, scaled by the singular values Σ.

3. Rotation

To better separate distinct oscillations for large D, we rely here on a modified varimax rotation of the ST-EOFs;15,41 see also  Appendix A for further algorithmic details.

4. Reconstruction

Dynamical behavior in X that is associated with a subset K{1,,κ} of EOFs can be reconstructed from Eq. (3) by

(6)

with K a diagonal matrix of size κ×κ, in which the elements Kkk equal unity if kK, and Kkk = 0 otherwise. Upon averaging along the skew diagonals of RK, i.e., over elements of the same time instance in Eq. (1), we finally obtain the reconstructed components (RCs).

Remark. We have chosen here to obtain the two matrices P and E from the SVD of X in Eq. (3). Alternatively, one could also start with the covariance matrix C=η1XX to first obtain E from the eigendecomposition of C, C=EΛE, and then P from a projection P=η1/2XEΣ1, according to Eq. (5). However, the latter projection requires computing the inverse of Σ, which becomes problematic in the presence of small eigenvalues in Λ=Σ2. Moreover, the calculation through straightforward eigendecomposition provides a more complete picture in the case of a rank-deficient covariance matrix, i.e., DM>N. In the latter case, one needs to calculate a reduced covariance matrix, C̃=η1XX, and its eigendecomposition, C̃=PΛP.23,42 Irrespective of these issues, there exist efficient and robust algorithms to compute the nonvanishing eigenelements directly from a reduced version of the SVD in Eq. (3).43 

The projection in Eq. (5) is optimal in the sense that a maximum amount of variance in X is captured by a minimal number of PCs. Optimality, though, is only given in a total least-mean-squares sense, and the PCs describe overall temporal behavior that is most common to all input time series.

In the joint M-SSA analysis of extensive atmospheric and oceanic datasets,44 for example, coupled ocean-atmosphere modes were most prominent, while separate M-SSA analyses of the oceanic and atmospheric datasets unveiled small but significant mismatches in the temporal behavior of the two fluid media. It was shown that a gradual decline of interannual variability in the atmospheric forcing was not accompanied by a similar decline of interannual variability in the ocean response; i.e., in some regions of the North Atlantic basin, the hypothesis of an ocean response to atmospheric forcing could be rejected with high confidence. Furthermore, separate M-SSA analyses over distinct time intervals showed that the spatio-temporal structure of interannual variability has changed over time.

Analyzing, however, extensive datasets by separate M-SSA analyses on subsets of the space-and-time domain can become quite cumbersome, and a more systematic approach is clearly preferable. In Sec. III C below, we discuss, therefore, a complementary version of the M-SSA algorithm that is flexible enough to reconstruct oscillatory behavior that undergoes structural changes in the space–time domain. The M-SSA oscillatory modes will then be analyzed in greater detail in Sec. III E, using a complex EOF analysis to systematically identify structural changes.

The projection in Eq. (5) averages over all input channels, while the T-EOFs of length N provide a comprehensive least-mean-squares picture of the temporal behavior. This property of M-SSA analysis arises from the decomposition of X in Eq. (3) into an outer product of a single set of T-EOFs in P and the channels' individual contributions in E.

Hence, in a typical situation, the window length M has to be chosen to cover at least once the longest oscillation period of interest. In the original M-SSA formulation, though, it was desirable to use M<N/3 or smaller in order to be able to identify temporal changes in the amplitude of the oscillation.45 While the complementary M-SSA algorithm still involves the same four steps 1–4 of Sec. III A, a simple restructuring of the trajectory matrix in step 1 helps circumvent this limitation; cf. also Fig. 1(b), in which a schematic diagram illustrates the main properties of the algorithm.

1. Embedding

The complementary M-SSA algorithm begins by concatenating the channel-wise transpose of Xd

(7)

into a new trajectory matrix Y of size M×DN.

2. Decomposition

This version of the M-SSA algorithm then continues, as before, with the SVD of the restructured trajectory matrix Y,

(8)

from which we obtain a new set of EOFs. In the following, we assume that M<N/2, and therefore M<N, so that U and V have both rank M, and the normalization factor η equals DN.

In the single-channel case with D = 1, we have X=Y, and the two SVDs in Eqs. (3) and (8) yield the same results, except that the two sets of EOFs interchange their respective roles; see also Fig. A1 in Ref. 10.

In the multi-channel case with D2, however, this restructuring in Y has far-reaching effects on the decomposition of the spatio-temporal variability. The new matrix of left-singular vectors U has now size M × M. To understand its role in the decomposition, remember that U is also the matrix of eigenvectors of the covariance matrix C=η1YY and, following Eq. (7), likewise of C=η1d=1DXdXd. Thus, U now describes spectral properties of the covariance structure that are common to all channels. In contrast to P, these eigenvectors U have now length M and take over the role of data-adaptive filters; i.e., the analog of sine and cosine functions.

Note that a similar idea has already been discussed in the context of Monte Carlo SSA hypothesis testing,23,46 in which the EOF test basis is determined from the eigendecomposition of the covariance matrix associated with the null hypothesis; to wit, the latter is estimated from the average over the entire set of all surrogate covariance matrices. Here, we likewise average over a set of covariance matrices in C, and determine a new EOF basis U that characterizes the common spectral properties of all input channels.

Originally introduced for single-channel SSA analysis,46 the Monte Carlo test for SSA turns out to be potentially rather misleading in the multi-channel setting of M-SSA.23 The EOF basis U describes only common spectral properties. Hence, oscillations of similar periods but with distinct spatio-temporal features, for example, could wind up as being subsumed in the same oscillatory pair of U. This would be the case, for instance, if the test were to associate the same significance level with both oscillations, although their corresponding eigenvalues might be quite different.23 

In the present paper, however, it is exactly this property of the EOF basis U that gives sufficient flexibility to the complementary M-SSA algorithm in associating oscillations of a similar period, but different spatio-temporal patterns. The corresponding spatio-temporal behavior of the oscillation is captured by the right-singular vectors V.

The matrix V has size DN×M and it is composed of D consecutive segments Vd of size N×M,

(9)

each of which is associated with a channel Xd in Y. In contrast to the single set of T-EOFs in matrix P that results from the SVD of Eq. (3), we have now D separate sets of T-EOFs in Eq. (8); see again Figs. 1(a) and 1(b).

3. Rotation

The separability of distinct oscillations is further improved by a modified varimax rotation of V, as in the standard M-SSA algorithm; see  Appendix A for further algorithmic details.

4. Reconstruction

Dynamical behavior in Y that is associated with a subset of EOFs can be reconstructed, in a manner analogous to Eq. (6), from

(10)

with R of size M×DN, while averaging along the skew diagonals of R—i.e., over elements that occur at the same time in Eq. (7)—yields the RCs of x.

Remark. Table I compares the properties of the standard M-SSA algorithm in Eq. (3) with those of the complementary M-SSA algorithm in Eq. (8). As already illustrated in Figs. 1(a) and 1(b), the main differences in the decomposition of the trajectory matrix between the two versions of the M-SSA algorithm can be summarized as follows:

TABLE I.

Properties of the standard M-SSA algorithm in Eq. (3) vs. the complementary M-SSA algorithm in Eq. (8).

Size ofStandard M-SSAComplementary M-SSA
Trajectory matrix X: N×DM Y: M×DN 
Temporal EOF P: N×κ V: DN×κ 
Filter EOF E: DM×κ U: M×κ 
Singular values Σ: κ×κ S: κ×κ 
Rank κ = min{N,DM} κ = M 
Size ofStandard M-SSAComplementary M-SSA
Trajectory matrix X: N×DM Y: M×DN 
Temporal EOF P: N×κ V: DN×κ 
Filter EOF E: DM×κ U: M×κ 
Singular values Σ: κ×κ S: κ×κ 
Rank κ = min{N,DM} κ = M 
  1. The standard M-SSA algorithm captures temporal behavior in the univariate T-EOFs of matrix P, which are common to all input channels, while the complementary M-SSA algorithm captures channel-wise temporal behavior in the multivariate T-EOFs of matrix V.

  2. The standard M-SSA algorithm projects the multi-channel time series onto a multivariate spectral filter given by the ST-EOFs of matrix E, while the complementary M-SSA algorithm projects all time series onto the univariate filter of matrix U.

To better understand the differences between the two versions of the M-SSA algorithm, we consider next the simple example of spatio-temporal harmonic oscillations. While the oscillation period is fixed, amplitudes and phases can vary in space and time.

Figure 2(a) shows the observed input signal, while the noise-free reference signal is shown in Fig. 2(b). The reference signal is composed of three different signals with varying spatio-temporal behavior, which are shown in Figs. 2(c)–2(e). The superposition of the three signals is meant to reflect structural changes in time: it is these changes we try to identify in the noise-contaminated signal in Fig. 2(a). The superimposed red noise is generated by AR(1) processes with parameter γ=0.9, driven by spatially correlated white noise processes whose covariance matrix is given by Wij=0.8|ij|.

FIG. 2.

Spatio-temporal oscillations in the presence of strong correlated noise. (a) Observed pattern; and (b) noise-free reference pattern with D = 50 channels of length N = 300. The reference pattern is composed of three different spatio-temporal oscillations in panels (c)–(e), with the oscillation period T = 20 in all three patterns. Superimposed on the reference signal is spatio-temporally correlated noise, with the signal-to-noise ratio set to 1.5.

FIG. 2.

Spatio-temporal oscillations in the presence of strong correlated noise. (a) Observed pattern; and (b) noise-free reference pattern with D = 50 channels of length N = 300. The reference pattern is composed of three different spatio-temporal oscillations in panels (c)–(e), with the oscillation period T = 20 in all three patterns. Superimposed on the reference signal is spatio-temporally correlated noise, with the signal-to-noise ratio set to 1.5.

Close modal

To cover at least one oscillation period, we set the window length to M = 30. In both M-SSA analyses, the standard and the complementary one, the leading pair of EOFs 1–2 captures oscillatory behavior of period T = 20, with the corresponding spatio-temporal reconstructions shown in Figs. 3(a) and 3(b), respectively. It turns out, though, that only the complementary M-SSA algorithm in Fig. 3(b) is able to capture the time-varying spatio-temporal structure in Fig. 2(b). The reconstruction of the standard M-SSA algorithm in Fig. 3(a), on the other hand, gives only a much fuzzier picture of the spatio-temporal structure of the oscillation, while many structural details are completely missed.

FIG. 3.

Reconstruction of spatio-temporal oscillations with the leading oscillatory pair RCs 1–2 of (a) the standard M-SSA algorithm and (b) the complementary M-SSA algorithm. The window length M = 30 and the leading 10 EOFs varimax rotated. Panels (c)–(e) show partial RCs from a reconstruction with the three leading EOFs of a subsequent complex EOF analysis of M-SSA EOFs 1–2 of panel (b). The complex EOF analysis in T-mode, with the 10 leading complex T-EOFs varimax rotated.

FIG. 3.

Reconstruction of spatio-temporal oscillations with the leading oscillatory pair RCs 1–2 of (a) the standard M-SSA algorithm and (b) the complementary M-SSA algorithm. The window length M = 30 and the leading 10 EOFs varimax rotated. Panels (c)–(e) show partial RCs from a reconstruction with the three leading EOFs of a subsequent complex EOF analysis of M-SSA EOFs 1–2 of panel (b). The complex EOF analysis in T-mode, with the 10 leading complex T-EOFs varimax rotated.

Close modal

Only the central part 20d40 of the signal is correctly reconstructed by the standard M-SSA in Fig. 3(a), while two further oscillatory pairs are necessary to reconstruct the two adjacent parts 1d20 and 40d50, respectively (not shown). This decomposition into spatially distinct patterns, though, is characteristic of the varimax rotation of the ST-EOFs in E, cf. Eqs. (3) and (4) and  Appendix A; it is less helpful, however, in identifying structural changes in the time domain.

In the complementary M-SSA algorithm, all input channels are projected onto the same set of univariate filters in U. Since the oscillation period is identical in the three oscillatory patterns of Figs. 2(c)–2(e), a single oscillatory pair in U suffices in Fig. 3(b) to reconstruct the complex spatio-temporal structure in Fig. 2(b). Despite the algorithm's excellence in identifying the complex spatio-temporal structure of the oscillation in the presence of strong correlated noise, it still lacks more specific information about the existence of structurally distinct patterns. In the present example, visual inspection can already provide useful information in this direction, but a more systematic approach is necessary for high-dimensional systems, as presented in Sec. III E.

For a more systematic decomposition of the spatio-temporal structure of the oscillation, recall that the complementary M-SSA algorithm in Eq. (8) captures temporal behavior in the multivariate T-EOFs of matrix V, cf. Eq. (9) and Fig. 1(b). Suppose we have identified oscillatory behavior in a pair of two T-EOFs, vp and vq, namely the p-th and q-th columns of V, respectively. The subsequent complex EOF analysis involves four main steps 1–4, similar to those in Secs. III A and III C; these steps are described in the following.

1. Embedding

Upon varimax rotation of V, cf.  Appendix A, the two T-EOFs are in phase quadrature and will therefore provide an expansion of the oscillation in the complex plane

(11)

According to Eq. (9), the complex-valued column vector z of length DN is likewise composed of D consecutive segments zd of length N, each of which is associated with a channel Vd in V.

Concatenating the channel-wise segments zd,

(12)

will finally provide the complex input signal Z of size N×D; see Fig. 1(c) for an illustration.

2. Decomposition

The complex EOF analysis then continues with the SVD

(13)

This SVD provides a factorization of Z into two unitary matrices, Uc and Vc, and a diagonal matrix Sc of real numbers, all three of which have rank κ=min{N,D}.47,48 The left-singular vectors Uc describe the temporal behavior of the oscillations and we refer to them as complex T-EOFs. The right-singular vectors Vc provide spatial information and describe the corresponding phase and amplitude relations for each of the input channels, and we refer to them as complex S-EOFs.

3. Rotation

Complex EOF analysis is known to be useful in identifying propagating and standing waves.49 By analogy with the varimax rotation of ST-EOFs in M-SSA, cf. Sec. III D, we seek to simplify the interpretation of the complex EOFs by a varimax rotation. In contrast to the M-SSA rotation, which attempts to simplify the spatial structure of the oscillations, there are two possibilities for achieving simpler oscillatory structures by varimax rotation in complex EOF analysis.50 

In the so-called S-mode,51 a simple structure rotation is applied to the complex S-EOFs and will result in spatially distinct clusters, the same as in the ST-EOF rotation for the M-SSA algorithm.

In the alternative T-mode, a rotation is applied to the complex T-EOFs and will result in distinct clusters in time. In the geosciences, the S-mode is more frequently used, although the T-mode is not uncommon either.52,53 Later on, we will see that these two modes provide complementary information about the system under study.

4. Reconstruction

To reconstruct dynamical behavior in the original input time series x from the two-stage decomposition of Eqs. (8) and (13), we first reconstruct the dynamical behavior Zk that is associated with the k-th complex EOF

(14)

Here, K is a diagonal matrix of rank κ, in which the k-th diagonal element equals unity and the other diagonal elements are set to zero. Following Eq. (11), the real and imaginary parts of Zk are associated with T-EOFs vp and vq, respectively; see once more Figs. 1(b) and 1(c), in which the association is illustrated by two dashed arrows. That is, by replacing the T-EOFs vp and vq in Eq. (10) by their counterparts in Zk, we obtain, again upon diagonal averaging, the corresponding partial RCs.

Note that k=1κZk yields a complete reconstruction of Z, i.e., of T-EOFs vp and vq, and the sum of all partial RCs from the subsequent complex EOF analysis likewise yields the original RCs from the M-SSA analysis.

In the above example of harmonic oscillations, we have already seen that the time-varying spatio-temporal structure in Fig. 2(b) is well captured by RCs 1–2 of the complementary M-SSA algorithm in Fig. 3(b). Figures 3(c)–3(e) show the resulting partial RCs as derived from a reconstruction with the three leading EOFs of a subsequent complex EOF analysis in the T-mode. In this analysis, the complex T-EOFs are varimax rotated, while the partial RCs give the desired reconstruction of the three reference patterns in Figs. 2(c)–2(e). A complex EOF analysis in the alternative S-mode (not shown) yields only spatially distinct patterns, similar to the standard M-SSA algorithm in Fig. 3(a).

In this paper, we discuss a two-stage approach that combines the M-SSA analysis with a subsequent complex EOF analysis. Namely, in the first step, the M-SSA algorithm identifies a spatio-temporal oscillation and expands it into the complex plane, cf. Fig. 1(b), while the complex EOF analysis then decomposes this spatio-temporal oscillation into structurally distinct patterns, cf. Fig. 1(c).

Complex EOF analysis is often used to recognize wave patterns, and it is commonly based on a complex expansion of the input signal via a Hilbert transform.47–49 Unless some bandpass filtering is performed beforehand, no particular time scale is attached to a pattern obtained in this way. The Hilbert transform approach is, moreover, subject to problems arising from end effects, which render this approach useless in the presence of strong trends. M-SSA, on the other hand, discriminates between different frequency peaks in a natural, data-adaptive way, and provides a robust decomposition into trends and different oscillatory patterns.

In the present work, we focus on structural changes in single oscillatory patterns that are embedded in possibly high-dimensional systems. To detect such changes, several authors have proposed implementing change-point detection algorithms that are limited, though, to the scalar time series and rely on sequential application of SSA.54–57 

In the context of M-SSA forecasting, Hassani and Mahmoudvand58 have concatenated the channel-wise trajectory matrix Xd into a full trajectory matrix X as well as Y, cf. Eqs. (2) and (7), respectively. These authors concluded that the forecasting performance on X, their VMSSA algorithm, is generally better than on Y, their HMSSA algorithm. This is not surprising, insofar as the VMSSA algorithm uses cross-channel information, and the forecasting procedure can benefit from extra information when the distinct channels share common frequencies. In the present work, we rely on a subsequent varimax rotation in both cases, which improves the identification of common oscillatory behavior and reduces the risk of spurious correlations.59 

To illustrate the insights into phase synchronization that are provided by the standard and the complementary M-SSA algorithm, respectively, we consider a chain of J = 20 Lotka-Volterra models60,61

(15)

These 20 models are diffusively coupled via their x–component, with c0 being the coupling strength. The position in the chain is given by j=1,,J and we assume free boundary conditions x0(n)=x1(n) and xJ+1(n)=xJ(n). Each model in the chain is an extended Rosenzweig-MacArthur version62 of a Lotka-Volterra predator–prey model, with a density-dependent logistic growth G(x)=rx(1xK1) and a Holling-type63 functional response F(x)=αx(1+αhx)1.

The Lotka-Volterra equations have long been used in economic theory as well.64,65 In a highly simplified way, the model captures adjustment delays in the socio-economic system, for example between production and demand or capital and labor. These delays can give rise to endogenous business cycle dynamics.27,66

With the parameters set to r = 0.5, K = 400, h = 0.3, ε=0.3, and μ=0.4, the system undergoes a Hopf bifurcation at α0.0194, after which a periodic solution arises and is stable. The amplitude and period of this limit cycle both increase with increasing α. In our simulation, we choose α=αj to vary linearly with j in the interval [0.0195,0.026], in which the period T varies between 22T25 time units.

The system (15) is integrated with an explicit Runge-Kutta (4,5) solver67 and an initial step size of 103. The solution is sampled at time intervals Δt=1, from which we discard the first 1000 samples to avoid transient behavior, while we keep the time series of length N=5000 samples. In our M-SSA analysis, the window length M = 30 is chosen to cover at least once the longest period of the oscillators.

In a similar experiment on a chain of coupled chaotic Rössler systems, the standard M-SSA algorithm has already proven to provide insightful information about the formation of synchronization clusters.15,24 In the present case of coupled Lotka-Volterra oscillators, synchronization of the distinct limit cycles manifests itself likewise via a clustering of the oscillators. As the coupling strength c increases, the number of clusters decreases: this is reflected by a decreasing number of singular values in the standard M-SSA algorithm in Fig. 4(a).

FIG. 4.

Synchronization for a chain of J = 20 coupled Lotka-Volterra oscillators (15). Spectrum of singular values from (a) the standard M-SSA algorithm, and (b) and (c) from the complementary M-SSA algorithm, with a subsequent complex-EOF analysis of ST-EOFs 1–2; here the latter analysis uses (b) the S-mode, and (c) the T-mode. The M-SSA window length is M = 30 and the 20 leading EOFs are varimax rotated.

FIG. 4.

Synchronization for a chain of J = 20 coupled Lotka-Volterra oscillators (15). Spectrum of singular values from (a) the standard M-SSA algorithm, and (b) and (c) from the complementary M-SSA algorithm, with a subsequent complex-EOF analysis of ST-EOFs 1–2; here the latter analysis uses (b) the S-mode, and (c) the T-mode. The M-SSA window length is M = 30 and the 20 leading EOFs are varimax rotated.

Close modal

In the complementary M-SSA algorithm, the leading oscillatory pair captures always more than 95% of the total variance, irrespective of the coupling strength (not shown); it is this pair that is then analyzed in a subsequent complex EOF analysis. In the S-mode analysis of Fig. 4(b), it is not surprising that the picture resembles qualitatively very well the picture of the standard M-SSA algorithm in Fig. 4(a), since both algorithms seek to simplify the spatial structure. The essential difference is in the representation of oscillatory behavior: an oscillatory pair of two lines in Fig. 4(a) corresponds to a single line of a complex oscillation in Fig. 4(b).

In the T-mode in Fig. 4(c), we likewise observe a cascade of vanishing singular values as the coupling strength increases, although the behavior at successive bifurcation points is somewhat less abrupt. To understand the small but visible differences between the results in the S- and T-mode in Figs. 4(b) and 4(c), respectively, we compare next the underlying spatio-temporal oscillations in greater detail.

Figure 5(a) shows a segment of the observed spatio-temporal behavior at coupling strength c = 0.1. At this coupling strength, the regime of global synchronization has not yet been reached, and the collective oscillation of the elements in the chain is interrupted by phase slips at the border that separates neighboring clusters.

FIG. 5.

Chain of purely periodic oscillators (15) at coupling strength c = 0.1. (a) Observed spatio-temporal oscillation of the x-components; and (b)–(d) partial RCs from a complex EOF analysis in S-mode. Parameter values and details given in Fig. 4.

FIG. 5.

Chain of purely periodic oscillators (15) at coupling strength c = 0.1. (a) Observed spatio-temporal oscillation of the x-components; and (b)–(d) partial RCs from a complex EOF analysis in S-mode. Parameter values and details given in Fig. 4.

Close modal

In the S-mode, the spectrum of singular values in Fig. 4(b) indicates two significant clusters at c = 0.1, while the corresponding partial RCs associated with the complex EOFs 1 and 2 in Figs. 5(b) and 5(c), respectively, show a clear separation into two spatially distinct oscillatory patterns. The combination of the two partial RCs in Fig. 5(d), on the other hand, captures most of the details of the observed spatio-temporal pattern in Fig. 5(a). Note that a similar two-cluster configuration is also found with the standard M-SSA algorithm at c = 0.1 in Fig. 4(a), while the corresponding RCs show patterns similar to that of the S-mode in Figs. 5(b)–5(d) (not shown).

In the T-mode, the spectrum of singular values in Fig. 4(c) likewise indicates two significant clusters at c = 0.1, but the corresponding partial RCs from a reconstruction with the complex EOFs 1 and 2 in Figs. 6(b) and 6(c), respectively, show a different picture. The temporal clustering in the T-mode now separates between epochs of global synchronization in Fig. 6(a), interrupted by epochs of antiphase behavior in Fig. 6(b). As in the S-mode, the combination of the two patterns in the T-mode captures most of the details of the observed spatio-temporal pattern, cf. Fig. 6(c). As the coupling strength increases further, epochs of global synchronization become more frequent, while the corresponding singular value in Fig. 4(c) gradually increases before the onset of complete synchronization at c0.23.

FIG. 6.

Same as Fig. 5, but with partial RCs from a complex EOF analysis in T-mode.

FIG. 6.

Same as Fig. 5, but with partial RCs from a complex EOF analysis in T-mode.

Close modal

In this study, we analyze macroeconomic data from the World Development Indicators (WDI) database of the World Bank.68 The annual datasets provide a comprehensive collection of global development data, from which we select five variables: Gross domestic product (GDP) at market prices, gross fixed capital formation (GDI, formerly gross domestic fixed investment), final consumption expenditure (CON), exports (EXP), and imports (IMP) of goods and services. All variables are in constant 2010 US$.

We restrict ourselves to the interval 1970–2015 of length N = 46 years, for which we have 104 economies with no missing values from at least one of the five macroeconomic indicators above. In this subset of the WDI dataset, each of these 104 economies is thus represented by at least one variable out of the five selected. The total number D = 336 of input time series for our analysis lies, therewith, between the (number of economies = 104) × (5 variables) = 520 and half this number. Additional estimates for the remaining countries and variables with missing values in the WDI dataset complete the global picture, as described in  Appendix B.

In our analysis of the macroeconomic data, we have chosen to separate the shorter-term fluctuations from the long-term trend in a single M-SSA analysis. In contrast to the common idea of first detrending the data,28–30 the single-step M-SSA analysis provides us with a more consistent separation into a permanent trend component and transitory fluctuations that are orthogonal to it.33 

Figure 7 shows the GDP time series of ten major economies, together with their reconstruction by the leading pair of RCs 1–2. This leading pair captures about 99% of the total variance, and it is associated with the growth trend component of the dataset.

FIG. 7.

Annual GDP data of ten major economies from 1970 to 2015 (red) and estimated trend component (black) captured by RCs 1–2 from a complementary M-SSA analysis of the WDI dataset. The window length is M = 12 years and all 12 EOFs are varimax rotated. The vertical bars in all panels indicate NBER-defined US recessions. The trend component captures 99% of the total variance. Note the different scales on the y-axis.

FIG. 7.

Annual GDP data of ten major economies from 1970 to 2015 (red) and estimated trend component (black) captured by RCs 1–2 from a complementary M-SSA analysis of the WDI dataset. The window length is M = 12 years and all 12 EOFs are varimax rotated. The vertical bars in all panels indicate NBER-defined US recessions. The trend component captures 99% of the total variance. Note the different scales on the y-axis.

Close modal

In classical spectral estimation methods, like the Fourier transform, the dominant character of this growth trend component can lead to a strong influence on neighboring low-frequency bands, due to leakage effects.69 For this reason, the dataset is often detrended first, prior to any spectral analysis.35,36 In the present analysis, though, a careful varimax rotation of the M-SSA EOFs helps reduce mixture effects and improves the separation of nearby frequencies.59 

Figure 8 shows the residuals after subtracting the trend component given by RCs 1–2 from the raw data. The success of a reasonable separation between the growth trend component and transitory fluctuations is clearly visible in the GDP trend residuals of the US economy. In this case, the downward fluctuations align very well with the official NBER-defined US recessions,70 while the stylized fact of an asymmetric business cycle—with the recession phase much shorter than the expansion phase—becomes strikingly apparent.71,72

FIG. 8.

GDP trend residuals (red) and reconstruction with RCs 3–4 (black) from the M-SSA analysis in Fig. 7. This RC pair captures 73% of the trend residuals' total variance, while the variance it captures in the GDP of each country is given in the legend of the corresponding panel (in %). Note the different scales on the y-axis.

FIG. 8.

GDP trend residuals (red) and reconstruction with RCs 3–4 (black) from the M-SSA analysis in Fig. 7. This RC pair captures 73% of the trend residuals' total variance, while the variance it captures in the GDP of each country is given in the legend of the corresponding panel (in %). Note the different scales on the y-axis.

Close modal

The ups and downs in the trend residuals of the other countries in Fig. 8 are not aligned that well with the ups and downs in the US trend residuals. This raises therefore the question of whether the transitory fluctuations of various countries are synchronized and, if so, whether the cluster configuration changes over time.

Figure 8 furthermore shows the reconstruction of the trend residuals with the leading oscillatory pair RCs 3–4, whose broad periodicity is of 7–11 year. The cyclic character of successive expansions and recessions is very well reproduced in many countries, while 73% of the residuals' variance is captured by this oscillatory pair.

Besides this main oscillatory mode, we observe two other oscillatory modes of 5–6-year and 3–4-year period in RCs 5–6 and RCs 7–8, respectively. These two modes capture 18% and 5% of the trend residuals' total variance, respectively, and represent not just simply harmonics of the main oscillatory mode, as we would expect from a Fourier decomposition of asymmetric cycles.

The second oscillatory mode of 5–6-year period is shown in Fig. 9, together with a reconstruction of the trend residuals with a combination of the two leading oscillatory modes, given by RCs 3–4 and 5–6. It turns out that the combination of the two modes provides a remarkably good fit to the apparently quite erratic behavior of the trend residuals.

FIG. 9.

Same as Fig. 8 but for a reconstruction with RCs 5–6 (black), as well as the sum of RCs 3–4 and RCs 5–6 (blue). The second RC pair captures 18% of the trend residuals' total variance, while the variance it captures in the GDP of each country is given in the legend of each panel (in %).

FIG. 9.

Same as Fig. 8 but for a reconstruction with RCs 5–6 (black), as well as the sum of RCs 3–4 and RCs 5–6 (blue). The second RC pair captures 18% of the trend residuals' total variance, while the variance it captures in the GDP of each country is given in the legend of each panel (in %).

Close modal

In the US GDP, the RCs 5–6 are characterized by a fairly constant amplitude throughout the whole observation interval. Groth et al.35 have already shown that an oscillatory mode of similar period plays an important role during US recessions, when the trajectory of the dataset stays closer to a suspected limit cycle, like the one in the non-equilibrium dynamic model (NEDyM) of Hallegatte et al.66 or in other endogenous business cycle models.27 

The remarkably persistent character of RCs 5–6 in the US economy is in contrast to a more transient behavior in several of the European countries in Fig. 9. In the latter, RCs 5–6 have a particularly high amplitude during the Great Recession of 2008–2009 and also indicate a smaller second recession afterwards. This behavior has often been called a double-dip recession in many European countries.73 An analysis of the phase relations in RCs 5–6 suggests a leading role of the US in generating this mode, i.e., the US economy leads the European countries by 3–12 month (not shown).

To understand the general characteristics of business cycle synchronization, we focus in Sec. V B on the main oscillatory mode with a near-periodicity of 7–11 year, while more detailed analyses of the other modes are left for future studies.

Note that oscillatory modes of similar periodicity have already been identified in previous M-SSA analyses of the US business cycle35 and of a subset of European countries.36 These studies, however, relied on a prior detrending of the raw data, which is subject to the criticism of business cycles being merely a spurious by-product of the detrending procedure.74 

Objective significance tests, though, have demonstrated that these modes cannot be generated by random shocks alone.35,36 The consistency in the oscillation period between M-SSA analyses of detrended data and raw data is an even stronger indicator for the existence of shared, universal properties of the underlying dynamics. Moreover, multi-annual cyclic behavior, such as the Juglar1 7–11-year cycles and the 3–4-year Kitchin2 cycles, has been discussed at some length before.3 

In the reconstruction of the trend residuals with RCs 3–4, a fairly complex structure of phase-and-amplitude modulations becomes apparent in Fig. 8. To simplify the interpretation, we continue with a subsequent complex EOF analysis of this main oscillatory mode, in which we try to identify structural changes over time.

Figure 10 shows the corresponding temporal patterns of the four leading complex EOFs, cf.  Appendix B. We have chosen the T-mode version of the analysis, in which we seek to identify distinct clusters in time via a varimax rotation of the complex T-EOFs. It turns out that we are indeed able to identify several oscillatory modes that characterize distinct time intervals. This cluster identification will be discussed in greater detail in Secs. V B 1V B 4.

FIG. 10.

Complex EOF analysis in T-mode of ST-EOFs 3–4 that corresponds to RCs 3–4 in Fig. 8. Shown are the complex-valued temporal patterns that characterize the spatio-temporal oscillations in the partial RCs, based on a reconstruction with the four leading complex EOFs. The real and imaginary parts are plotted in blue and red, respectively, while the envelope is plotted in green. The variance captured by each mode is given in the legend of each panel (in %).

FIG. 10.

Complex EOF analysis in T-mode of ST-EOFs 3–4 that corresponds to RCs 3–4 in Fig. 8. Shown are the complex-valued temporal patterns that characterize the spatio-temporal oscillations in the partial RCs, based on a reconstruction with the four leading complex EOFs. The real and imaginary parts are plotted in blue and red, respectively, while the envelope is plotted in green. The variance captured by each mode is given in the legend of each panel (in %).

Close modal

Figure 11 displays rescaled versions of the GDP trend residuals and their corresponding RCs 3–4 for the ten major economies from Fig. 8. In addition, the partial RCs from a reconstruction with the complex EOFs 1–4 are shown in Figs. 11(b)–11(e).

FIG. 11.

GDP trend residuals of ten major economies (in color) and their reconstructions (black) with (a) RCs 3–4, as reproduced from Fig. 8; and (b)–(e) partial RCs 1–4 from the complex EOF analysis in Fig. 10. For ease of comparison, the GDP trend residuals have been rescaled to the interval [100,100].

FIG. 11.

GDP trend residuals of ten major economies (in color) and their reconstructions (black) with (a) RCs 3–4, as reproduced from Fig. 8; and (b)–(e) partial RCs 1–4 from the complex EOF analysis in Fig. 10. For ease of comparison, the GDP trend residuals have been rescaled to the interval [100,100].

Close modal

1. Complex EOF 1

The leading oscillatory mode in complex EOF 1 reaches its amplitude maximum during the Great Recession of 2008–2009, and it remains high until the end of the time interval, cf. Fig. 10(a). In terms of overall impact, the leading complex EOF 1 captures more variance than all the remaining complex EOFs together. Such dominance is in agreement with the view that this recession was the most severe global one during the postwar period.75 The persistently high amplitude after the recession is, moreover, consistent with the observation that the world economy is still struggling with post-crisis adjustments.73 

The corresponding partial RCs in Fig. 11(b) provide a good visual approximation of the trend residuals in many countries during this time interval, a fact that emphasizes the global character of this crisis. For the US economy, the good fit to the GDP trend residuals throughout the full observation interval, together with the ups-and-downs in phase with the NBER-defined recessions, suggest that this mode plays an important role in its dynamics.

The global importance of the US economy on world economic activity raises, therefore, the question of the role that it plays in the generation of this recurrent pattern elsewhere, too. For the other countries, though, we get a more complicated picture, in which the phase and the amplitude of this mode vary from country to country.

To get a better picture of the global pattern, we next plot in Fig. 12 the corresponding phase and amplitude relations for each of the countries on a world map. For each country, the relations among the variables' phase and amplitude are shown in a polar coordinate system, with the two-letter country code at the origin.

FIG. 12.

World map of phase and amplitude relations in complex S-EOF 1. For each country, the relations among the variables' phase and amplitude are shown in a polar coordinate system, with the two-letter country code at the origin. The corresponding variable codes and colors of the pointers are given in the small compass inset at the lower left. Estimates for variables with missing values are indicated by transparent pointers. Phase differences are given with respect to US GDP in a clockwise manner; i.e., positive and negative values indicate a phase that leads or lags the US GDP, respectively. The land area of each country is proportional to its maximum amplitude over all of its variables.

FIG. 12.

World map of phase and amplitude relations in complex S-EOF 1. For each country, the relations among the variables' phase and amplitude are shown in a polar coordinate system, with the two-letter country code at the origin. The corresponding variable codes and colors of the pointers are given in the small compass inset at the lower left. Estimates for variables with missing values are indicated by transparent pointers. Phase differences are given with respect to US GDP in a clockwise manner; i.e., positive and negative values indicate a phase that leads or lags the US GDP, respectively. The land area of each country is proportional to its maximum amplitude over all of its variables.

Close modal

For the 104 economies and those variables with no missing values, the phase and amplitude can be directly derived from the complex S-EOFs in Vc; they are indicated by opaque pointers. The map also shows estimates of the phase and amplitude relations for time series with missing data; these are indicated by transparent pointers. Details of the estimation procedure can be found in  Appendix B.

To guide the reading of the complex world map of phase and amplitude relations, we have furthermore chosen to rescale the land area of each country in proportion to its importance in the oscillatory mode, i.e., proportional to the maximum amplitude of its variables. To calculate this so-called cartogram, we have used the ScapeToad toolkit;76 this toolkit applies a diffusion-based algorithm to find a good balance between the correct size of the land area and low distortion of map regions.77 

In Fig. 12, the prominent role that the US economy plays in the generation of this 7–11-year oscillatory mode becomes immediately apparent in several of its aspects. In terms of amplitude, we see that the US economy contributes a large part to the variance of the mode. This part is about the same order of magnitude as all the European Union countries taken together, and also comparable to the role of China, Japan and Russia taken together. In terms of phase, we see that almost all countries and their variables lag behind the US GDP. The typically stronger link of the United Kingdom (GB) with the US is reflected in a smaller phase difference between the two economies, while other European economies, such as Germany (DE), France (FR), Spain (ES), and Italy (IT), show a larger phase difference.

Furthermore, we see that in the US, the investment sector plays a leading role in this mode, i.e., the GDI (turquoise pointer) leads the GDP (dark blue pointer) by 3–6 months. This is consistent with the situation at the onset of the financial crisis in 2008, when large US investment banks deteriorated most rapidly.75 

On the other hand, one notices that the impact on the German economy was strongest in exports (red pointer), while the GDP impact was smaller and occurred only considerably later. However, the apparently smaller impact on the German GDP arises from the fact that the RCs 3–4 capture a generally smaller part of the German GDP trend residuals during the Great Recession; see again Fig. 8. The otherwise strong impact on the GDP is instead captured by RCs 5–6, as per Fig. 9; this higher-mode impact is consistent with a double-dip recession in many European countries.

In contrast to the strong synchronized decline of economic activity in many countries across the world during the Great Recession, Fig. 12 also indicates that the GDP of China (CN) is in phase opposition to the US GDP. This anti-correlation is interesting as it points to positive effects of the Great Recession on the Chinese economy, and it is also visible in a strong increase of the corresponding partial RCs during this epoch in Fig. 11(b).

2. Complex EOF 2

The second oscillatory mode, captured by complex EOF 2, reaches its maximum amplitude at around 1987, cf. Fig. 10(b). The corresponding partial RCs in Fig. 11(c) provide a good visual fit to the GDP trend residuals of Japan during this epoch. This good fit is interesting insofar as Japan's economy did suffer from an economic bubble around 1987.

The importance of this mode for the dynamics of the Japanese economy is clearly visible in the corresponding map of phase-and-amplitude relations in Fig. 13. The US economy, too, plays a major role in this mode, and the US GDP leads the Japanese GDP by slightly more than a quarter of a cycle, i.e., about 3 years. The corresponding partial RCs of the US GDP, cf. Fig. 11(c), again appear in phase with the NBER-defined recessions, although their relative importance in the explanation of the US GDP trend residuals is much smaller in comparison with the leading oscillatory mode in Fig. 11(b).

FIG. 13.

Same as Fig. 12, but for complex S-EOF 2.

FIG. 13.

Same as Fig. 12, but for complex S-EOF 2.

Close modal

European economies, on the other hand, show a less consistent picture of phase relations in Fig. 13. While the strong link of the United Kingdom to the US is again reflected in a small phase difference between the two economies, the German economy seems to follow Japan in this mode.

3. Complex EOF 3

The next oscillatory mode is captured by complex EOF 3 and it reaches a maximum amplitude around the year 2000, cf. Fig. 10(c), although the amplitude increases again toward the end of the observation interval. The corresponding partial RCs in Fig. 11(d) provide a good visual fit to the GDP trend residuals of China, while this mode's importance for the dynamics of the Chinese economy is also reflected in the phase-and-amplitude map of Fig. 14.

FIG. 14.

Same as Fig. 12, but for complex S-EOF 3.

FIG. 14.

Same as Fig. 12, but for complex S-EOF 3.

Close modal

This mode also seems to capture several consecutive events that have affected the world economy at the turn of the millennium. On the one hand, we see in Fig. 14 that countries in Southeast Asia—like Indonesia (ID), Thailand (TH), and Singapore (SG)—and in the Far East, like South Korea (KR), show a phase lead of 2.5 years prior to the maximum around the year 2000, a lead that could be attributed to the Asian financial crisis in 1997. The latter crisis is known to have started in Thailand, a timing that is reflected in Thailand's phase lead of about 2.5 years in the figure. The figure shows, moreover, that the Japanese (JP) economy was also affected, without any significant phase lag, maybe because of its strong export links with these Asian economies.

The US economy, on the other hand, lags these Asian economies by about two years. The corresponding partial RCs of the US GDP in Fig. 11(d) again appear in phase with the NBER recessions, and they show a strong decline during the 2001 recession. Figure 14 also suggests later effects of this mode in many European economies, with a particularly strong impact on the German one.

4. Complex EOF 4

Finally, the present synchronization analysis yields an oscillatory mode in complex EOF 4 that reaches its maximum amplitude during the time interval 1970–1980, as seen from Fig. 10(d). Among the ten major economies shown in Fig. 11(e), the partial RCs provide a good visual fit to the GDP trend residuals of Germany during this interval. It is remarkable that the strong downward swing in the partial RCs agrees in its timing with the 1979 oil shock.

Despite the small visual importance of the partial RCs for the US GDP in this mode, the corresponding phase-and-amplitude map in Fig. 15 indicates that the absolute effects of this mode on the US economy are still quite large. It appears that its impact was particularly strong on US consumption (purple pointers in the figure), while the impacts on German and Japanese consumption, but also on that of other countries, like Mexico (MX), become visible as well. With respect to the chronology of the 1979 oil shock, it is interesting to note that the mode is characterized by a large lead in Iran's (IR) exports, a lead that is consistent with a preceding strong reduction in Iran's oil production (not shown).

FIG. 15.

Same as Fig. 12, but for complex S-EOF 4.

FIG. 15.

Same as Fig. 12, but for complex S-EOF 4.

Close modal

In this section, we have tried to assess the evidence for synchronized business cycle activity across countries, regions, and the world. Relying on the complementary M-SSA algorithm introduced in Sec. III C, we were able to identify a major oscillatory mode of 7–11-year periodicity that already captures 73% of the trend residuals' variance, while a second oscillatory mode of 5–6-year period captures another 18%, for a combined total of 91%.

The combination of these two modes thus provides a remarkably good reconstruction of the apparently rather erratic behavior of the trend residuals. Despite the appealing simplicity that the low-order approximation by just two such modes involves, the details of the corresponding spatio-temporal structure are still quite intricate. In a subsequent complex EOF analysis, we have therefore tried to break down this intricate structure into simpler components.

In a complex EOF analysis in T-mode, as outlined in Sec. III E, it turns out that the behavior of the leading oscillatory mode with the 7–11-year period can be essentially decomposed into four distinct patterns. Each of the four is more or less associated with a certain subinterval of time in Fig. 10. This clustering in time gave birth to the heuristic idea of the so-called “snapshots,” which characterize the state-space evolution over time.47 The individual snapshots, though, exhibit a rather complex structure of phase and amplitude relations, which is consistent with fairly sudden structural changes between the snapshots. These results flesh out, therefore, the rather complex nature of business-cycle synchronization as time evolves.78 

The idea of snapshots associated with distinct time intervals is only partially valid, though, according to a more careful inspection of Fig. 10. Note, for example, the rather persistent character of the oscillation captured by the leading complex EOF in Fig. 10(a), which has a non-vanishing amplitude throughout the full observation interval. This oscillatory mode seems, therewith, to provide substantial evidence for an objective, quantifiable world business cycle, as suggested, for instance, by Kose et al.8 The mode's strong increase toward the end of the time interval we analyzed is likely to be driven by an increasing global integration of markets, an integration that may also have led to a globally synchronized recession in many regions of the world, cf. Fig. 12.

The present analysis was conducted with no particular connection to the arguments about the causal mechanisms of the cycles that we detected, like the endogenous vs. exogenous nature of business cycles. Be that as it may, the phase and amplitude relations in Fig. 12 render a remarkably diverse picture of the distinct countries' state-space evolution and of the lag-and-lead ordering in their macroeconomic variables. While some highly idealized models—like the NEDyM model,66 for instance—have had remarkable success in reproducing certain stylized facts of business cycles in an isolated economy, the generalization of these results to groups of countries and the inclusion of possible synchronization among them remains a subject of ongoing research.

As part of a better understanding of economic synchronization, the complex EOFs 2–4 of lower variance can also be seen as smaller transient shocks, superimposed on the world business cycle of complex EOF 1. Several of the features discussed in connection with the snapshots of Figs. 13–15 appear to be associated with sector-specific and region-specific synchronization effects.

In addition to the analysis in T-mode, we have furthermore analyzed the M-SSA ST-EOFs 3–4 using a complex EOF analysis in S-mode (not shown). Such an S-mode analysis tends to bring out regional and sectorial clusters, cf. Sec. III E. In using it, we were able to identify three main oscillatory modes, each of which is dominated by the temporal behavior of the US, China, and Japan, respectively. The leading US mode captures 66% of the variance, while the two other modes capture only 11% each. The corresponding phase-and-amplitude map of the leading S-mode (not shown) resembles in many details that of the leading T-mode in Fig. 12, and thus provides further evidence for the leading role of the US economy in the generation of a world business cycle.

Beside these three S-modes, there was no further indication of smaller, supranational clusters in this leading mode of 7–11-year period, for example on the European Union level. This finding is consistent with other analyses, in which global factors play a major role in the business cycle dynamics in most countries, while region-specific factors play a more minor role.8,79 On the other hand, the question of whether other higher-frequency modes reflect more region-specific factors will be left for future studies.

In the present paper, we have studied common dynamical properties of business-cycle fluctuations across countries and macroeconomic indices. In a large sample of over 100 economies that represent economic regions from all around the world, we were able to identify shared mechanisms and common spectral properties of business cycle activity.

To identify shared oscillatory modes, we relied here on the advanced spectral methodology of multivariate singular spectrum analysis (M-SSA; Ghil et al.,10 Alessio,37 and references therein). This methodology has already proven its efficiency and accuracy in the spatio-temporal analysis of large datasets and in the analysis of phase synchronization.

In the present paper, we have furthermore introduced a modification of the M-SSA approach to identify changes in the spatio-temporal structure of oscillatory modes. We proposed herein a complementary M-SSA algorithm that is flexible enough to reconstruct oscillatory behavior with structural changes in the space and time domain. A subsequent complex EOF analysis49 helps the systematic decomposition of the complex spatio-temporal structure into simpler oscillatory patterns.

We have discussed two variants of simple-structure rotation for complex EOFs. In the so-called S-mode, a rotation is applied to the complex S-EOFs and results in spatially distinct clusters, in a manner that is similar to varimax M-SSA analysis.15,23 In the alternative T-mode, a rotation is applied to the complex T-EOFs and yields a complementary view of distinct clusters in time. In a numerical study of a chain of coupled oscillators, we have shown that—in the case of intermittent synchronization, in which global synchronization is not yet reached—the T-mode rotation separates time intervals of global synchronization from those of asynchronous, out-of-phase behavior.

These algorithmic developments and their numerical results were applied next to our macroeconomic data. The key finding is the identification of a common mode of business cycle activity that is shared by many economies and their individual macroeconomic indicators. The leading mode of variability, with its 7–11-year period, has a persistent component that supports the existence of a world business cycle. Superimposed on this persistent component, we were able to identify further components that capture substantial fractions of variability over distinct time intervals. These lower-variance modes can be linked to several major events in the world economic activity of the postwar era.

It is a pleasure to thank a reviewer for an admirably thorough and constructive report, which greatly helped clarify the exposition. The research herein was supported by Grant No. N00014-16-1-2073 from the Multidisciplinary University Research Initiative (MURI) of the Office of Naval Research and by National Science Foundation Grant No. OCE-1243175 (A.G. and M.G.), as well as by the KIC-WINnERS project of the European Union (A.G.) and the Energy and Prosperity Chair of the Institut Louis Bachelier (M.G.).

Groth and Ghil15 have proposed a modified version of the varimax algorithm80 for the M-SSA eigenvectors, in order to simplify their dynamical interpretation and to reduce mixture effects. In this version, only the spatial variance between the channels Ed of E in Eq. (4) is maximized. Let edk(m) be the m-th row element of the k-th column in segment Ed. Prior to the calculation of the varimax criterion in each rotation step, one computes the participation index

(A1)

of channel d to the k-th ST-EOF.

The varimax algorithm then seeks to maximize the quadratic functional

(A2)

subject to the normalization hd=k=1Sπdk. The varimax algorithm's original simplicity of pairwise rotations of eigenvectors is also maintained in the maximization of VM. Alternatively, a closed matrix formulation of the algorithm for the simultaneous rotation of all S eigenvectors by iterative SVD decompositions has been proposed.41 Whichever algorithm is used to maximize VM, the resulting RCs typically tend to yield a unimodal, narrowband power spectral density that is clearly associated with a unique frequency.59 

In the complementary M-SSA algorithm, we likewise seek to improve the separability of distinct oscillations by a subsequent varimax rotation of V. Let vdk(n) be the n-th row element of the k-th column in segment Vd. The modified varimax algorithm then seeks to maximize the functional (A2) for the participation index

(A3)

In the standard M-SSA algorithm, a helpful tool for the understanding of the spatio-temporal dynamics of the reconstructed oscillatory behavior is that of a phase-composite analysis,40 in which the SVD of the corresponding RCs is calculated. A phase is then defined as the argument of (i) the leading PC, expanded into the complex plane via its first derivative;81 or of (ii) the leading two PCs in phase-quadrature.44 In either case, the oscillatory behavior is sufficiently well represented in a two-dimensional subspace, in which the oscillation is described as the outer product of an instantaneous and a spatial phase.44 

In the complementary M-SSA analysis of Sec. III C, however, the RCs can exhibit more complex oscillatory behavior, as seen, for example, in Fig. 3(b). A subsequent complex-EOF analysis, cf. Sec. III E, then yields a further decomposition of the complex oscillatory behavior into a set of oscillations with simpler spatio-temporal structures. According to Eq. (13), this simplification is understood in the sense that the oscillations can be described by the outer product of complex T-EOFs and S-EOFs.

Although this procedure is reminiscent of the phase-composite analysis—with the complex S-EOFs already providing spatial phase information—the complex T-EOFs have reduced length N only. Since Eq. (13) provides a factorization of Z, the complex T-EOFs do not contain direct phase information within the window length M.

It is only the RCs that capture the phase of the time series in a well-defined least-squares sense, and a possible solution to this problem could be to first calculate the corresponding partial RCs from Eqs. (10) and (14), and then the SVD.

A mathematically more elegant solution, without the need for an additional SVD, starts by replacing T-EOF vp and vq in (10) with the real and imaginary part of the k-th complex T-EOF in Uc, respectively. We obtain, accordingly, two matrices Rp and Rq, each of size M×N, while averaging along the skew diagonals yields the real and imaginary part of a new complex vector u of length N. The oscillation in the k-th partial RC is finally given by the outer product of u and the corresponding k-th complex S-EOF in Vc, as seen, for example, in Figs. 10 and 12, respectively.

To obtain estimates of the spatial phase and amplitude for additional time series with missing values, we start by setting all missing values to zero. Next, we project the filled-in time series onto the oscillatory pair of “filter” EOFs in the matrix U obtained by the complementary M-SSA algorithm. This projection yields estimates for vp and vq, which are then projected onto Uc from the complex-EOF analysis to finally obtain estimates of Vc, as done, for example, in order to obtain Fig. 12.

1.
C.
Juglar
,
Des Crises Commerciales Et Leur Retour Périodique En France, En Angleterre, Et Aux Etats Unis
(
Guillaumin
,
Paris
,
1862
).
2.
J.
Kitchin
, “
Cycles and trends in economic factors
,”
Rev. Econ. Stat.
5
,
10
16
(
1923
).
3.
A. F.
Burns
and
W. C.
Mitchell
,
Measuring Business Cycles
(
National Bureau of Economic Research
,
New York
,
1946
).
4.
G.
Soros
, “
The crisis & what to do about it
,”
Real-World Econ. Rev.
48
,
312
318
(
2008
);
G.
Soros
, “
The crisis & what to do about it
,”
The New Paradigm for Financial Markets: The Credit Crisis of 2008 and What It Means
(
Public Affairs
,
2008
).
5.
P.
Romer
, “
The trouble with macroeconomics
,”
Am. Econ.
(to be published).
6.
M.
Artis
and
W.
Zhang
, “
International business cycles and the ERM: Is there a European business cycle?
,”
Int. J. Finance Econ.
2
,
1
16
(
1997
).
7.
B.
Süssmuth
, “
National and supranational business cycles (1960-2000): A multivariate description of central G7 and Euro15 NIPA aggregates
,”
Working Paper 658
(
CESifo
,
2002
).
8.
M. A.
Kose
,
C. M.
Otrok
, and
C. H.
Whiteman
, “
International business cycles: World, region, and country-specific factors
,”
Am. Econ. Rev.
93
,
1216
1239
(
2003
).
9.
J.
de Haan
,
R.
Inklaar
, and
R.
Jong-A-Pin
, “
Will business cycles in the Euro area converge? A critical survey of empirical research
,”
J. Econ. Surv.
22
,
234
273
(
2008
).
10.
M.
Ghil
,
M. R.
Allen
,
M. D.
Dettinger
,
K.
Ide
,
D.
Kondrashov
,
M. E.
Mann
,
A. W.
Robertson
,
A.
Saunders
,
Y.
Tian
,
F.
Varadi
, and
P.
Yiou
, “
Advanced spectral methods for climatic time series
,”
Rev. Geophys.
40
,
1
41
(
2002
).
11.
R.
Vautard
and
M.
Ghil
, “
Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series
,”
Physica D
35
,
395
424
(
1989
).
12.
M.
Ghil
and
R.
Vautard
, “
Interdecadal oscillations and the warming trend in global temperature time series
,”
Nature
350
,
324
327
(
1991
).
13.
D.
Broomhead
,
R.
Jones
, and
G.
King
, “
Topological dimension and local coordinates from time series data
,”
J. Phys. A
20
,
L563
(
1987
).
14.
M.
Paluš
and
I.
Dvorák
, “
Singular-value decomposition in attractor reconstruction: Pitfalls and precautions
,”
Physica D
55
,
221
234
(
1992
).
15.
A.
Groth
and
M.
Ghil
, “
Multivariate singular spectrum analysis and the road to phase synchronization
,”
Phys. Rev. E
84
,
036206
(
2011
).
16.
S.
Boccaletti
,
J.
Kurths
,
G.
Osipov
,
D.
Valladares
, and
C.
Zhou
, “
The synchronization of chaotic systems
,”
Phys. Rep.
366
,
1
101
(
2002
).
17.
A.
Pikovsky
,
M.
Rosenblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
2003
).
18.
G. V.
Osipov
,
J.
Kurths
, and
C.
Zhou
,
Synchronization in Oscillatory Networks
(
Springer
,
Berlin/Heidelberg
,
2007
).
19.
M. G.
Rosenblum
,
A. S.
Pikovsky
, and
J.
Kurths
, “
Phase synchronization of chaotic oscillators
,”
Phys. Rev. Lett.
76
,
1804
1807
(
1996
).
20.
G.
Osipov
,
B.
Hu
,
C.
Zhou
,
M.
Ivanchenko
, and
J.
Kurths
, “
Three types of transition to phase synchronization in coupled chaotic oscillators
,”
Phys. Rev. Lett.
91
,
024101
(
2003
).
21.
I. Z.
Kiss
,
Q.
Lv
, and
J. L.
Hudson
, “
Synchronization of non-phase-coherent chaotic electrochemical oscillations
,”
Phys. Rev. E
71
,
035201
(
2005
).
22.
K.
Pukenas
, “
Detecting phase synchronization in coupled oscillators by combining multivariate singular spectrum analysis and fast factorization of structured matrices
,”
J. Vibroeng.
16
,
2624
2630
(
2014
).
23.
A.
Groth
and
M.
Ghil
, “
Monte Carlo singular spectrum analysis (SSA) revisited: Detecting oscillator clusters in multivariate datasets
,”
J. Clim.
28
,
7873
7893
(
2015
).
24.
L. L.
Portes
and
L. A.
Aguirre
, “
Enhancing multivariate singular spectrum analysis for phase synchronization: The role of observability
,”
Chaos
26
,
093112
(
2016
).
25.
L. G.
Arnold
,
Business Cycle Theory
(
Oxford University Press
,
2002
).
26.
R. G.
King
and
S. T.
Rebelo
, “
Resuscitating real business cycles
,”
Working Paper 7534
(
National Bureau of Economic Research
,
2000
).
27.
C.
Chiarella
,
P.
Flaschel
, and
R.
Franke
,
Foundations for a Disequilibrium Theory of the Business Cycle
(
Cambridge University Press
,
2005
).
28.
R. J.
Hodrick
and
E. C.
Prescott
, “
Postwar U.S. business cycles: An empirical investigation
,”
J. Money Credit Bank
29
,
1
16
(
1997
).
29.
F.
Canova
, “
Detrending and business cycle facts
,”
J. Monet. Econ.
41
,
475
512
(
1998
).
30.
M.
Baxter
and
R. G.
King
, “
Measuring business cycles: Approximate band-pass filters for economic time series
,”
Rev. Econ. Stat.
81
,
575
593
(
1999
).
31.
U.
Woitek
, “
The G7-countries: A multivariate description of business cycle stylized facts
,” in
Dynamic Disequilibrium Modelling: Theory and Applications
, edited by
W.
Barnett
,
G.
Gandolfo
, and
C.
Hillinger
(
Cambridge University Press
,
1996
), pp.
283
309
.
32.
A.
Dickerson
,
H.
Gibson
, and
E.
Tsakalotos
, “
Business cycle correspondence in the European Union
,”
Empirica
25
,
49
75
(
1998
).
33.
M.
de Carvalho
,
P. C.
Rodrigues
, and
A.
Rua
, “
Tracking the US business cycle with a singular spectrum analysis
,”
Econ. Lett.
114
,
32
35
(
2012
).
34.
L.
Sella
and
R.
Marchionatti
, “
On the cyclical variability of economic growth in Italy, 1881–1913: A critical note
,”
Cliometrica
6
,
307
328
(
2012
).
35.
A.
Groth
,
M.
Ghil
,
S.
Hallegatte
, and
P.
Dumas
, “
The role of oscillatory modes in U.S. business cycles
,”
OECD J. Bus. Cycle Meas. Anal.
2015
,
63
81
.
36.
L.
Sella
,
G.
Vivaldo
,
A.
Groth
, and
M.
Ghil
, “
Economic cycles and their synchronization: A comparison of cyclic modes in three European countries
,”
J. Bus. Cycle Res.
12
,
25
48
(
2016
).
37.
S. M.
Alessio
,
Digital Signal Processing and Spectral Analysis for Scientists: Concepts and Applications
(
Springer
,
2016
).
38.
D. S.
Broomhead
and
G. P.
King
, “
Extracting qualitative dynamics from experimental data
,”
Physica D
20
,
217
236
(
1986
).
39.
D. S.
Broomhead
and
G. P.
King
, “
On the qualitative analysis of experimental dynamical systems
,” in
Nonlinear Phenomena and Chaos
, edited by
S.
Sarkar
(
Adam Hilger
,
Bristol, England
,
1986
), pp.
113
144
.
40.
G.
Plaut
and
R.
Vautard
, “
Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere
,”
J. Atmos. Sci.
51
,
210
236
(
1994
).
41.
L. L.
Portes
and
L. A.
Aguirre
, “
Matrix formulation and singular-value decomposition algorithm for structured varimax rotation in multivariate singular spectrum analysis
,”
Phys. Rev. E
93
,
052216
(
2016
).
42.
M. R.
Allen
and
A. W.
Robertson
, “
Distinguishing modulated oscillations from coloured noise in multivariate datasets
,”
Clim. Dyn.
12
,
775
784
(
1996
).
43.
E.
Anderson
,
Z.
Bai
,
C.
Bischof
,
S.
Blackford
,
J.
Demmel
,
J.
Dongarra
,
J. D.
Croz
,
A.
Greenbaum
,
S.
Hammarling
,
A.
McKenney
, and
D.
Sorensen
,
LAPACK Users' Guide
, 3rd ed. (
Society for Industrial and Applied Mathematics
,
Philadelphia, PA
,
1999
).
44.
A.
Groth
,
Y.
Feliks
,
D.
Kondrashov
, and
M.
Ghil
, “
Interannual variability in the North Atlantic ocean's temperature field and its association with the wind stress forcing
,”
J. Clim.
30
,
2655
2678
(
2017
).
45.
R.
Vautard
,
P.
Yiou
, and
M.
Ghil
, “
Singular-spectrum analysis: A toolkit for short, noisy chaotic signals
,”
Physica D
58
,
95
126
(
1992
).
46.
M. R.
Allen
and
L. A.
Smith
, “
Monte Carlo SSA: Detecting irregular oscillations in the presence of colored noise
,”
J. Clim.
9
,
3373
3404
(
1996
).
47.
R. W.
Preisendorfer
and
C. D.
Mobley
,
Principal Component Analysis in Meteorology and Oceanography
(
Elsevier
,
Amsterdam
,
1988
).
48.
I. T.
Jolliffe
,
Principal Component Analysis
, 2nd ed. (
Springer
,
New York
,
2002
).
49.
J. D.
Horel
, “
Complex principal component analysis: Theory and examples
,”
J. Clim. Appl. Meteorol.
23
,
1660
1673
(
1984
).
50.
M. B.
Richman
, “
Rotation of principal components
,”
Int. J. Climatol.
6
,
293
335
(
1986
).
51.
R. B.
Cattell
,
Factor Analysis
(
Harper and Row
,
New York
,
1952
).
52.
M. A.
Salles
,
P. O.
Canziani
, and
R. H.
Compagnucci
, “
The spatial and temporal behaviour of the lower stratospheric temperature over the Southern Hemisphere: The MSU view. Part II: Spatial behaviour
,”
Int. J. Climatol.
21
,
439
454
(
2001
).
53.
R.
Huth
,
C.
Beck
,
A.
Philipp
,
M.
Demuzere
,
Z.
Ustrnul
,
M.
Cahynová
,
J.
Kyselý
, and
O. E.
Tveito
, “
Classifications of atmospheric circulation patterns
,”
Ann. N. Y. Acad. Sci.
1146
,
105
152
(
2008
).
54.
V.
Moskvina
and
A.
Zhigljavsky
, “
An algorithm based on singular spectrum analysis for change-point detection
,”
Commun. Stat. Simul. Comput.
32
,
319
352
(
2003
).
55.
T.
Idé
and
K.
Inoue
, “
Knowledge discovery from heterogeneous dynamic systems using change-point correlations
,”
in Proceedings of the 2005 SIAM International Conference on Data Mining
(Society for Industrial & Applied Mathematics SIAM,
2005
), pp.
571
575
.
56.
Y.
Mohammad
and
T.
Nishida
, “
On comparing SSA-based change point discovery algorithms
,” in
Proceedings of the 2011 IEEE/SICE International Symposium on System Integration SII
(Institute of Electrical and Electronics Engineers IEEE,
2011
), pp.
938
945
.
57.
N.
Itoh
and
N.
Marwan
, “
An extended singular spectrum transformation (SST) for the investigation of Kenyan precipitation data
,”
Nonlinear Process. Geophys.
20
,
467
481
(
2013
).
58.
H.
Hassani
and
R.
Mahmoudvand
, “
Multivariate singular spectrum analysis: A general view and new vector forecasting approach
,”
Int. J. Energy Stat.
1
,
55
83
(
2013
).
59.
Y.
Feliks
,
A.
Groth
,
A. W.
Robertson
, and
M.
Ghil
, “
Oscillatory climate modes in the Indian Monsoon, North Atlantic and Tropical Pacific
,”
J. Clim.
26
,
9528
9544
(
2013
).
60.
A. J.
Lotka
,
Elements of Physical Biology
(
Williams & Wilkins
,
Baltimore
,
1926
).
61.
V.
Volterra
,
Leçons Sur La Théorie Mathématique De La Lutte Pour La Vie
(
Gauthier-Villars
,
Paris
,
1931
).
62.
M. L.
Rosenzweig
and
R. H.
MacArthur
, “
Graphical representation and stability conditions of predator-prey interactions
,”
Am. Nat.
97
,
209
223
(
1963
).
63.
C. S.
Holling
, “
The components of predation as revealed by a study of small-mammal predation of the European pine sawfly
,”
Can. Entomol.
91
,
293
320
(
1959
).
64.
R.
Goodwin
, “
A growth cycle
,” in
Socialism, Capitalism and Economic Growth
, edited by
C.
Feinstein
(
Cambridge University Press
,
1967
), pp.
54
58
.
65.
P. A.
Samuelson
, “
Generalized predator-prey oscillations in ecological and economic equilibrium
,”
Proc. Natl. Acad. Sci.
68
,
980
983
(
1971
).
66.
S.
Hallegatte
,
M.
Ghil
,
P.
Dumas
, and
J.-C.
Hourcade
, “
Business cycles, bifurcations and chaos in a neo-classical model with investment dynamics
,”
J. Econ. Behav. Organ.
67
,
57
77
(
2008
).
67.
L. F.
Shampine
and
M. W.
Reichelt
, “
The MATLAB ODE suite
,”
SIAM J. Sci. Comput.
18
,
1
22
(
1997
).
68.
See http://databank.worldbank.org the WDI database.
69.
C. W. J.
Granger
, “
The typical spectral shape of an economic variable
,”
Econometrica
34
,
150
161
(
1966
).
70.
See http://www.nber.org/cycles.html for the US business cycle reference dates.
71.
D. E.
Sichel
, “
Business cycle asymmetry: A deeper look
,”
Econ. Inquiry
31
,
224
236
(
1993
).
72.
J.
Morley
and
J.
Piger
, “
The asymmetric business cycle
,”
Rev. Econ. Stat.
94
,
208
221
(
2012
).
73.
United Nations
,
World Economic Situation and Prospects 2013
(
UN
,
New York
,
2013
).
74.
C. R.
Nelson
and
H.
Kang
, “
Spurious periodicity in inappropriately detrended time series
,”
Econometrica
49
,
741
751
(
1981
).
75.
International Monetary Fund
,
World Economic Outlook: Crisis and Recovery
, World Economic and Financial Surveys (
IMF
,
Washington, DC
,
2009
).
76.
See http://scapetoad.choros.ch for the ScapeToad toolkit.
77.
M. T.
Gastner
and
M. E. J.
Newman
, “
Diffusion-based method for producing density-equalizing maps
,”
Proc. Natl. Acad. Sci.
101
,
7499
7504
(
2004
).
78.
Note that this concept of snapshots is quite distinct from the one used in the theory of nonautonomous and random dynamical systems and of their pullback and random attractors; see, for instance, Ghil, Chekroun, and Simonnet,82 Ghil83, and references therein.
79.
M. J.
Artis
, “
Is there a European business cycle
?,”
Working Paper 1053
(
CESifo
,
2003
).
80.
H.
Kaiser
, “
The varimax criterion for analytic rotation in factor analysis
,”
Psychometrika
23
,
187
200
(
1958
).
81.
V.
Moron
,
R.
Vautard
, and
M.
Ghil
, “
Trends, interdecadal and interannual oscillations in global sea-surface temperatures
,”
Clim. Dyn.
14
,
545
569
(
1998
).
82.
M.
Ghil
,
M. D.
Chekroun
, and
E.
Simonnet
, “
Climate dynamics and fluid mechanics: Natural variability and related uncertainties
,”
Phys. D
237
,
2111
2126
(
2008
).
83.
M.
Ghil
, “
The wind-driven ocean circulation: Applying dynamical systems theory to a climate problem
,”
Discr. Cont. Dyn. Syst. A
37
,
189
228
(
2017
).