Hongxia Chen and Jimin Ye

(School of Mathematics and Statistics,Xidian University,Xi’an 710126,China)

Abstract:Recently, several approaches have been proposed to discover the causality of the time-independent or fixed causal model. However, in many realistic applications, especially in economics and neuroscience, causality among variables might be time-varying. A time-varying linear causal model with non-Gaussian noise is considered and the estimation of the causal model from observational data is focused. Firstly, an independent component analysis (ICA) based two stage method is proposed to estimate the time-varying causal coefficients. It shows that, under appropriate assumptions, the time varying coefficients in the proposed model can be estimated by the proposed approach, and results of experiment on artificial data show the effectiveness of the proposed approach. And then, the granger causality test is used to ascertain the causal direction among the variables. Finally, the new approach is applied to the real stock data to identify the causality among three stock indices and the result is consistent with common sense.

Keywords:time-varying causal model; independent component analysis (ICA); granger causality test; causality inference

1 Introduction

Inferring causalities between variables is a primary goal in many applications. Literatures are mainly concerned with two different causality frameworks: the underlying causal structure and the formalized causal model. Under appropriate assumptions, the conditional independence relied on detection of the latent causality between variables had been solved[1-2]. Formalized causal models often give a more useful and accurate tool to establish the causality and it links up the time sequences and causal find. On the basis of the Granger[3], Granger causality in time series falls into the framework of causality discovery combined with the temporal constraint that the effect cannot precede the cause[4]. Several formalized causal models with some particular assumptions had been established[5-10], as well as the models with non-Gaussian noise[11-13], and methods to find causal ordering in these causal models had been proposed, in which there existed an ICA based method[14-15]. However, these formalized causal models are often time-invariant or the analytic data is restricted to be in equilibrium states.

In most realistic circumstances, particularly in economics, neuroscience and weather analysis, the causality may alter with time. Once one always makes use of an invariable causal model, it may cause misunderstanding about causal relations. By introducing time information, a time-varying linear causal model (TVLC) with changing coefficients between the observed processes has been established, where the noise term of the model is assumed to be Gaussian variable[16]. Note that in most practical applications, the non-Gaussian noise is widely applied, TVLC with the non-Gaussian additive noise is considered in this paper, and the main attention is placed on how to estimate the model accurately from data as well as ascertain the causal directions among the variables.

Taking into account the time instantaneous and lagged effects, there are two kinds of methods for estimating the time-varying models. One uses adaptive filters or sliding windows[17-18], due to the assumption that the lagged causal coefficients and instantaneous causal coefficients are independent of the time, while causal effects vary over time these approaches might cause large errors of estimation. The other inverts the model estimation to non-parametric estimation problem via Gaussian Process prior[16], under the assumption that the noise term is a random variable which is Gaussian with independent identically distributed(i.i.d.).

Regarding the TVLC with the additive non-Gaussian noise, in the paper, an original ICA based two-step approach is raised to solve the new model. The model is separated into two parts, one consists of the lagged coefficients and the instantaneous coefficients simultaneously, it is essentially a convolution model, thus, in the first stage, the mixture of the lagged coefficients and the instantaneous coefficients is estimated via deconvolution. The other part contains the instantaneous coefficients and noise term only, it is just ICA model. In second stage, the instantaneous coefficients are estimated with ICA method. The lagged coefficients are estimate eventually using estimated results of two stages.

The primary contents of this paper consist of two aspects. The TVLC model to interpret the time-varying causal effects between the observational data is generalized to the case with the additive non-Gaussian noise. An ICA based two stages method is proposed to figure out the time-varying coefficients of the model, and the causal direction between the variables is ascertained through the granger causality test.

2 Definition of Time-varying Linear Causal Model

We first concisely retrospect usual formalized causal model[2]ahead of presentation the time-varying formalized causal model. In normal shape, a formalized causal model is consisted of a series of following formulas

xi=fi(pai,si),i=1,2,…,N

(1)

By including the time effect, the functional causal model (1) can be extended to a time-dependent functional causal model, which explains time-varying causal influences among the observed variables. As shown in Fig.1.

Fig.1 Time-varying causal’s graph

The temporal information is counted as a corporate reason to other observational variables. Dotted lines express the effects of timeTto the observed variables. What one mostly interested in is the causal structure amongx1,x2andx3which is represented by solid lines.

Assume thatx(t)is a multidimensional time sequences,x(t)=(x1(t),x2(t),…,xN(t))Twith a limited dimensionN. In light of the work[16]and to capture time-varying causal relations, we propose thatx(t) obeys the following data generating models:

(2)

