Multiscale Granger causality

TL;DR

This paper extends Granger causality to analyze directed interactions across multiple temporal scales in complex systems, using state space models to improve estimation accuracy and reveal multiscale information transfer patterns.

Contribution

It introduces a theoretical framework for multiscale Granger causality using state space models, addressing practical issues like filtering and downsampling.

Findings

State space approach yields higher estimation accuracy than traditional AR models.

Multiscale GC reveals meaningful information transfer patterns in climate data.

Method effectively captures dynamics across multiple time scales.

Abstract

In the study of complex physical and biological systems represented by multivariate stochastic processes, an issue of great relevance is the description of the system dynamics spanning multiple temporal scales. While methods to assess the dynamic complexity of individual processes at different time scales are well-established, multiscale analysis of directed interactions has never been formalized theoretically, and empirical evaluations are complicated by practical issues such as filtering and downsampling. Here we extend the very popular measure of Granger causality (GC), a prominent tool for assessing directed lagged interactions between joint processes, to quantify information transfer across multiple time scales. We show that the multiscale processing of a vector autoregressive (AR) process introduces a moving average (MA) component, and describe how to represent the resulting ARMA process using state space (SS) models and to combine the SS model parameters for computing exact GC values at arbitrarily large time scales. We exploit the theoretical formulation to identify peculiar features of multiscale GC in basic AR processes, and demonstrate with numerical simulations the much larger estimation accuracy of the SS approach compared with pure AR modeling of filtered and downsampled data. The improved computational reliability is exploited to disclose meaningful multiscale patterns of information transfer between global temperature and carbon dioxide concentration time series, both in paleoclimate and in recent years.

Figures

Equations

$$F_{i\to j}=\ln \frac{{\lambda }_{j|jk}}{{\lambda }_{j|ijk}}.$$

$$X_{n+1}$$

$$Y_n$$

$$Z_{n+1}$$

$$Y_n$$

$$\mathbf{P}$$

$$-(\mathbf{AP}{\mathbf{C}}^T+\mathbf{\Theta })(\mathbf{CP}{\mathbf{C}}^T+\mathbf{\Psi }{)}^{-1}(\mathbf{CP}{\mathbf{A}}^T+{\mathbf{\Theta }}^T),$$

$$\mathbf{\Phi }$$

$$\mathbf{K}$$

$$Y_n^{(jk)}={\mathbf{C}}^{(jk)}{Z}_n+E_n^{(jk)}$$

$$Y_n=\sum _{k=1}^p{\mathbf{A}}_k{Y}_{n-k}+U_n,$$

$${\bar{y}}_{m,n}=\sum _{l=0}^q{b}_l{y}_{m,n\tau -l}.$$

$${\tilde{Y}}_n$$

$${\bar{Y}}_n$$

$${\tilde{Y}}_n=\sum _{k=1}^p{\mathbf{A}}_k{\tilde{Y}}_{n-k}+\sum _{l=0}^q{\mathbf{B}}_l{U}_{n-l}$$

$$\tilde{\mathbf{C}}=\left[\begin{array}{cccccc}{\mathbf{A}}_1& \cdots & {\mathbf{A}}_p& {\mathbf{B}}_1& \cdots & {\mathbf{B}}_q\end{array}\right],$$

$$\tilde{\mathbf{A}}=\left[\begin{array}{ccc}\tilde{\mathbf{C}}& & \\ {\mathbf{I}}_{M(p-1)}& 0_{M(p-1)\times M(q+1)}& \\ 0_{M\times M(p+q)}& & \\ 0_{M(q-1)\times Mp}& {\mathbf{I}}_{M(q-1)}& 0_{M(q-1)\times M}\end{array}\right],$$

$$\tilde{\mathbf{K}}={\left[\begin{array}{cccc}{\mathbf{I}}_M& 0_{M\times M(p-1)}& {\mathbf{B}}_0^{-T}& 0_{M\times M(q-1)}\end{array}\right]}^T,$$

$${\mathbf{\Theta }}_{\tau }$$

$${\mathbf{\Xi }}_{\tau }$$

