Causal discovery or identification has been one of the few most fundamental goals in machine learning and economics (e.g., Eichler, 2013). It is preliminary for causal inference. In the discovery process, in general, an automatic algorithm iteratively identifies the latent causations among the data variables. However, the latency could hide a deep and complicated data structure over which we merely observe a superficially precedent connection between any two variables. The pairwise precedence relations indeed constitutes a directional network with various strength and uncertainty for each edge. The bilateral relations, however, are not in a logical and consistent order. This paper applies a network equilibrium (Hu and Shapley, 2003) and defines a new type of causality. The solution

The views expressed herein are those of the author and should not be attributed to the IMF, its Executive Board, or its management. integrates direct and indirect influences among the data, accrues the real causal relations, and filters out the noisy ones. Defining causality is a profound convoluted question with many possible answers which do not satisfy everyone (Granger, 1980). There are dozens of approaches in the literature. They generally capitalize on the conditional probability under interventions or the significance of regression coefficients (e.g., Granger, 1969 and1980;Pearl, 1995). However, conditional probability involves a parametric probability distribution and the precision in estimating them from the data. In a regression, the results depend entirely on the appropriate selection of variables, and causal factors that are not incorporated into the regression model cannot be represented in the output. Adding or dropping a covariate, for example, could make a significant effect insignificant; it could also make an insignificant one significant. Moreover, the estimated causal relations are not transitive. For example, if x significantly Granger causes y and y significantly Granger causes z; however, x may not Granger cause z. In this causal chain, y could merely be a pass-through. Direct but artificial ones often marginalize indirect causal factors. In studying climate change, for example, solar activities have had little trending vibration in the past hundred years and seem to have no impact on the up-trending global warming. The aggregate impact over a thousand years, however, may be no longer negligible. Also, indirect impacts passing through other determinants to the change have played an essential role in the last hundred years. A regression model, however, only averages the direct effect in the estimate sample, and does not automatically add up the small effects. In this paper, we propose a new approach that attempts to resolve the above issues. Our research is based on two presumptions. First, a causal effect could act on sequential occasions, which could occur over a long time and happen multiple times. Thus, an aggregate and accrual effect would be suitable for measurement. Secondly, indirect causal reason passes through one or more intermediaries and gains causal strength through the spillovers from these intermediaries. The device we use is a sequence of fictitious vector autoregressive (VAR) models in which the causal relations are under investigation. The shocks to the VAR are natural direct interventions; responses to the shocks are plausible effects. The limit variance decomposition aggregates both direct and indirect, both short-run and long-run effects. The decomposition quantifies the bilateral impacts among all variables. It defines a network with directed edges over which the bilateral impacts are the flow volumes. Then we define an equilibrium over the network, which resolves potential conflicts in the bilateral influences or plausible causations. In the literature, Leamer (1985) suggested that variance decomposition may contain causal interpretation and stated that Sims (1980) had already implied that. The equilibrium has a global and a local version. On a global scale, it produces a unique solution when the network is highly connected. The vector solution quantifies the causality distribution for the movement of the time series variables. When there are both exogenous and endogenous variables, the values in the solution for any endogenous one would be statistically insignificant. A hypothesis test can formally conduct this. Similarly, we also test the disconnection among the variables after element-wisely removing these endogenous variables from the fictitious VAR model. The local version of causality distribution attributes each determinant a fair share of contribution to an effect variable. The precision of these hypotheses tests is evidenced in our simulation studies. We also apply the distribution to study climate changes and find which factors are more important than others. The rest of the paper is organized as follows: Section 2 introduces a bilateral interaction network among a vector of time series data by a vector autoregressive model. Section 3 derives the valuation formula, which measures the amount of causality each time series contains; the solution is the global causality distribution which applies to the whole vector. Section 4 studies a causality distribution that measures causal variables' contributions to any effect variable. Section 5 provides algorithms to identify the causal relations. Section 6 uses the simulation studies to test the effectiveness of the new approach and conducts an empirical study for climate change. Finally, Section 7 concludes with further comments. Our exposition is self-contained, and the proofs are in the Appendix. We study a group of covariance stationary time series, y 1 , y 2 , ..., and y n (i.e., their means, variances, and covariances remain unchanged over time), which can be adequately described by a linear model with lagged values of themselves. We index the time series by N = {1, 2, • • • , n}. The data generating process, however, is a black box; we attempt to identify which are endogenous, which are exogenous, and which cause which variables. For convenience, we stack the time