wheresi(t),i=1,2,…,Nare non-Gaussian variable that mutually independent and also independent of the reasonsxj(t-p),xk(t) andt.ai,j,p(t) andbi,k(t) represent the time-varying lagged causal coefficients and the instantaneous causal coefficients respectively, and they are assumed to change over time smoothly and slowly. Note that different from the existed Gaussian noise in other model,si(t)，i=1,2,…,Nin Eq. (2) are non-Gaussian variables and not limited to specific distributions.

3 Model Estimation

The part of model estimation includes the determination of causality and the estimation of the model’s coefficients. Although our model belongs to causal graph model in essence, it is still a regression model. The causal graph model can infer not only from cause to result, but also from result to cause. In this paper, according to the model definition, the causality diagram which only infers from cause to result is considered, and then the granger causality test is used to ascertain the causality.

3.1 Determination of Causality

In so far as the time sequences, the causality of granger between two variablesXandYdefines as follows. In the condition of the bygone information ofXandYare known, the forecast ofYis superior to that only depends on the past information ofY. In other words, variableXis conducive to interpret prospective variation of variableY. Then variableXis considered to be the granger causal ofY.

The time sequences must be steady is a precondition of the granger causality test. So it is necessary to deal with the time series’ stationarity first by unit root test. A common method of unit root test is the augmented Dickey-Fuller test. If the sequences have a roots of unity, then it is not a stationary sequences. The granger causality test can be proceeded after proving that the time series is stationary process. Both the unit root test and the granger causality test can be done through the program named Eviews7.

The specific theory of the granger causality test is as follows, as for two stationary series {xt} and {yt} satisfy the regression model:

(3)

wherecis the constant term, andnis the time lagged term. Testing the variation ofxis not the reason ofyis equivalent to theFtest of the statistical null hypothesis

H0:β1=β2=…=βn=0

The statistical test value is:

where RSS1is the summation of squared regression errors of formula (3),RSS0is the summation of squared regression errors of formula (3) under null hypothesis,Nis the sample size. StatisticsFfollow standardFdistribution. Once theFinspection value is larger than the threshold of standardFdistribution, then refuse the original hypothesis, and signifies thatx’s variation is the cause ofy’s. Or accept the original hypothesis and signifies that the variation ofxis not the reason ofy[19].

3.2 Determination of Model’s Coefficients

In this part, a two stage approach based on ICA is proposed to estimate the coefficients of TVLC model (2). The TVLC model (2) written in matrix form as below.

(4a)

or

(4b)

wherex(t)=(x1(t),x2(t),…,xN(t))T,Iis theN×Nidentity matrix,B(t) with elementsbi,k(t) is a matrix, which implies the instantaneous causal relation and could be arranged to be a strict lower triangular,Ap(t) with elementsai,j,p(t) is a matrix,s(t) with elementssi(t) is a vector.

Denote

n(t)=(I-B(t))-1s(t)

then, Eq.(4) turns into:

(5a)

n(t)=x(t)-z(t)

(5b)

Then the estimation of the causal coefficientsAp(t),B(t) and sourcess(t) in model (4) can be estimated in two steps according to Eq.(5).Cp(t) can be estimated through the Eq.(5a), andB(t) can be estimated by combine the Eq.(5b) andn(t)=(I-B(t))-1s(t). In summary, the estimations ofAp(t) andB(t) are obtained. In the following, we describe two parts in detail.

Step1Estimate mixture coefficientsCp(t) of the lagged coefficients and the instantaneous coefficients.

It is easy to see that Eq. (5a) is a convolution model. Thus the problem of the model’s coefficients estimation convert to a deconvolution problem. In order to solve the problem, find an inverse filterw(t) is necessary, and then recover the input datax(t) from the output dataz(t). A deconvolution algorithm based on minimum entropy[20]is applied to estimate each column of the model’s coefficientsCp(t).

The algorithm is described as follows.

(1)To initializew(0)=1,i=1.

(2)To compute iterativelyx(t)=w(t)(i-1)*z(t).

(4)To computew(i)=C-1b(i), whereCis the autocorrelation matrix of seriesz(t).

Else, leti=i+1 and then go back to (2).

The estimation of the mixture coefficientsCp(t) are obtained. Due to the mixture of the lagged coefficients and the instantaneous coefficients, the instantaneous coefficients are estimated next and the estimation of the lagged coefficients are obtained finally.

Step2Estimate the instantaneous coefficientsB(t) from the noise terms.

In Eq. (4), let:

n(t)=(I-B(t))-1s(t)

(6)