$${\mathbf{\Xi }}_1$$

$$y_{1,n}$$

$$y_{2,n}$$

$$x_{1,n}$$

$$x_{2,n}$$

$$z_{1,n}$$

$$z_{2,n}$$

$$y_{1,n}$$

$$y_{2,n}$$

Full text

††preprint: APS/123-QED

I Introduction

Granger causality (GC) is a powerful tool for assessing directional interactions from time series data according to the notion of time lagged influence first proposed by Wiener Wiener (1956) and then formalized by Granger Granger (1969) and Geweke Geweke (1982, 1984) in the framework of vector autoregressive (AR) modeling of stochastic processes. Since its formulation, GC has gained increasing popularity and is nowadays ubiquitously employed in several scientific fields ranging from econometrics to social and climate sciences, neuroscience and physiology Freeman (1983); Smirnov and Mokhov (2009); Bressler and Seth (2011); Porta and Faes (2015). The great success of this measure comes from its conceptual simplicity, data-driven nature, and relative ease of implementation. An additional appealing property of GC is the principled interpretation of its generalized probabilistic formulation, which is closely related to the the information-theoretic concept of transfer entropy Barnett et al. (2009).

Many processes in physics, biology and other fields exhibit dynamics spanning multiple temporal scales Ivanov et al. (1999); Kang et al. (2009); Valencia et al. (2009); Chou (2011); Wang et al. (2013). The multiscale properties of an observed stochastic process can be explored first resampling at different temporal scales the originally measured realization of the process, and then assessing the dynamical complexity of the rescaled series Costa et al. (2002). This approach has been followed with great success to quantify the multiscale behavior of the individual dynamics of scalar processes Kang et al. (2009); Valencia et al. (2009); Chou (2011); Wang et al. (2013). However, its extension to the multiscale computation of the information transfer between processes, though attempted in empirical studies Lungarella et al. (2007); Paluš (2014), is far less straightforward. In fact, the multiscale evaluation of lagged influence measures such as the GC is severely complicated by theoretical and practical issues Barnett and Seth (2015); Solo (2016). These issues arise from the rescaling procedure, which essentially consists in a filtering step eliminating the fast temporal scales (classically performed by averaging Costa et al. (2002)) followed by a downsampling step that coarse-grains the time series around the selected scale. The filtering step leaves theoretically unchanged the GC values, but degrades severely their estimation affecting reliability, stability and data demand Barnett and Seth (2011). The downsampling step is even more problematic, as it alters GC values in a way that was unknown until very recently Barnett and Seth (2015); Solo (2016) and impacts consistently detectability and accuracy of GC estimates. For these reasons, research on the data-driven inference of the multiscale structure of coupled processes, though holding a great potential, is still largely undeveloped.

In order to provide a formal extension of GC to multiscale analysis, here we introduce for the first time an analytical frame for its computation on linear multivariate stochastic processes subjected to averaging and downsampling. We exploit the theory of state space (SS) models Aoki and Havenner (1991); Barnett and Seth (2015); Solo (2016) to yield exact GC values for for coupled processes observed at different time scales. The high computational reliability of the associated multiscale GC estimator is demonstrated in simulated AR processes and then used to explore, spanning a wide range of time scales, the patterns of information transfer between anthropogenic emissions and global temperatures.

II Granger Causality

II.1 Definition