The bilateral interventions are directional, and all together make a directed network. Figure 1 is an example in which "→" is for one-way directional influence, and "↔" is for two-way directional influence in other figures. In the network, nodes are the variables, and the directional edges are the intervention relations.
We study a few basic relations which may have similarities to the concepts in the theory of discrete-time Markov chains. For example, an exogeneity class corresponds to an irreducible class while the endogeneity class corresponds to the transient class. Familiarity with the Markov chain theory is not required but helpful. Exogeneity Classes: There are a few levels of exogeneity (e.g., Engle et al., 1983). We use the popular version: a class of variables is exogenous to all others if by an unrestricted vector autoregression, however, the coefficients from other variables could still be significant in the equations of y 7 and y 8 due to the sampling errors. Conversely, y 7 and y 8 may also be significant in other equations in the VAR model. Thus, the plausible causal linkage may not be reliable from the estimated coefficients alone. In a hierarchy, one causal factor may not directly influence effect ones; also, all direct causations have a one-way direction. In Figure 3, there are three levels of hierarchy. {y 1 , y 2 } is the exogeneity class but has no direct causal effect to the third level. After removing {y 1 , y 2 }, both {y 3 , y 4 } and {y 5 } are de facto exogeneity classes. In a VAR model, restrictions should be placed such that some coefficients are zeros to ensure one-way directions are in place. Also, through a sequence of VARs, conditional exogeneity could be finally figured out. In a circular causal structure, there exist one or more circular causal chains. In Figure 4, there are two circular causal chains: {y 9 } → {y  Periodic: This is a special circular chain, and periodicity is also a class property. The circular causal chain has a minimal period d if the chain returns to the same variable after traveling multiple d periods. Figure 5 has a period 6. The number of intermediaries for a circular chain is the same for any member to recurrent to itself. The endogeneity class may contain one or more subclasses which are exogenous to the other variables in the endogeneity class. In Figure 6, {y 5 , y 6 , y 7 , y 8 , y 9 } is the endogeneity class. Itself has a subclass {y 5 , y 6 } which is exogenous to {y 7 , y 8 , y 9 }. To find the contribution of y 5 to y 7 , for example, we put y 1 , y 2 and y 3 as exogenous in the VAR of y 5 , • • • , y 9 , and ignores y 4 . The general unrestricted VAR (GUVAR) with p lags is of the form

where y .t , ν, and u t are n dimensional column vectors, and A i 's are n × n unknown coefficient matrices. As usual, we assume that the residual u t is serially uncorrelated and has a constant variance-covariance Σ over time t. However, the model ( 1) is fictitious and has not differentiated the exogenous and endogenous variables in the list, as well as the causal relations among them. It is a mechanism to identify them. Restrictions could be placed on the short-run coefficient matrices to reflect the dependence among the component time series. When there are both exogenous and endogenous variables in the list y .t , in theory, the endogenous ones have no short-run impact on the exogenous ones. Therefore, after some rearrangements, the true VAR(p) has the pattern for all the coefficient matrices

For two or more isolated classes of exogenous variables, the coefficient matrices of the true VAR(p) have the pattern of block diagonal as

In practice, however, we are uncertain if there are any exogenous variables. Also, an estimated coefficient matrix of (1) would have a nonzero upper right block, due to the sampling errors. The VAR(p) in (1) can be rephrased to a VAR(1) structure (4) by stacking y 1 , y 2 , • • • , y n into a large vector Y ,

where

, and

Besides, A is a np × np matrix, I n is the n × n identity matrix, Y , V and U are np dimensional column vectors. Clearly, U t is also serially uncorrelated and has variance-covariance cov(U t ) = J ΣJ where J = I n 0 . . . 0 is a n × np matrix. In general, some elements of a block of an estimated coefficient matrix may be statistically insignificant; but some may be significant. Even all elements in the block are insignificant, we can still not tell if the block is zero statistically. For the linear engine (1) of the covariance-stationary y .t , the driving forces are the residuals or shocks u t , either positive or negative. Without the forces, the series remain as a constant vector and no cause-effect relations could be detected. Though human could also engineer shocks (now called interventions), we ignore the mechanism how they are generated. With a given u t , the whole series vector responds from time t in which a tiny slice of causal effect would be hardly noticeable. The aggregated effects through a sequence of u t could be enormous. However, a linear regression generally averages, not aggregates, the effects from the innovations. We analyze the variances of the shock series where the regression amplifies or mitigates these effects. The shock series is already embedded in the data y .t through the dynamics (1). Given Y t , the h-step ahead vector from (4) is

(

The conditional variance depends on the length of h, not on the value of Y t . A s J ΣJ(A s ) , if it exists.
Then we have a steady-state variance which satisfies the Lyapunov equation

The solution to the Lyapunov equation solves the following linear equation (cf, e.g., Hamilton, 1994, Equations 10.2.13 and 10.2.18)

where ⊗ is for the Kronecker product of matrices. Therefore,

This also implies the covariance of vec(X) has an even more complicated form. Let L be a lower triangular matrix obtained by a Cholesky decomposition of

L j L j where L j be the jth column of L. Then the contribution of y j to X is

Element-wise ratios between (8) and ( 7) provide y j 's contribution in percentage to the limit variance of y .t . The covariance matrix Σ could also be written as Σ = (LQ)(LQ) for any orthogonal matrix Q. However, LQ has an alignment issue with the order in y .t . A common practice is to order the components of y .t from the most exogenous to the least exogenous (e.g., Sims, 1980). This way, the shock from the first component series of y .t affects all others contemporaneously and the L 1 vector designates the contemporary effects; the shock from the last component series only affects the series itself contemporaneously where L n has only one non-zero element. However, (8) includes not only the direct contemporaneous impact. For a finite h, we take the first n elements from Y t+h in (5) to get cov(y

where Φ s = JA s J . Thus, the variance of the h-step forecast of variable

2 where e i is the i-th column of I n . Therefore, the amount of forecast variance of variable y i , cov(y i,t+h |y t ), accounted for by the shocks to variable y j is given by ω ij,h ,

The asymptotic distribution of ( 9) could be derived (e.g., Lütkepohl, 2005, page 111); however, it does not account for the identity n j=1 ω ij,h = 1 and the nonnegativity of (9). Thus, the distribution cannot be used in the usual way for the hypothesis test of ω ij,h = 0. Theorem 1 shows the existence of ω ij,h in the infinite horizon of h. We let the unconditional covariance

which can be exactly calculated from ( 7) and ( 8). Unlike a causal inference, the discovery of causal relation should be invariant in some way. It is a YES or NO question. It does not rely on the unit the variables use or the sign and size of the coefficients in regression. A causal relation between y i and y j does not merely depend on the coefficients and their standard errors in the regression of y i and y j . It also relies on the causal structure in y .t and the global uncertainty in the modeling. The formulas (7), (8), and (9) integrate both the coefficient matrices and the uncertainty. Meanwhile, (2) and ( 3) are two common relational structures in y .t . The following two theorems indicate the invariance of patterns and linearity. Theorem 2. If all the coefficient matrices A i have the pattern (2), then Ω has the same pattern. Proof. See Appendix A2. 3), then Ω is a block lower triangle matrix. The pattern for the coefficient matrices A i restricts the short-run direct linkages. If these restrictions carry on in the long run, then the pattern preserves in Ω as it integrates both short-run and long-run. The block diagonal pattern ( 3) is a special case of (2); thus, Corollary 1 follows. However, Ω is not necessary a block diagonal; the order of the exogeneity classes matters. If class i precedes class j in the VAR, then i has a directional impact on j in Ω; but not vice versa. As an example of the linearity invariance, the emission of carbon dioxide is a causal factor to global warming, measured either in the Celsius scale or the Fahrenheit scale. The invariance also makes it possible to compare the causality contributions across the variables in y .t . When the variables instantaneously affects each other, we could write a general structural VAR (SVAR) form as

Then

and

Taking variance on both sides of ( 12) and letting t → ∞, we obtain the Lyapunov equation for X,

and its solution is

where Ã = (I n -A 0 ) -1 JA . The proportion contributed by y j is

The element-wise ratios between ( 14) and ( 13) normalizes the contributions. Theorem 4. If y .t is covariance stationary and ||(I n -A 0 ) -1 || ≤ 1, then both ( 14) and ( 13) exist and do not depend on any value of y .t . Proof. See Appendix A4. The bilateral influence matrix Ω from the general unrestricted VAR ( 1) is a Markov transition matrix, upon which one may define irreducibility, ergodicity, and periodicity. However, each element of the matrix contains a heterogeneous uncertainty and, in an empirical analysis, no element in the matrix would be zero due to the sampling errors. Of course, many elements would be statistically insignificant when the uncertainty is concerned. Ignoring the uncertainty, therefore, all variables in y .t are artificially irreducible, ergodic, and aperiodic. There are many inconsistencies in Ω if the transition has a causal interpretation. The transitivity does not hold in the variance decomposition of Ω. Secondly, indirect causal relations are not properly linked. Also, spurious causal linkages are not removed. Lastly, there exists heterogeneous uncertainty in Ω. To overcome the inconsistencies, we consider the t-step influence Ω t where Ω t (i, j) is the y j 's indirect influence over y i , passing through t -1 intermediaries including both y i and y j . The t-step matrix make indirect influences direct. In Figure 3, the first level has a direct influence to the third level in Ω 2 .
The matrix also integrates any intermediate effects as

Sun's causal effects on global warming, for example, could accumulate for many years and repeatedly through many intermediates (e.g., oceanic evaporation and ozone). Finally, for an empirical Ω, the indirect influences eventually become stable over time, lim

where 1 n is the n × 1 vector with all ones, π i ≥ 0, and

Therefore, the average direct and indirect influence also have the same limit,

In the limits, the ith column remains the same across all rows. It is y i 's causality to the movement of the whole vector y .t . So we call π i to be y i 's global causality. The row vector π = (π 1 , π 2 , • • • , π k ) also satisfies the following counterbalance equilibrium (Hu and Shapley, 2003) over the directional network Ω:

In the counterbalance, on the one hand, y i 's causality π i to the movements of y .t derives from its effect on all components of y .t . On the other hand, it also distributes its π i to others which cause y i 's changes. In the distribution of π i , the fair division is through the ith row of Ω, which is the percentages of y i 's variability explained by each individual time series. So y j gets π i ω ij from π i . Similarly, π j also distributes π j to others and y i gets π j ω ji . Looping on all j in N , we get the following equation (16) and (15). So, the ith column of Ω stipulates the weights with which y i collects its causality distribution. The causality distribution π is a mixed cooperative and non-cooperative solution. It is noncooperative because i∈N π i = 1; one series' gain means another's loss. It is also cooperative because (16) implies that a larger π j improves π i more, as long as ω ji > 0 remains the same. Also by (16), increasing ω ji also enhances π i . A long-run of competitions and cooperations finally reaches the ultimate solution π. The distribution π could also come from an endogenously weighting. Ignoring the spillover effects, the ith row of Ω is a proper assignment of responsibility or reason for y i 's movement alone; but it is not a proper measurement for the whole vector of y .t . We assign a number π i to weigh the ith row where important variable gets more weight and less important one gets less weight. Then the weighted assignment of responsibility is πΩ, which also measures each variable's contribution to the movement of the vector y .t . So π is proportional to πΩ. Setting n i=1 π i = 1, we get the equation (15). Another application of the counterbalance equilibrium can be found in a sorting algorithm for big data (Hu, 2020). Instead of using Ω in the equilibrium, one could also apply the frequency decompositions (Geweke 1982; Kaminski et al. 2001) or the generalized variance decomposition (Lanne and Nyberg, 2016). Both these two decompositions also have a unit sum for each row. An iterative algorithm for (15) is

for a given nonnegative vector π (0) such that it has a unit sum. The limit of π (t)  would be π for an empirical Ω. The general unrestricted VAR in (1) is highly artificial as it ignores the causeeffect relations within the vector y .t . Consequently, this may result in overwhelmingly many insignificant elements in Ω. One solution is to iteratively identify the causal relations and then re-specify the VAR model using already identified relations. For example, hypothesis tests could be used to make the VAR from general to specific and from unrestricted to restricted. By speculating a possible relation after observing some signals, one generally sets up a hypothesis on the relation. Finally, he or she extracts the evidence from the data to approve or reject the hypothesis. The distribution π mitigates the uncertainty in Ω by its spillover effects from the third parties. Indeed, the spillover effects take two further actions: confirmation of real comparative causality and off-setting of noisy ones. Consequently, π i could still be insignificant even if many elements of Ω's ith column are statistically significant. If y i has a real comparative causality in the directional network, then Ω has consistent large values in the i's column, especially for variables with large values in π. If the advantage is a noise, on the other hand, then the sampling errors could largely account for the positive entires in the ith column. To see the effects of a small variation of Ω on π, let us first shock both Ω(j, i) and Ω(k, i) by δ. This can be seen as y i has some real improvement. Secondly, we shock Ω(j, i) by δ and Ω(k, i) by -δ. This can be a case of pure white noise. By Theorem 5, the effects add up when two shocks have the same sign while they offset each other when the shocks are opposite. The small-sized derivative in Theorem 5 also indicates that π i 's standard error is small compared to those for the elements in Ω's ith column. In the theorem, we use three notations for simplicity: the matrix Z i is the transpose of Ω with the ith row and the ith column removed; the column vector α ji takes the jth row from Ω and then drops its ith element; and π -i is the row vector when π i is removed from π. Theorem 5. For any i, j ∈ N ,

Proof. See Appendix A5. The filtration could also be explained by the chains of direct influences. The unit sum in each row of Ω places a zero-sum restriction on the noise of each row. So some noise are positive while others are negative.
For a large t,

which is a sum of n t-1 tiny particles. Each particle is the product of t elements from Ω in which some contain positive noise and others contain negative noise. The summation mitigates the noise when the noise offsets each other. Adding up the small pieces could make a difference; averaging them, as in a regression, is still small. While mutually causal relations exist in a set of data, the one-way cause-effect linkage is also pervasive. Say, y 1 has a direct influence on y 2 . But this does not imply that y 2 is necessarily essential to y 1 's variability. Generally, we may put all direct or indirect influencing and influenced data into one huge set and analyze the influences between the classifications inside the huge set. In Figure 7, for example, classes 1 and 2 influence class 3; however, the latter has no impact on classes 1 and 2. One question is how much is the relative causal importance of classes 1 and 2 to class 3. That is, which class is more important to class 3? And by how much? For example, when applied to climate changes, there are many causal factors accounting for the changes. How much is attributed to carbon dioxide emissions? How about the changes in the Sun, the emissions from volcanoes, and the variations in Earth's orbit? Fig. 7: Causal-transient class. For any y i in T , some y j has causal effect on y i but not vice versa. So y i 's influence on the whole set of times series is transient and temporary, and it dies out gradually. On the other hand, an exogeneity class accrues causality from the causal-transient y i . Let T e be the restriction of Ω to the endogeneity class T (see Figure 8). Without loss of generality, say, T has m members, namely,

So T e is the sub-matrix of the last m columns and the last m rows of Ω. As usual, let 1 n be the n × 1 vector of ones, and 0 n the n × 1 zero vector. Also, let I n be the n × n identity matrix. T t e = 0 m 0 m , the m × m zero matrix. Proof. See Appendix A6. Corollary 2. I m -T e is invertible; also, if T = ∅, then y .t contains at least one exogeneity class; and finally, π i = 0 for all i ∈ T . Proof. See Appendix A7. By the above theorem, the existence of T does not affect the causality distribution π cs . So we may remove T from y before we calculate the distribution π cs , and the removal reduces the computational cost. This reduction is extremely important when y .t has a large number of causal-transient members. For Ω in . For C 1 , the causality distribution is π c1 = ( 1 2 , 1 2 ); and for C 2 , π c2 = ( 4 19 , 15  19 ). However, when there are two or more exogeneity classes in y .t , the solution to the counterbalance equilibrium (15) is not unique. One remedy is to assign some specific causality quota to each class. and we assign C s with a causality quota q s where 0 ≤ q s ≤ 1 and k s=1 q s = 1, then there exists a unique causality distribution to solve the equilibrium π = πΩ and the quota restriction j∈Cs π j = q s for all 1 ≤ s ≤ k.

Proof. See Appendix A8. For the above example of Fig. (10), if we assign the quota q 1 = .4 and q 2 = .6, then π satisfies that

The unique solution is π 1 = π 2 = 1 5 , π 3 = 12 95 , and π 4 = 9 19 . They are actually [.4π c1 , .6π c2 ]. For any i ∈ T , there may be some j ∈ T such that j directly influences i, i.e., ω ij > 0. But since j has no global causality in y .t , its influence on i must come from members in N \ T . In some sense, j acts as messenger only, passing the effects from N \ T to i. This shows that ω ij is not enough to characterize the causation of j on member i. To find the causal strength on an element in T , let C be an exogeneity class and denote by µ c i the causality by C on the member i ∈ T . Conditional on the first move in Ω, µ c i satisfies the following the steady-state equilibrium of inhomogeneous equations

Here j∈C ω ij is the causality absorbed directly by C from i while the causality absorbed indirectly by C from i is j∈T ω ij µ c j . Thus, ( 18) is the counterbalance equation with external influences; the players are the elements in T , and the external influences come from C.

We can approximately approach a solution to (18) by the iteration of

with the starting value µ (0) = 0 m . Theorem 8 establishes the existence and uniqueness of the class's causality on individuals in T ; its corollary provides a computational method for the solution. Theorem 8. There exists a unique solution µ c i to (18) for any exogeneity class C. And the solution is nonnegative. Proof. See Appendix A9. (19) and

Proof. See Appendix A10. As the transient members has no causality share in y .t , the exogneity classes undoubtedly have total ultimate influence over these endogeneity members. Theorem 9 confirms this. To illustrate the class causality on transient members in Figure 10, by 20,

.
Thus,

As indicated above, µ c i is the total causality of the exogeneity class C over the causal-transient member i ∈ T . Now how to re-distribute µ c i among the members inside C? By Theorem 10, this redistribution is through the global causality distribution within the class. Across all exogeneity classes, the total shares of distribution add up to 1, claimed in Corollary 5. Therefore, the shares of

defines a local causality distribution for i ∈ T . A large share contributes more to the movement of y i than a small share, and the distribution varies with i ∈ T . It reduces to the global causality distribution when there is only one exogeneity class. This is similar to the responsibility index defined by Shapley (1994), which measures the ultimate controlling power no matter how long the chain of commands is. In the above example, the local causality distribution is

Many relations could be set by specifications rather than by estimations. As we all know, for example, global warming on the surface of the Earth should have no causal effect to the sunspots. Also, the rooster's crow has no causal effect to the sunrise, even though they are highly correlated. In this section, we test the causal data structure in y .t by bootstrapping the data. We order the time series in the VARs using previously calculated π values, from the largest to the least. By ( 16), variable with a large π i is likely to influence more data directly. When y i is endogenous, then in theory ( 2) is the pattern for all coefficient matrices and π i is necessarily zero. Thus, the test hypothesis is

This signals an experiment design which simulates many independently identically distributed (i.i.d.) π i 's. The sample mean πi and its standard deviation σi help to decide if π i is statistically zero or not; the other statistics, such as sample median and quantiles, are useful in choosing the sample size for bootstrapped datasets. Under the null hypothesis H 0 , the z-statistic πi /σ i has an approximate standard normal distribution, by the central limit theorem. At the α% significance, therefore, we identify y i as endogenous if the test statistic z i < z 1-α , where z 1-α is the critical value. Otherwise, y i is classified as exogenous. Pick a large integer γ as the number of bootstrapped datasets; Take an random ordering of variables for the first VAR; πi ←-0 and σi ←-0; for k ← 1 to γ do Bootstrap the original data and estimate the VAR;

Calculate Ω by ( 7)-( 8) and calculate π; πi ←-πi + πi and σi ←-σi + π 2 i ; Reorder the variables according to the values in π, from the largest to the least; end λi ←-λi/γ and σi ←-σi/γ; return zi ←-λi/σi. Algorithm 1: Hypothesis Test of π i = 0 versus π i > 0. Algorithm 1 details these steps, where σi is for the estimated standard deviation of π i under H 0 . The algorithm actually calculates z i for all i = 1, 2, • • • , k. In observance of the interdependence among data in y .t , we use the algorithm to drop at most one endogenous time series at a time. The dropped one has smallest t-statistic and its t-statistic is less than the cutoff value. After it is dropped from y .t , we run the algorithm again using the updated y .t . We continue this step until no z-statistic is less than the cutoff value. After removing all endogenous time series from y .t , the rest are independent exogeneity classes. For example, y 1 and y 9 in Figure 2 are from disjoint classes. Without loss of generality, we assume that y .t contains no endogenous variables and also assume that the estimated global causality πj > πi , i.e., y j precedes y i in the VAR. By Corollary 1, when the coefficient matrices have the pattern (3), no causal shock would travel from y i to y j if they are in disjoint classes. Given π (0) = e i , therefore, π (t) j = 0 for all t > 0 in (17). This leads to the following hypothesis:

Pick a large integer γ as the number of bootstrapped datasets; Take an random ordering of variables for the first VAR; πj ←-0 and σj ←-0; for k ← 1 to γ do Bootstrap the original data and estimate the VAR;

Calculate Ω by ( 7) and ( 8) ; Calculate π by (17) given π (0) = e i ; πj ←-πj + πj and σj ←-σj + π 2 j ; Reorder the variables according to the values in π, from the largest to the least; end λj ←-λj/γ and σj ←-σj/γ; return zj ←-λj/σj. Algorithm 2: Testing (y i , y j ) from different exogeneity classes. Algorithm 2 calculates π using (17) and the starting value π (0) = e i . If the returned z j statistic is less than the cutoff value, we accept the hypothesis H 0 . When an exogenous variable y i has no causal influence on an endogenous variable y j , such as y 7 and y 6 in Figure 2, we also have π (t) j = 0 in (17) with π (0) = e i . Therefore, Algorithm 2 can test no causality from an exogeneity class to an endogenous variable. When a class C has a period d > 1, then for any i ∈ C, lim

Therefore, starting from π (0) = e i , the iteration ( 17) could be used to test the hypotheses of d = 1, by slightly modifying Algorithm 2. After finding all exogeneity classes, we further investigate the structure of the endogeneity class by adjusting the VAR model.
The new VAR ignores the exogenous variables which have no causal impact on the endogeneity class. It uses the other exogeneity classes as the exogenous variables, which have no equations in the VAR. In Figure 6, for example, the new VAR has {y 5 , y 6 , y 7 , y 8 , y 9 } as endogenous variables and {y 1 , y 2 , y 3 } as the exogamous variables. The new VAR is to test if {y 5 , y 6 } is exogenous to {y 7 , y 8 , y 9 } and if y 9 and {y 7 , y 8 } mutually affect each other. Algorithms 1 and 2 could apply to the modified VARs. This progress could go on and on until all substructures are fully identified. Specifically, for example, if y j affects y i contemporaneously, then A 0 = e i e j multiplied by an unknown coefficient in (10). We compare the VAR (1) and SVAR (10) models. If π i has a significant difference, then we say the contemporaneous causal effect is effective. Let πi be the global causality distribution for (10). Then, the test hypothesis is

Algorithm 3 implements the sample mean and the sample standard error of πi -π i under H 0 . If the z-score is within the cutoff value, then we accept H 0 and claim that the instantaneous causality is ineffective. This section applies simulation experiments to study the performance of our new estimation methods.The estimated causal structures are compared with the true ones from which data are generated. We also list the number of errors the estimation makes. In terms of exact identification of the actual models, the new approach has a precision of 83% for the causal structures in Figures 23456. For not identified cases, some causal relations are correctly identified, bringing the accuracy to 87.9%. On average, each estimated figure also includes one and half spurious causal relations. This section also conducts empirical research which estimates the determinants' contributions to climate changes. We find that the natural factors (the Sun and the Earth) have a majority share. The emission of carbon dioxide also plays a significant role. of errors. Finally, we aggregate the discrepancy statistics of all 1,000 datasets for each of the structures in Figures 2 to 6. In Table 1, the column of Exact Identification lists the cases when the estimated structure is the same as the figure; no causal relation is missed or added to the estimated structure. However, the estimated structure could miss one or more causal relations. These are listed in the column of Number of Omissions. The estimation could add one or more additional causal relations which do not exist in the figure; they are in the column of Number of Additions. The last column lists the percentage of correctly identified causal relations. For example, there are ten causal relations in Figure 2 and 91% of them are correctly identified in the estimation. Table 1 summarizes the discrepancy statistics for each of the five figures in its 1,000 simulated datasets. The average precision is around 83% when comparing the estimated figures and the figures to be estimated. About 4.5% of cases make a slight mistake that either fails to find only one true causal relation or takes a wrong causal one. Over 13% of cases make a big mistake that fails to identify at least four true causal relations or takes at least four wrong causal ones. Geological records show that there have been some significant variations in the Earth's climate over the past hundreds of years since the industrial revolution. These have been caused by many factors, either natural or human. Understanding the contributions by the causal factors helps policy-makers make the right decisions to resolve the problem. Determinants of climate change have been well studied in the literature (e.g., Crowley, 2000; Pasgaard and Strange, 2013; Stern and Kaufmann, 2014), but viewpoints are also diversified, partially due to the different interpretations of causality. We, therefore, use the commonly agreed determinants in the study. These come from three categories: the Sun, the Earth, and human activities. The categories are not exclusive; for example, carbon dioxide (CO2) could come from factories, humans' breaths, or volcanoes. Some are exogenous to the effects, such as the Sun's radiations; others are endogenous to global warmings, such as the levels of CO2 and real GDP per capita (RGDP PC). There also exist ambiguous causal relations, such as that between the ocean currents and volcanoes, for which we run the identification algorithms. Table 2 lists 15 time series for the VARs. Because of the sample size of the data, many factors are not included, such as atmospheric aerosols. We use data for the years between 1900 to 2019. The data come from GitHub.com, NASA (National Aeronautics and Space Administration), the USGS (United States Geological Survey), and the World Bank. Different time series have their own cyclic and trending patterns. The duration of the sunspot cycle, for example, is around eleven years, while El Niños and La Niñas generally occur about every two to seven years. Also, extrapolation of the estimated relation to an infinite horizon could be subject to many disastrous venerabilities. For example, the Earth's cool and warm periods cycle roughly every 100,000 years, caused by changes in Earth's orbit around the Sun. For this reason, we use (9) with h = 120 as the elements of Ω, Ω 120 for the global causality distribution, and µ (120) in Corollary 3 as the class causality. With this finite horizon, π i may not be zero for an endogenous series, which may have a nonzero local causality on another endogenous series. For the VARs, we choose 20 as the lag length, which covers the cycles of the Sun, ocean currents, and economic business cycles. However, the sample size 120 is not enough to afford 15 time series for all 20 lags.
Therefore, based on lag exclusion tests, the VARs finally pick lags 1, 2, 5, 7, 12, and 20. Table 2 summarizes the local causality distribution for global warming. Causal identification algorithms show that Solar irradiances, Earth's Milankovitch cycles and volcanic eruptions are three independent exogenous variables, and the other series are exogenous. The three exogeneity classes account for 39.6% of global warming. This percentage mostly comes from indirect effects, such as the greenhouse. Combined with the El Niño and part of CO2 and Methane, the natural factors make up more than half of the responsibility for global warming. Among the greenhouse gases that mitigate the infrared radiation from the Earth, carbon dioxide emission has about 50% share. Lastly, changes in population and vegetation cover also contribute 12%. The endogenous factors interact with each other in loops. For example, increased atmospheric temperature evaporates more water which makes the green-house thicker. Conversely, a thicker greenhouse increases the temperature, which evaporates more water. In this loop, each factor gets its causality distribution, which answers the question of Thurman and Fisher (1988). A policy-maker could capitalize on the multiplier effect in the loop. If he or she decreases the level of CO2 by 10%, for example, the direct effect is a deduction of global warming by 1.67%. Nevertheless, the spillover effects make the deduction even further. This paper provides a game-theoretic framework to study a fundamental issue in statistical learning and econometrics. When observing the bilateral influences over a group of time series, we set up a directional network by a matrix. However, we could not simply apply a causal interpretation to the matrix as transitivity is invalid, noisy causal relations are not filtered, and transient causality is spuriously embedded. By introducing a counterbalance to the matrix, the solution mitigates the noise, makes transitivity to a class property, and zeros out the impact from causal-transient data. The set of time series is decomposed into a few disjoint classes, either exogenous or endogenous to others. Within each exogeneity class, a counterbalance could also be set up, and its solution quantifies the causality for the movement of the entire class. We also set up a similar counterbalance for the endogeneity class. This balance is subject to external influence from exogeneity classes. The solution explains the reason and responsibility of the movement of any particular endogenous variable. However, the equilibrium is set at the long-run level, not at the short-run dynamics, because a real causal factor could affect over a long time and involve many intermediaries. We convert the causal identification into a hypothesis testing statistical decision problem by estimating the expected effect and uncertainty. A few advantages could explain the method's superior performance over other methods. First, by admitting the indirect influences from the third parties, the limit bilateral impacts are consistent and resolve internal conflicts. Secondly, the formula and estimation procedures are simple and easy to implement. Also, as illustrated in the simulation studies, the performance is not perfect but good enough for a small-scale analysis. Next, by normalizing the sum of rows in Ω, we lay all the noise of Ω in a Procrustean bed, allowing only zero-sum for the noise in each row. The procrustean ruled noise is further offset in the sum of infinitely linked chains of causality in which both positive and negative noise is present. Finally, the causal inference called the causality distribution for an effect variable is a byproduct of the causal discovery. The distribution fairly allocates the responsibility of the effect to the source factors. In practice, for example, one could apply it to taxation based on the causality distribution of the pollution, climate change, income, or social inequality. One could extend the solution from different angles, which are worth further development. First, linear functions may be too simple in a complicated data analysis, though linear methods could be tried first before moving on to nonlinear alternatives. However, formulating a stochastic matrix like Ω is still possible; the matrix should summarize the directional influence between any two time series, even in a nonlinear setting. The h-step forecast error variance decomposition ( 9) is a slightly different alternate, especially when h is the lag length of the VAR or the sample size of the data. Secondly, the algorithms have space to improve the accuracy. The 83% precision for the basic causal structures could imply a muchdiscounted accuracy for a complex situation with many variables. The causal relation may be tested for any pair of time series, ignoring the computational cost and the class property. Thirdly, the transition matrix Ω captures the variability of the time series vector. Positive and negative changes, however, could have asymmetric roles in the data movement process. Next, not all data structures work well in these algorithms, for example, for extremely weak exogeneity. Of course, the results also depend on the lag length selection and other specifications of the VARs. Finally, in the algorithms, we capitalize on the bootstrapping for the standard errors of the π. The data re-sampling technique not only involves more computational cost but also brings new sampling errors. A direct way is to calculate the standard errors using ( 7)-( 8) or ( 13)-( 14) to avoid costs and new errors. Most economists at one time or another have probably found themselves in an unconscious position of mixing up correlation with causation, though not because they consciously know the differences. At such times it may be convenient to have an illustration at hand to show the limitations of regression, rather than a list of causality definitions or inferences. What, then, is causality? The answer, it appears, is that any argument carried out with sufficient precision is mathematical logic. VAR or regression modeling is not necessary, but consistency within any causal chain is. Seeing through the appearance to perceive the essential determinants is another necessity. Lastly, the argument matches people's common sense of causality, such as those pictured in Figures 2 to 6.
For these considerations, the foregoing causality distributions may serve as two practical solutions. By the Lyapunov Theorem, the solution to (6) exists if and only if the time series y .t is asymptotically stable. So the limit of the denominator in (9) exists as h → ∞, and the limit is an entry of X. The numerator also has a limit because it is increasing with h and bounded by the denominator. The ratio is unconditional because both the numerator and denominator of ( 9) are already unconditional. ♦

First, say, the pattern of (2) has the sizes of

That is, there are n-m variables in the exogeneity class and m in the endogeneity class. Clearly, the product of any two matrices with the pattern also has the same pattern. Secondly, the moving averaging representation of y t+h is

where the matrix φ s can be recursively deduced from

with the starting values φ 0 = I n and φ s = 0 n 0 n for any s < 0. By the mathematical induction, clearly, φ s has the same pattern for any s = 0, 1, • • • . Next, Σ's Cholesky component L j for any endogenous y j has zeros in its first n -m elements as y j is after these n -m exogenous variables in the VAR. With the pattern of φ s , φ s L j also has zeros in its first n -m elements because φ s L j is a linear combination of the last m columns of φ s . Therefore, (φ s L j )(φ s L j ) has a zero block in its upper right corner. Finally, Ω has the same pattern because

and its contribution by ỹj is

In the diagonals of ( 23) and ( 22), the only changes from the original VAR in (1) are the first elements, both multiplied by β 2 . Therefore, the ratios between the two diagonals remain the same as ω ij for all i = 1, 2, • • • , n. ♦

Using (11), we write

When we make a small perturbation ∆Ω to Ω, the new causality distribution π + ∆π satisfies the counterbalance equation of

subject to ∆Ω1 n = 0 n and ∆π1 n = 0. After subtracting π = πΩ from (24), we get ∆π[I n -Ω -∆Ω] = π∆Ω and its first-order approximation

Let us increase ω ji by ∆ω ji and calculate the effect of the change on π. The elements of ∆Ω are all zeros except the jth row. After the increase of ω ji , the sum of the jth row is 1 + ω ji . To maintain the unit sum of the jth row, divide the jth row by 1+ω ji . Compared to their original values, thus, the other elements in the row decrease roughly proportionally. Thus, the linear approximation of the jth row of ∆Ω is

where β is the row vector delimited by the parenthesis. Thus, ∆Ω = ∆ω ji e j β + o(∆ω ji ) where e j is the jth column of I n . Dividing (25) by ∆ω ji and letting ∆ω ji → 0, we have the derivative of π with respect to ω ji :

or, in the usual column representation of unknowns,

To solve dπ dωji from (26) and the identity dπ dωji 1 n = 0, we write the augmented matrix for dπ dωji as

In the (n + 1) × (n + 1) matrix, the sum of the last n rows is a zero row vector. We delete the i + 1st row to drop the collinearity. We also move the ith column to the first without changing the order of other columns; this movement reorders the ith variable in dπ . The result is the augmented matrix for d(πi,π-i) dωji ,

to the left side of (28) to get the following new augmented matrix

all its elements are non-negative. Therefore, the first equation in (29) implies that

Also, the second equation of (29) implies that

By the Perron-Frobenius theorem, the largest eigenvalue (in magnitude) of T e is a real positive number, say ρ, and its eigenvector v is a positive vector. So, T e v = ρv. Also, 0 ≤ min i∈T j∈T T e (i, j) ≤ ρ ≤ max i∈T j∈T T e (i, j) ≤ 1. We claim that ρ < 1. Otherwise, if ρ = 1, we define the set of elements which have the largest value in v:

Clearly, the set Ψ v is not empty. For any i ∈ Ψ v , for the equality

to hold in T e v = ρv, we must have v j = v i whenever T e (i, j) > 0, or j ∈ Ψ v whenever T e (i, j) > 0. Thus, Ψ v is an independent exogeneity class because T e (i, j) = 0 for any i ∈ Z and j ∈ Z. We get a contradiction to the causaltransient class T . Therefore, ρ < 1. We write T e as the Jordan normal form T e = F J e F -1 where F is an invertible matrix. The Jordan matrix J e has T e 's eigenvalues on the diagonal and 1s or 0s on the superdiagonal. Thus, lim First, if there exists a vector x ∈ R m such that (I -T e )x = 0 m , then x = T e x and, as a consequence, x = T t e x for all t = 1, 2, • • • . Letting t → ∞ gives that x = 0 m by Theorem 6. Therefore, I m -T e is invertible. Next, if T = ∅ and y .t contains no exogeneity class, then N itself is a causaltransient coalition and Ω t converges to a zero matrix as t → ∞.
So Ω t 1 n → 0 n . However, this contradicts to Ω t 1 n = 1 n for any t = 1, 2, • • • . Finally, by π = πΩ and the causal-transiency of T , π i = j∈T π j T t e (j, i) for all t ≥ 0. Letting t → ∞ results in π i = 0, by Theorem 6. ♦

In Ω = In the iteration (19), µ (t) is an increasing vector as t → ∞ because µ (t+1) -µ (t) = T e (µ (t) -µ (t-1) ) = • • • = T t e (µ (1) -µ (0) ) = T e µ (1) ≥ 0 m . Next, clearly, µ (0) ≤ 1 m . If we assume that µ (t) ≤ 1 m for some t > 0, then by (19), By the mathematical induction, µ (t) ≤ 1 m for all t. Therefore, lim t→∞ µ (t) exists. The limit satisfies (18) and is between 0 m and 1 m . If (18) has two different solutions µ c and δ c , without loss of generality, say, τ = max j∈T {µ c j -δ c j } > 0. Let

Then for any i ∈ Z, for the equality

to hold if and only if µ c j -δ c j = τ whenever ω ij > 0. That is, j ∈ Z whenever ω ij > 0 and ω ik = 0 whenever k ∈ Z. Therefore, Z is an exogeneity class and we get an contradiction to the causal-transient class T . This contradiction indicates that τ can not be positive and the solution to (18)

Theorem 10. If π c j is j's causality within the class C, then lim t→∞ Ω t (i, j) = π c j µ c i for any i ∈ T and j ∈ C.