wheren(t)=(n1(t),n2(t),…,nN(t))T. Because of the assumption in models that ?t,si(t),i=1,2,…,Nare mutually independent, and ?i,si(t) is i.i.d. noise term. Thus (Ι-B(t))-1can be seen as a mixture matrix and Eq.(6) can be considered as a time-varying mixture ICA model. The observed variablen(t) is known, which is computed asx(t) in Eq. (4) minus its estimationz(t) in Eq. (5), we need to estimate mixture matrix (I-B(t))-1ands(t). For convenience, denoted (I-B(t))-1byA(t), then the unknown variables can be estimated by the following method[21].

Suppose the general mixing process is

x(t)=A(t)s(t)

where the component ofs(t) areNindependent random variables. In the case thatAis an unknownN×Ntime independent matrix, this is also called blind signal segregation. Performing blind signal segregation means to seek the matrixJso as toy(t)=Jx(t) gives a reconstruction ofsi(t). Let

C0=E{xxT}=AK(0)AT

C(τ)=E{x(t)x(t-τ)T}=AK(τ)AT

are the two-point correlation matrix at equal times and 2nd order cumulative matrix at time lagτrespectively. A possible practical method for obtainingJis to primarily fulfill a principal component analysis onC0and compute the whitening matrixV, and then to find the orthogonal matrixEwhich diagonalizesCz(τ)=E{z(t)z(t-τ)T}, wherez(t)=Vx(t) is the whitened mixture, thenA=VEandK(τ) are estimated.

Regarding toA(t) is time-varying matrix, which is varying smoothly and relatively slowly, i.e., in any time period with scaleT, there being a time sizeTso as to the hybrid matrix is nearly constant,i.e., the norm of dA/dtis small compared to 1/T. There-in-after, a quantityA’s average means an average on the time interval [t-T,t], and can be computed as:

Hence, we get:

C0=E{x(t)x(t)T}=

AT(t-t')

(7)

C(τ)=E{x(t)x(t-τ)T}=

AT(t-τ-t')

(8)

On condition that the time correlation ofAis feeble, we use an estimation: for every time period [t-T,t], one computes a matrixJtsuch thatJtachieves signal segregation in mean on that time period.

For the sake of obtaining a well estimation of the hybrid matrix for a prescribed time interval [t-T,t], we estimateA(t-t') with a linear expansion int'fort'from 0 toT. Then we will compute two matrices instead of computing one matrixJ. The assumption on the slowly change of the hybrid matrix signifies, at 1st order,A(t-t') written as below fort'in 0 andT

(9)

E{x(t)x(t-τ)T}=

(10)

And get

(11)

where

NoteK(0)(τ)=K(τ),C(0)=C0.We can see thatC(τ) is dissymmetric for nonzeroτfrom Eq.(11). From the aforementioned equation, it shows in the process of estimatingC0andC(τ), a new unknown variableK(τ) is involved. To address this issue, we presume that under the integrals, the source terms can be replaced by their average,

(12)

(13)

the cumulants. We obtain the below equations:

(14)

(15)

(16)

After above two steps, all causal coefficients include the lagged causal coefficients and the instantaneous causal coefficients of original model are obtained. At the same time, the causality between variables are also determined. That is to say, the model estimation is finished. The effectiveness of the proposed method will be illustrated in simulation section.

4 Simulations

1)Simulation 1. 2000 data are generated from the below equation which contains the time delayed and contemporary causal effect, and possesses slickly and slowly varying coefficients:

a1,1,1(t)=0.3·(sin(t)+1.1)
a1,2,1(t)=0.2·cos(t)+0.05
a2,1,1(t)=0.2·sin(t)+0.1

a2,2,1(t)=0.5·(cos(t)+0.05)

b2,1(t)=0.2·cos(t),si(t)～U[0, 1],i=1, 2

The proposed TVLC model is used to match the data and we first obtain the estimation of the mixture of the lagged coefficientsai,j,1(t),i,j=1, 2 and the instantaneous coefficientsb2,1(t). Then we estimate the instantaneous coefficientsb2,1(t) through the mixture of noise terms and recover the sources through the proposed ICA based method. Finally, the estimation of the lagged coefficients are obtained. According to the estimated coefficients, we compare the true value and predicted value ofxi(t), Fig.2 shows the comparison between them with time scale 20. The comparison between the sources and its estimations with time scale 20 are shown in Fig.3. It is well known that there are two indeterminacies in ICA. The indeterminacy of signs appears in the second figure in Fig.3. We also show the estimates value ofxi(t) with time scale 100 in Fig.4 and the recovery of sources with time scale 100 in Fig.5. For more time scale, the figure is too difficult to show the estimated effectiveness. Thus they are not shown out. From figures, the proposed method is obviously effective. Whether the comparison between the true value and predicted value ofxi(t) or the comparison between the sourcessi(t) and its estimation, the result has a little bit of errors generated in the progress of the matrix computing. However, the overall trend are well.

Fig.2Thecomparisonbetweenthetruevalueandpredictedvalueofxi(t)withtimescale20

Fig.3Thecomparisonbetweenthesourcesanditsestimationswithtimescale20

Fig.4Estimatesvalueofxi(t)withtimescale100

Fig.5Therecoveryofsourceswithtimescale100

Finally, we use granger causality test to ascertain the causality betweenx1(t) andx2(t) by EViews7. For this generated model, the contemporary causality between variablesx1andx2is known from its form. The model shows that variablex1is the instantaneous cause to variablex2. Thus, the granger causality is used here to verify the relationship between two variables. First judging the stationary ofxi(t) using unit root test. The statistical value shown in Table 1 is larger than the critical value and that meansxi(t) is not a stationary series. After the first difference of the seriesxi(t), the statistical value is smaller than the critical value and means that the 1st difference series are stationary series (as shown in Table 2). Appling granger causality test to the final series, for the model has one delay, we get thatx1is the grange cause tox2which agrees with our generated model. The statistical result is shown in Table 3.

Table 1 The stationary test statistical value of xi(t)

Table 2 The stationary test statistical value of the first difference of xi(t)

Table 3 The granger causality test’s result

2) Real Data Test. The daily adjust closing rates of three stock indices are chosen from 03/16/2005 to 07/31/2014, which are US’s NASDAQ, UK’s FTSE, and Japan’s N225. The three exponents are the main stock indices all over the world, so we have a great interest in the causal relations among them. We analyze the return series of the stock indices.