To lay the groundwork for multiscale GC computation and introduce notations, we start resuming the calculation of GC for general multivariate processes Granger (1969); Geweke (1984). Let us consider a discrete-time, stationary vector stochastic process composed of M real-valued zero-mean scalar processes, $Y_{n}=[y_{1,n}\cdots y_{M,n}]^{T}$, $-\infty<n<\infty$, and assume the process $y_{j}$ as the target and the process $y_{i}$ as the driver (the remaining $M-2$ processes form the vector $Y_{k}$, where $k=\{1,\ldots,M\}\backslash\{i,j\}$). Then, denoting the present and the past of vector and scalar variables respectively as $Y_{n}$, $y_{n}$, and $Y_{n}^{-}=[Y_{n-1}^{T}Y_{n-2}^{T}\cdots]$, $y_{n}^{-}=[y_{n-1}y_{n-2}\cdots]$, GC from $y_{i}$ to $y_{j}$ (conditional on $Y_{k}$) quantifies the extent to which $y_{i,n}^{-}$ improves the prediction of $y_{j,n}$ above and beyond the extent to which $y_{j,n}$ is predicted by $y_{j,n}^{-}$ and $Y_{k,n}^{-}$. This definition is assessed in the time domain performing a regression of the present of the target on the past of all processes, yielding the prediction error $e_{j|ijk,n}=y_{j,n}-\mathbb{E}[y_{j,n}|Y_{n}^{-}]$, and on the past of all processes except the driver, yielding the prediction error $e_{j|jk,n}=y_{j,n}-\mathbb{E}[y_{j,n}|y_{j,n}^{-},Y_{k,n}^{-}]$ ($\mathbb{E}$ is the expectation operator). The prediction error variances resulting from these ”full” and ”restricted” regressions, $\lambda_{j|ijk}=\mathbb{E}[e_{j|ijk,n}^{2}]$ and $\lambda_{j|jk}=\mathbb{E}[e_{j|jk,n}^{2}]$ are then combined to yield GC from $y_{i}$ to $y_{j}$ as Geweke (1984)

[TABLE]

The measure (1) is the log-likelihood ratio for the two linear regressions associated with the projections $\mathbb{E}[y_{j,n}|Y_{n}^{-}]$ and $\mathbb{E}[y_{j,n}|y_{j,n}^{-},Y_{k,n}^{-}]$ Barnett and Seth (2015). According to the axiomatic definition of transfer entropy and its equivalence (up to a factor 2) to Eq.(1) for Gaussian variables Barnett et al. (2009), the GC can be interpreted as the rate of ”information transfer” from driver to target.

II.2 State-Space Formulation

Following the derivations of a recent work by Barnett et al. Barnett and Seth (2015), now we move to describe the computation of GC for state state space (SS) processes. The very well known linear SS representation of an observed multivariate process $Y$ is given by Anderson and Moore (1979)

[TABLE]

where $X$ is the state (unobserved) process, and $W$ and $V$ are zero-mean white noise processes with covariances \Xi$$\equiv$$\mathbb{E}[W_{n}W_{n}^{T}] and \Psi$$\equiv$$\mathbb{E}[V_{n}V_{n}^{T}], and cross-covariance \Theta$$\equiv$$\mathbb{E}[W_{n}V_{n}^{T}].

The SS process has an equivalent representation, referred to as ”innovations form” SS (ISS), evidencing the innovations $E_{n}=Y_{n}-\mathbb{E}[Y_{n}|Y_{n}^{-}]$, i.e., the residuals of the linear regression of $Y_{n}$ on its infinite past $Y_{n}^{-}$, whose covariance matrix is \Phi$$\equiv$$\mathbb{E}[E_{n}E_{n}^{T}]. The ISS representation, which is typically associated with Kalman filtering, is characterized by the state process $Z_{n}=\mathbb{E}[X_{n}|Y_{n}^{-}]$ and by the Kalman Gain matrix $\mathbf{K}$:

[TABLE]

The SS and ISS representations share the state and observation matrices $\mathbf{A}$ and $\mathbf{C}$, and differ in the noise matrices ($\mathbf{\Xi},\mathbf{\Psi},\mathbf{\Theta}$) and ($\mathbf{K},\mathbf{\Phi}$). To find the ISS parameters ($\mathbf{A},\mathbf{C},\mathbf{K},\mathbf{\Phi}$) from the SS parameters ($\mathbf{A},\mathbf{C},\mathbf{\Xi},\mathbf{\Psi},\mathbf{\Theta}$) it is necessary to solve a so-called discrete algebraic Ricatti equation (DARE), formulated in terms of the state error variance matrix $\mathbf{P}$:

[TABLE]

from which $\mathbf{K}$ and $\mathbf{\Phi}$ are obtained as

[TABLE]

Then, GC can be computed from the ISS parameters as follows Barnett and Seth (2015). The error variance of the full regression is simply the $j-th$ diagonal element of the innovation covariance, $\lambda_{j|ijk}=\mathbf{\Phi}(j,j)$. The error of the restricted regression is obtained by forming a $submodel$ that excludes the driver process, i.e. considering a state space model with state equation (3a) and observation equation

[TABLE]

where the superscript (a) denotes selection of the rows with indices $a$ of a matrix. Of note, this technique represents one of the key elements in the derivation of a rigorous formalism for defining the information flow between the states of both discrete stochastic mappings and continuous time stochastic systems Liang (2016). Here, we have that the submodel (3a, 6) is an SS model with parameters ($\mathbf{A},\mathbf{C}^{(jk)},\mathbf{K}\mathbf{\Phi}\mathbf{K}^{T},\mathbf{\Phi}(jk,jk),\mathbf{K}\mathbf{\Phi}(:,jk)$), which can be converted to an ISS model with innovation covariance $\mathbf{\Phi}^{R}$ solving the DARE (4,5), so that the restricted error variance becomes $\lambda_{j|jk}=\mathbf{\Phi}^{R}(j,j)$. This shows that GC can be computed numerically from the ISS parameters ($\mathbf{A},\mathbf{C},\mathbf{K},\mathbf{\Phi}$) of an observed process $Y$.

III Multiscale Granger Causality

In this Section we develop our framework for the multiscale computation of GC for linear multivariate processes. Here we consider the most common operalization of GC, i.e. that grounded on the AR representation of multivariate processes Granger (1969); Barnett et al. (2009); Bressler and Seth (2011); Porta and Faes (2015):

[TABLE]

where $p$ is the model order, $A_{k}$ are $M\times M$ matrices of coefficients, and $U_{n}=[u_{1,n}\cdots u_{M,n}]^{T}$ is a vector of $M$ zero mean Gaussian innovation processes with covariance matrix \Sigma$$\equiv$$\mathbb{E}[U_{n}U_{n}^{T}]. To study the observed process $Y$ at the temporal scale identified by the scale factor $\tau$, we apply to each constituent process $y_{m},m=1,\ldots,M$,the following transformation which performs a weighted average of $q$ consecutive samples of the process:

[TABLE]

This rescaling operation corresponds to transform the original process $Y$ through a two-step procedure that consists of the following filtering and downsampling steps, performed respectively with a filter of order $q$ and a rate of downsampling equal to $\tau$. The filtering and downsampling steps yield respectively the processes $\tilde{Y}$ and $\bar{Y}$ defined as:

[TABLE]

The change of scale in (8) generalizes the averaging procedure originally proposed in Costa et al. (2002), which sets $q=\tau-1$ and $b_{l}=1/\tau$. In this study we identify the $b_{l}$ as the coefficients of a linear lowpass filter with cutoff frequency set at $f_{\tau}=1/2\tau$ to avoid aliasing in the subsequent downsampling step. Here, we develop a FIR filter of order $q$ using the window method and implementing a Hamming window Oppenheim and Schafer (1975). The design of a causal filter which performs one-side filtering was chosen on purpose to avoid that past and future samples mix up in the filtering process with potentially harmful consequences on the evaluation of causality. The use of one-side filtering is in agreement with the adoption of a Euler forward scheme, rather than a central differencing scheme, in the definition of causal measures of information flow (see, e.g., Liang (2014, 2015)). Moreover, the use of a FIR filter achieves better elimination of the fast temporal scales with respect to the averaging procedure commonly adopted in multiscale complexity analysis Costa et al. (2002); in our preliminary work Faes et al. (2016) we have shown indeed that the use of simple averaging may induce the detection of spurious causal influences over uncoupled directions.

Substituting (7) in (9a), the filtering step leads to the process representation:

[TABLE]