We first consider all possible pairs and test whether the granger causality test is able to find plausible causal direction. Only when the causality among three stock variables is dissymmetric can the model in Section 2 be established and the coefficients can be estimated. The granger causality test commands that the series must be stationary, so we first test the stability of the three stock indices series by the unit root test. The unit root test’s results are shown in Table 4. Three series are all stationary. The granger causality test can be done to judge the right causal order between stocks. Through Eview7 the result of granger causality tests are shown in Table 5. According to the results, the test result of pair N225 and FTSER show that FTSER effect N225, pair NASDAQ and FTSER shows NASDAQ effect FTSER, and pair NASDAQ and N225 shows NASDAQ effect N225. Thus, the proposed approach favors the following causal relation NASDAQ→FTSE→N225. This really accords with the sequence on account of time lags: the time domains consistent with NASDAQ, FTSE and N225 are UTC-5, UTC, UTC+9, respectively. That is to say the granger causality test result satisfies the condition of the TVLC model. We then fit our TVLC model and estimate the various causal effects. Some of the estimated lagged causal coefficients are shown in Fig.6 and the contemporary causal influence are shown in Fig. 7.

Table 4 The unit root test results of three stock indices

Table5Theresultofgrangercausalitytestsofthreestockindices

Null HypothesisF-StatisticProb.N22R does not Granger Cause FTSER5.624 150.017 8FTSER does not Granger Cause N22R391.222 003.E-80NASR does not Granger Cause FTSER229.006 003.E-49FTSER does not Granger Cause NASR2.297 520.129 7NASR does not Granger Cause N22R804.727 005E-151N22R does not Granger Cause NASR0.020 270.886 8

The variation tendencies in Fig.6 and Fig.7 show that the causal influence of one-day delayed from FTSE to N225 is fairly little, while other lagged effects are slightly significantly. In contrast, the one-day delayed causal effects of stock itself are more significant than the causal effects from one stock to the other. The contemporary causal influences from NASDAQ to FTSE, and from FTSE to N225 are evident, which are concomitant with the effect because of the time lag. The causal coefficients become much larger fromT2toT3, this may be the result caused by the financial crisis of 2008.

Fig.7Instantaneouscausalintensities’estimatedresultonshares

5 Conclusions

In existing models, either the coefficients are fixed or the noise term is Gaussian, while in the proposed TVLC model, the coefficients are time-varying and the noise is non-Gaussian. In economics, neuroscience and climate analysis, the causality may alter with time and the noise is non-Gaussian, the proposed TVLC model has more widely use in practice and therefore has great research value. A two-step approach based on ICA to identify the time-varying causal influences is proposed by reserving the temporal information as a corporate reason. Under the appropriate assumptions, the method has a well estimated result. This work links the ICA based method with causality discovery and can be further studied.