where $\mathbf{B}_{l}=b_{l}\mathbf{I}_{M}$ ($\mathbf{I}_{M}$ is the $M\times M$ identity matrix). Hence, the change of scale introduces a moving average (MA) component of order $q$ in the original AR$(p)$ process, transforming it into an ARMA$(p,q)$ process. Then, exploiting the close relation between ARMA and SS models Aoki and Havenner (1991), the process (10) is turned into an ISS model by defining the state process $\tilde{Z}_{n}=[Y_{n-1}^{T}\cdots Y_{n-p}^{T}U_{n-1}^{T}\cdots U_{n-q}^{T}]^{T}$ that, together with $\tilde{Y}_{n}$, obeys the state equations (3) with parameters ($\tilde{\mathbf{A}},\tilde{\mathbf{C}},\tilde{\mathbf{K}},\tilde{\mathbf{\Phi}}$), where

[TABLE]

[TABLE]

[TABLE]

and where $\tilde{\mathbf{\Phi}}=\mathbf{B}_{0}\mathbf{\Sigma}\mathbf{B}_{0}^{T}$ is the covariance of the innovations $\tilde{E}_{n}=\mathbf{B}_{0}U_{n}$. Moreover, the downsampled process $\bar{Y}_{n}$ can be put in ISS form directly from the ISS formulation of the filtered process $\tilde{Y}_{n}$: exploiting a recent result (theorem III in Solo (2016)), we find that $\bar{Y}_{n}=\tilde{Y}_{n\tau}$ has an ISS representation with state process $\bar{Z}_{n}=\tilde{Z}_{n\tau}$, innovation process $\bar{E}_{n}=\tilde{E}_{n\tau}$, and parameters ($\bar{\mathbf{A}},\bar{\mathbf{C}},\bar{\mathbf{K}},\bar{\mathbf{\Phi}}$), where $\bar{\mathbf{A}}=\tilde{\mathbf{A}}^{\tau}$, $\bar{\mathbf{C}}=\tilde{\mathbf{C}}$, and where $\bar{\mathbf{K}}$ and $\bar{\mathbf{\Phi}}$ are obtained solving the DARE (4,5) for the SS model ($\bar{\mathbf{A}},\bar{\mathbf{C}},\mathbf{\Xi}_{\tau},\tilde{\mathbf{\Phi}},\mathbf{\Theta}_{\tau}$) with

[TABLE]

The overall procedure for multiscale analysis is depicted in Fig. 1: filtering with cutoff $f_{\tau}$ the AR($p$) process $Y$ yields an ARMA($p,q$) process, which is equivalent to an ISS process; the subsequent downsampling yields a different SS process, which in turn can be converted to the ISS form solving the DARE. Thus, both filtered and downsampled processes are described by ISS models, whose parameters can be used to compute GC by forming a submodel in which the target is observed without considering the driver process (eq. (6)) and solving the DARE for this submodel. This procedure allows analytical computation of GC measures for multiscale (filtered and downsampled) processes, which is illustrated in the following for simulated and real time series.

IV Simulation Study

IV.1 VAR process with time delayed causal interactions

Theoretical analysis and simulations are first performed for the bivariate AR process with equations:

[TABLE]

with iid noise processes $u_{1,n},u_{2,n}\sim\mathcal{N}(0,1)$. The parameters in (12) are set to generate autonomous dynamics with strength $c_{ii}$ and lag $d_{ii}$ for each scalar process $y_{i}$, and causal interactions with strength $c_{ij}$ and lag $d_{ij}$ from $y_{j}$ to $y_{i}$ ($i,j=1,2$). We consider two parameter configurations: unidirectional interaction at lag 2 from $y_{1}$ to $y_{2}$, obtained setting $c_{12}=0$ and $c_{21}=0.5,d_{21}=2$, where also autonomous dynamics are generated for $y_{1}$ but not for $y_{2}$ ($c_{11}=0.5,d_{11}=1,c_{22}=0$); bidirectional interactions with different lags and strengths ($c_{12}=0.75,d_{12}=2;c_{21}=0.5,d_{21}=7$) in the presence of autonomous dynamics for both processes ($c_{11}=c_{22}=0.5,d_{11}=d_{22}=1$). First, we study the exact values of multiscale GC obtained from the true AR parameters. The theoretical trends depicted in Figs. 2 and  3 (black solid lines) document from the perspective of SS modeling the invariance of GC under filtering, already proven in Barnett and Seth (2011). The behavior of the information transfer across multiple temporal scales is thus shaped by the downsampling step, revealing the tendency of GC to peak at scales corresponding with the lag of the imposed causal interactions: maximal information transfer is found at $\tau=2$ for $F_{1\rightarrow 2}$ in the unidirectional scheme (Fig. 3a,b), and at $\tau=7$ for $F_{1\rightarrow 2}$ and $\tau=2$ for $F_{2\rightarrow 1}$ in the bidirectional scheme (Fig. 3c). The behavior is general, in the sense that it was observed also for different parameter configurations.

Next, we test reliability of multiscale GC estimates obtained from finite length realizations of (12) (a) following a naïve approach whereby GC is computed performing full and restricted regressions on the filtered and downsampled time series, and (b) performing AR identification on the original time series and then applying the new proposed framework to the estimated AR parameters [inthisstudyARmodelswereidentifiedthroughordinaryleastsquaresandusingtheBayesianInformationCriterion(BIC)tosetthemodelorderasseene.g.in]marple1987digital. Application of the two approaches to 100 process realizations of 500 points is depicted in Figs. 2a, 3a and in Figs. 2b, 3b respectively, and indicates the need of state-space analysis: while the computation of GC after filtering and downsampling returns strongly biased and highly variable estimates, the new framework yields accurate detection of GC across multiple scales.

IV.2 VAR process with multiscale causal interactions

As a second simulation example, we consider the case of coupled stochastic processes displaying multiscale structure and scale-dependent causal interactions. Specifically, we consider the bivariate process $Y_{n}=[y_{1,n}y_{2,n}]^{T}$ obtained as the instantaneous mixing of pairs of scalar processes taken from the two bivariate AR processes $X_{n}=[x_{1,n}x_{2,n}]^{T}$ and $Z_{n}=[z_{1,n}z_{2,n}]^{T}$ :

[TABLE]

In the bivariate processes $X$ and $Z$, autonomous dynamics are set for the subprocesses $x_{1}$ and $z_{1}$ according to Eqs. (13a, 13c); the autonomous rhythms are obtained by placing a pole with modulus $\rho_{x_{1}}=\rho_{z_{1}}=0.95$ in the complex plane representation of each individual subprocess, and the phase of the poles is varied to obtain slow oscillations for $x_{1}$ ($\phi_{x_{1}}=0$) and faster oscillations for $z_{1}$ ($\phi_{z_{1}}=0.47$). Moreover, causal interactions are imposed from $x_{1}$ to $x_{2}$ and from $z_{1}$ to $z_{2}$ according to Eqs. (13b, 13d); in this study, the variances of the uncorrelated innovations are set to $\lambda_{u_{1}}=0.25,\lambda_{u_{2}}=0.5,$ for the process $U$, and to $\lambda_{w_{1}}=1,\lambda_{w_{2}}=0.5,$ for the process $W$. Then, the mixing obtained with Eqs. (13e, 13f) is such that the bivariate process $Y$ exhibits causal interactions from $y_{2}$ to $y_{1}$ visible at small time scales for the faster oscillations, as well as causal interactions from $y_{1}$ to $y_{2}$ visible at larger time scales for the slower oscillations.

To test the ability of our framework to detect these multiscale behaviors, we performed VAR identification on realizations of 1000 data points of the observed process $Y$ and then computed the GC for temporal scales ranging from 1 to 15; the model order was optimized using the BIC criterion, and the order of the lowpass FIR filter was set to $q=6$. The results of the analysis are reported in Fig. 4. The exemplary realizations shown in Fig. 4(a) display multiscale patterns characterized by a slow rhythm (cycle of $\sim 70$ points) superimposed to faster oscillations (cycle of $\sim 13$ points); the multiscale GC analysis reveals that the direction of interaction is from $y_{2}$ to $y_{1}$ for the fast oscillations ($F_{y_{2}\to y_{1}}$ is maximum at low time scales), and from $y_{1}$ to $y_{2}$ for the slower rhythm ($F_{y_{1}\to y_{2}}$ emerges at higher time scales when fast oscillations are filtered out). These results are confirmed by the analysis extended to several process realizations reported in Fig. 4(b), which indicates that GC peaks at scale $\tau=2$ along the direction $y_{2}\to y_{1}$, and at scales $\tau=4$ and $\tau=8$ along the direction $y_{1}\to y_{2}$, thus detecting the multiscale patterns of bidirectional interaction imposed in the simulation.

V Practical Application

As a practical application, we consider the multiscale analysis of GC between carbon dioxide concentration ($CO_{2}$) and global temperature ($GT$). It is widely considered that the raise of $CO_{2}$ is a main cause of global warming Booth et al. (2012), although the validity of such causal relation is still under debate. The problem of understanding the causes of climate change is usually tackled by numerical experiments using Global Climate Models Delworth et al. (2006) which aim at catching the complexity of climate dynamics. However, data-driven approaches, as GC, are also fruitful in assessing cause-effect relationships between temperature and external forcings. In Kodra et al. (2011) it has been shown that $CO_{2}$ Granger causes temperature, based on data from 1860 to 2008, partly from ice cores, and analyzing second differences of both $CO_{2}$ and $GT$. Similar conclusions were found in Attanasio (2012), using GC, in Stips et al. (2016) by estimating the time rate of information flowing from one time series to the other, and in Smirnov and Mokhov (2009) using a physical approach.

Here, we first analyze the global land-ocean temperature index GIS and $CO_{2}$ concentration CO (2) measured at monthly resolution from March 1958 to February 2017. The measured time series, lasting 708 data points, are shown in Fig. 5(a). To fulfill stationarity criteria, we de-trended the two series applying an L1 norm filter. The analyzed time series, normalized to zero mean and unit variance, are shown in Fig. 5(b). We applied the proposed framework to compute GC along the two directions of interaction for time scales ranging from 1 to 100 years. To test the statistical significance of the estimated multiscale patterns of causality, the analysis was performed both for the original time series and for a set of 100 pairs of uncoupled time series sharing the autocorrelation and amplitude distribution of the original series; these surrogate series are generated using the iterative amplitude-adjusted Fourier transform (IAAFT) algorithm Schreiber and Schmitz (1996).

The results depicted in Fig. 5(c) show that the GC along both directions is not distinguishable from the corresponding surrogate counterparts at $\tau$ equal one year, i.e. when standard GC analysis not encompassing multiple time scales is performed. On the other hand, the multiscale approach reveals, at longer time scales ($>$ 10 years), that GC is significantly higher than the surrogate threshold along both directions of interaction, thus showing the need of a multiscale approach to put in evidence this mutual interdependency between $CO_{2}$ and $GT$. Moreover we remark that the GC along the direction $CO_{2}\to GT$ is characterized by a higher value of the statistics computed on the original time series, but also by higher values for the surrogate time series, compared with the direction $GT\to CO_{2}$.

Next, we dramatically change the time scales and turn to consider paleoclimatological data, so as to analyze the GC between $GT$ and $CO_{2}$ concentration on the Vostok Ice Core data from 400,000 to 6,000 years ago, extended by the EPICA Dome C data which go back to 800,000 years ago ICE . The two time series, which are sampled with non-uniform time spacing, are shown in Fig. 6(a). Here, we studied the data resampled to an uniform time spacing of 729.77 years, corresponding to a time series length of 1095 points, and after normalization to zero mean and unit variance (Fig. 6(b)); results of multiscale GC analysis did not change substantially if the original non-uniformly sampled time series were considered, or applying slightly different uniform resampling. In Kang and Larsson (2014), empirical evidences for the existence of Granger causal influences along both directions $CO_{2}\to GT$ and $GT\to CO_{2}$ have been found after correcting for deterministic trends on the same data. Here, applying our framework for multiscale causality analysis we obtained the GC curves reported in Fig. 6(c). We find that, at paleolithic time scales, the GC $GT\to CO_{2}$ is highly significant and peaks around 1000 and 10000 years. This result may be related to the lags between Antarctic deglacial warming and $CO_{2}$ increase reported in Caillon (2003), and also confirms the good evidences reported on the fact that higher global temperatures do promote a rise of greenhouse gas levels Scheffer et al. (2006). The opposite causal influence from $CO_{2}$ to GT is much less pronounced and exceeds the IAAFT threshold for statistical significance only at very small time scales.

Summarizing, our results show that carbon dioxide and temperature changes are interdependent at multiple time scales, with a predominance of $GT\to CO_{2}$ effects at paleolithic scales, and the presence of bidirectional causal interactions between $GT$ and $CO_{2}$ at the time scales of modern climate. These results support the expectations that changing temperatures could be held responsible for changes in greenhouse gas concentrations on paleolithic time scales, while during the last 60 years the effect of human activities becomes evident as anthropogenic radiactive forcings are seemingly driving the global temperature changes. These causal relationships between $CO_{2}$ and global warming have been recently demonstrated in Stips et al. (2016): in that work, the use of the rigorous formulation of information flow provided by Liang (2014) led to evidence a clear unidirectional nature for the causal relation $GT\to CO_{2}$ in paleoclimatological data, and for the causal relation $CO_{2}\to GT$ in modern climate data; in the same work, the application of the standard GC index to modern climate data suggested the presence of bidirectional effects $CO_{2}\to GT$ and $GT\to CO_{2}$, thus pointing to some ambiguity in the assessment of a predominant direction of interaction using GC. Our results agree with this interpretation, as we do not find a prevalent causal direction using the classical GC index computed at the smallest time scale, and the use of surrogate data indicates the lack of statistical significance (Fig. 5). Nevertheless, the analysis performed at higher time scales reveals the existence of significant GC $CO_{2}\to GT$ and, for the first time to our knowledge, a nontrivial GC $GT\to CO_{2}$ also in modern climate. Although this result needs to be confirmed by the implementation of more robust measures of information flow, it may be of great relevance for climate studies as it is indicative of a positive feedback which will increase the effect of anthropogenic emissions on global temperatures.

VI Conclusions

The present study makes the first step toward the theoretical understanding of multiscale causal relations between coupled stochastic processes, and opens the way to the reliable estimation of these relations starting from simple AR identification. This will likely boost new impetus for research in the area of data-driven causality analysis, both in physics and in a wide variety of applicative fields. The proposed framework is flexible enough to encompass more general model representations that may unveil important multiscale features of coupled processes. For instance, integrating the standard AR representation with fractional integrated (FI) innovation modeling Sela and Hurvich (2009) would be straightforward as ARFI models have an SS representation, and would easily lead to assess multiscale GC in the presence of long-range correlations.

The proposed setting provides also the basis to expand the applicability of multiscale GC to nonstationary and nonlinear SS processes Kitagawa (1987), and to formalize exact computation of cross-scale information transfer within and between multivariate processes Paluš (2014), thus opening new avenues of research in the evaluation of causal interactions among coupled processes. Of particular interest in this context is the recent formalization of the notion of information flow based on first principles, rather than axiomatic postulates or empirical proposals, implemented in Liang (2016). The latter work completes the rigorous formalism introduced in Liang and Kleeman (2007a, b) and provides a well-principled alternative to the operational implementation of GC, and of transfer entropy intended as its non-parametric generalization, which are known to be complicated in many practical settings to an extent that spurious causalities may be revealed (e.g., in the presence of unobserved variables, measurement noise, or inappropriate time resolution) Smirnov (2013); Hahs and Pethel (2011); Nalatore et al. (2007); Chicharro and Panzeri (2014). Hence, the availability of a rigorous derivation of the information flowing among the components of discrete time stochastic mapping, provided in Liang (2016) and extended therein to continuous time stochastic mappings and to deterministic systems, certainly constitutes a firm basis for the design of a more faithful analysis of causality between dynamical system components operated across multiple temporal scales.