Nonlinear Granger Causality Analysis for Multivariate Time Series Using HSIC-GL Model
-
摘要:
因果分析是数据挖掘领域重要的研究课题之一.由于传统的Granger因果模型难以准确识别多变量系统的非线性因果关系，本文提出一种基于Hilbert-Schmidt独立性准则（Hilbert-Schmidt independence criterion，HSIC）的组Lasso模型的Granger因果分析方法.首先，利用HSIC将输入样本和输出样本映射到再生核Hilbert空间，克服了传统的Granger因果模型不能应用于非线性系统的缺陷.然后，建立具有组Lasso约束的回归模型，对多变量及其组派生变量进行因果分析，并采用贝叶斯信息准则进行模型选择，避免了人为设置滞后阶数和正则化参数.最后，根据HSIC-GL模型的回归系数和显著性检验结果，实现了多变量时间序列之间的非线性因果分析.通过对非线性和混沌系统的仿真实验，验证了该方法的有效性.最后将其应用于沈阳空气质量指数（AQI）和气象时间序列的因果关系分析.
-
关键词:
- 多元时间序列 /
- Granger因果 /
- Hilbert-Schmidt独立性 /
- 准则 /
- 组Lasso
Abstract:Causality analysis is an important research topics in the field of data mining, but traditional Granger causality models have difficulty accurately identifying the nonlinear causality of multivariable systems. We propose a novel Granger causality analysis method based on the HSIC and group Lasso (HSIC-GL) model. Firstly, we use the Hilbert-Schmidt independence criterion (HSIC) to map the input and output samples into the Hilbert space of the reproducing kernel, which overcomes the inability to apply the traditional Granger causality model to nonlinear systems. Then, we establish a regression model with group Lasso constraints, which implements a causality analysis between multivariate and group-derived variables. The Bayesian information criterion is used for model selection, which prevents the artificial setting of the lag order and regularization parameters. Lastly, based on the regression coefficients and the results of significance tests of the HSIC-GL model, a nonlinear causality analysis is performed on the multivariable time series. The effectiveness of the proposed method is verified by the results of simulations of nonlinear and chaotic systems. We successfully applied this method to the air quality index and meteorological time series in Shenyang, China.
-
0. 引言
多元时间序列广泛存在于工业、气象、医学[1-3]等多个领域. 时间序列分析利用时间序列数据、应用数理统计等方法来挖掘系统潜在信息，揭示系统未来发展规律，是研究系统演化机理、建立系统模型的重要手段. 相较于单变量时间序列，多元时间序列往往包含更加丰富的系统信息，能够更加准确地揭示系统演化规律. 但是，多变量系统变量间的影响关系更加复杂，数据中难免会存在着与预测对象不相关的无关变量以及对预测对象作用相同或相似的冗余变量. 无关变量和冗余变量不仅会增加预测模型建立的难度，还会延长训练时间，对模型的预测效果产生负面的影响[4]. 因此，研究时间序列变量间的影响关系，为模型选择合适的输入变量具有重要的研究意义.
因果关系分析已经广泛应用于揭示多变量系统间的相互影响关系，可以有效识别复杂系统中无关变量和冗余变量. 在多变量时间序列的建模预测中，因果关系解释了因变量对于目标变量的影响关系. 通过因果关系分析可以有效地剔除无关和冗余变量，为模型选择合适的输入，达到建立准确的预测模型、提高预测精度的效果.
Granger因果[5]是一种常用的因果关系分析方法. 它通过建立线性自回归(vector autoregressive，VAR)模型来揭示变量间的相互影响关系. Grange因果基于可预测思想：对于两个时间序列X和Y，如果X历史信息的加入有助于减少Y的预测误差，则可以说存在X→Y的Granger因果关系. Granger因果一经提出就受到了学者们的广泛关注[6]. 但由于其只适用于二变量线性因果关系分析，具有很大局限性，因此国内外学者提出了大量改进模型[7]. 针对多元时间序列因果关系分析，Geweke等向VAR模型中加入条件变量，提出了条件Granger因果(conditional Granger causality，CGC)模型[8]，可以有效区分直接和间接因果关系. Siggiridou等在条件Granger因果模型的基础上，采用Backward-in-time方法对每个变量的滞后阶数进行限制选择，提出了mBTS-CGCI(modify backward-in-time selection CGCI，mBTS-CGCI)方法[9]. 有效降低了VAR的模型阶数，可以实现高维时间序列的因果关系分析. 然而，由于条件Granger因果模型涉及大量参数估计，面对大规模数据集时存在较大困难. 另外，在模型求解时利用最小二乘方法求解得到的解并不稀疏. 而Lasso回归通过添加l1范数惩罚项，收缩回归系数，可以产生稀疏的变量选择结果. 因此Arnold等[10]应用Lasso回归结合Granger因果模型，提出了Lasso-Granger因果模型(Lasso-GC)，该方法通过建立一个回归模型就可以实现系统所有变量的因果关系分析，极大地降低了计算复杂度. 此后，Lasso-GC及其改进模型被广泛应用于因果关系分析[11].
随着研究的深入，学者们发现非线性系统广泛应用于各个领域. 而基于VAR模型的线性Granger因果模型可能无法准确获取非线性系统变量间的内在影响关系. 因此，研究学者也提出了大量非线性Granger因果改进模型. Ancona等基于Granger因果关系的预测思想，提出基于径向基函数的非线性Granger因果模型[12]，实现了二变量的非线性因果关系分析. Marinazz等[13]提出基于核方法的Granger因果(kernel Granger causality，KGC)模型. 该模型应用核函数将原始数据进行非线性映射，在再生核Hilbert空间(reproducing kernel Hilbert space，RKHS)中进行Granger因果关系分析. 除了基于核方法的非线性Granger因果模型，还有基于Copula函数的Granger因果模型[14]、基于神经网络的Granger因果模型[15]等. 此外，有学者基于信息理论，提出了条件熵[16]、基于混合嵌入的偏互信息(partial mutual information from mixed embedding，PMIME)[17]等因果关系分析方法. 它们都能有效地分析非线性系统中变量间的相互作用关系. 但这些基于信息理论的因果分析方法涉及概率密度函数的计算，当数据维度增加或样本量较大时，大规模数据的概率密度函数计算将变得异常困难，并且其计算量会成倍增加. 因此，随着数据维度及规模的增加，如何高效识别复杂多变量系统中的非线性因果关系，具有极大的挑战和研究价值.
传统的Granger因果及其扩展模型只利用变量的衍生规律来解释因果关系，往往忽视了冗余信息. 另外，传统的非线性因果模型，如核Granger因果模型，仅仅将输入变量进行非线性映射，没有对目标变量进行非线性映射，从而可能导致缺乏同等变换的物理意义解释. 因此，为了解决上述问题，本文提出了一种新的基于Hilbert-Schmidt独立性准则和群组Lasso(HSIC-GL)的Granger因果模型，可以有效地识别多变量系统的非线性因果关系. 该模型利用HSIC将输入变量和目标变量同时进行非线性映射来实现Granger因果模型在非线性系统中的应用，借助HSIC实现类似于最大相关最小冗余的效果[18]，可以有效识别无关和冗余变量. 另外，当输出中涉及大量相似特征时，基于HSIC关联测度的特征选择模型可能会保留冗余信息[19]. 因此，建立HSIC-GL模型解决上述问题. 该模型可以通过建立一个回归方程实现多变量因果分析. 此外，本文使用贝叶斯信息准则(Bayesian information criterion，BIC)来确定模型阶数和群组Lasso的正则化参数，排除了人为因素的干扰. 最后，根据因果邻接矩阵并运用显著性检验框架实现因果关系分析. 实验结果也验证了本文方法的有效性和优越性.
1. 算法介绍
针对样本集的统计独立性评价问题，本节将介绍Hilbert-Schmidt独立性准则. 然后，在此基础上，介绍HSIC-GL模型.
1.1 Hilbert-Schmidt独立性准则
HSIC[20]是一种评价两个样本集统计独立性的核独立性准则. Gretton等[20]在RKHS中引入互协方差算子，并将HSIC定义为互协方差算子的Hilbert-Schmidt范数. 设X和Y为两个随机变量，从X和Y的概率密度函数中提取样本(X，Y). 定义两个非线性映射φ：X→F，ψ：Y→G，F和G分别表示X和Y的RKHS，X和Y对应的核函数分别为
(1) 
(2) φ：X→F和ψ：Y→G的互协方差算子Cxy：G→F定义：
(3) 式中，⊗表示张量积，μx=Exφ(x)，μy=Eyψ(y)，Ex，Ey，Exy表示数学期望. HSIC[20]定义为Hilbert-Schmidt范数的平方：
(4) 式中，Exx′yy′表示(x，y)和(x′，y′)的联合数学期望. 对于一个数据对Z={(xi，yi)|i=1，…，n}，HSIC的经验估计[20]为
(5) 式中，tr(·)表示矩阵的迹运算，K，L∈Rn×n分别为Kij=k(xi，xj)和Lij=l(yi，yj)的核矩阵.
1.2 HSIC-GL模型
Yamada等[18]在Lasso回归模型的基础上，提出了一种基于HSIC的非线性特征选择模型. 然而，在对回归模型的研究中，Yuan等[21]指出Lasso模型虽然在计算上是高效的，但其针对选择单个输入变量，而无法选择出具有解释意义的组派生变量. 并且当直接应用于回归模型时，Lasso倾向于根据单个输入变量的强度，而不是根据一组输入变量的强度进行选择，往往会导致选择出多余的变量. 群组Lasso在经验风险最小化的同时，对组内和组间采用不同的惩罚约束来进行变量选择，可以有效解决上述问题.
分析一个包含d维特征、样本数为n的矩阵X∈Rd×n与目标变量为Y∈Rn的相关性. 首先将d维特征分为G个组，dg表示第g组的特征维数. 然后，将原始输入和输出样本映射到RKHS，并得到Gram矩阵K∈Rn×n和L∈Rn×n，最后建立HSIC-GL的目标函数，如式(6)所示：
(6) 式中，‖·‖Frob表示矩阵的Frobenius范数，为矩阵所有元素平方和再开方.
式中的第一项可以进行展开求解，如式(7)所示：
(7) 式中，HSIC(L，L)是一个常数，与分析结果无关. 如果第k个特征Xk与Y具有很强的相关性，则HSIC(K(k)，L)和αk都取得较大值，而式(7)取得最小值. 反之如果Xk与Y相互独立，则HSIC的值接近于0，对应的系数αk也接近于零.
2. 基于HSIC-GL的Granger因果关系分析
本节将详细描述基于HSIC-GL的Granger因果分析模型(HSIC-GL-GC)，并推导出群组Lasso的求解算法. 结尾还将给出HSIC-GL-GC的伪代码和复杂度分析.
2.1 平稳性分析
Granger因果模型是针对平稳时间序列的因果关系分析方法. 因此，在进行建模之前，需要对数据进行平稳性分析. 平稳性分析最常用的方法为单位根检验，如果时间序列具有单位根，则表明该时间序列是非平稳的，反之则为平稳时间序列. 本文采用ADF (augmented Dickey-Fuller)[22]检验方法来检验时间序列平稳性，ADF如式(8)所示：
(8) 式中，α是一个常数，β是趋势项系数，m是自回归模型的阶数，εt是误差项. 在零假设下进行单位根检验，如果接受零假设H0∶ δ=0，则表明时间序列是非平稳的. 如果时间序列是非平稳的，则需要采用差分方法实现平稳化. 反之如果拒绝原假设，而接受备择假设H1∶ δ < 0，则表明时间序列是平稳的，可以直接进行Granger因果关系分析.
2.2 模型选择
在进行模型拟合时，模型阶数对拟合精度影响较大. 如果VAR模型中解释变量的最大滞后阶数p太小，则可能存在残差自相关，导致参数估计不一致. 适当增加p可以解决上述问题. 但是，如果p过大，则会导致估计参数过多，自由度严重降低，直接影响模型参数估计的有效性[23]. 另外，在群组Lasso模型中，正则化参数λ对变量的选择结果影响很大，如果正则化参数值过大，可能会错误剔除相关特征. 反之，如果正则化参数值过小，则可能保留不相关或冗余的特征. 因此，本文采用BIC[23]来自动选择合适的模型阶数和正则化参数λ，以消除人为选择的不确定性干扰. BIC如式(9)所示：
(9) 式中，L表示似然函数，k是模型有效参数的个数，n是样本个数. 其中，对于线性回归模型y=WTX+ε，即P(y|X，θ)~N(y|WTX，σ2)，其似然估计函数为
(10) 式中，
(11) 式中，第一项表示模型拟合的优良性，第二项表示对模型复杂度的惩罚力度，可以有效刻画构建模型的准确度.
2.3 HSIC-GL-GC求解
考虑一组长度为n的矩阵X∈Rn×d，Y∈Rn×1. 本文通过建立如式(6)所示的HSIC-GL模型，然后求解群组Lasso获得稀疏的变量选择结果.
本文采用高斯核函数对X和Y进行非线性映射，如下所示：
(12) 
(13) 式中，σx，σy为核宽，x′，y′为核中心. 通过式(12)、式(13)将其进行非线性映射，然后应用HSIC得到Gram矩阵，并利用Shooting算法[21]求解如式(6)所示的群组Lasso模型. 文[21]指出式(6)的第一项可以利用矢量化算子改写为如式(14)：
(14) 式中，vec(·)表示矢量化算子. 然后，根据Shooting算法[21]，式(6)解的充分必要条件如式(15)和式(16)所示.
(15) 
(16) 其中，(vec(K(g)))Tvec(K(g))=I. 由式(15)可以求得
(17) 式中，Sg=(vec(K(g)))T(vec(K)-vec(K)α-g)，α-g=[α1T，…，αg-1T，0，αg+1T，…，αGT]. 然后，通过迭代即可求得全部解，并根据该解得到因果邻接矩阵.
2.4 显著性检验
在因果关系分析中进行显著性检验能够有效消除偶然因素导致的虚假因果，使结果更具可信度. 进行显著性检验时，首先对原始时间序列Xk(k=1，2，…，d)进行随机移位，得到置换时间序列，然后进行置换检验，分析原始序列的因果关系与置换序列的因果关系的检验统计性. 最后，根据显著检验结果推断原始序列因果关系的有效性. 置换时间序列按式(18)构造.
(18) 式中，s表示移位因子. 显著性检验的完整性步骤如下：
步骤1 计算原始时间序列的因果关系，得到因果邻接矩阵.
步骤2 根据式(18)对原始时间序列移位，获得100组不同置换时间序列.
步骤3 分别对这100组置换时间序列进行因果关系分析，获得每一组置换时间序列的因果邻接矩阵.
步骤4 提出H0假设：置换时间序列的因果关系邻接矩阵与原时间序列因果关系邻接矩阵相同. H1假设：置换时间序列的因果关系邻接矩阵与原时间序列因果关系邻接矩阵不同.
步骤5 计算pvalue. 如果pvalue>0.95，即原始时间序列的因果系数大于95%置换时间序列因果系数，则拒绝显著性水平为α=0.05的原假设，证明了原始时间序列因果关系的有效性. 反之接受原假设，则表明得到的原始因果关系为虚假因果关系.
2.5 总结与讨论
首先，总结HSIC-GL-GC的伪代码，然后对其计算复杂度进行了分析. 基于HSIC-GL的Granger因果关系分析伪代码如算法1所示.
算法1 HSIC-GL-GC 输入：时间序列{X，Y}，其中X=[X1，X2，…，Xd]T∈Rd×n，Y=[y1，y2，…，yn]∈R1×n；
输出：因果关系邻接矩阵.
1. 平稳性检验：
2. for p=1：pmax(pmax是设定的最大滞后阶数)
3. 利用HSIC将X和Y映射到RKHS，分别得到对应的Gram矩阵K和L；
4. for λ=λmin：λmax
5. 根据式(6)构建回归模型；
6. 根据计算BIC；
7. end for
8. 根据计算得到的最小BIC值确定模型阶数和惩罚系数，建立模型并计算α；
9. end for
10. for k=1：d
11. if pvalue(k)>0.95，即通过显著性检验
12. 存在Xk→Y的因果关系；
13. else if
14. 不存在Xk→Y的因果关系；
15. end if
16. end for
17. 输出：因果关系分析结果.在实际应用中，算法的计算复杂度是一个重要的评价标准，它决定了算法在实际中的应用普及程度. 本节分析了HSIC-GL-GC的时间复杂度，并与mBTS-CGCI和Lasso-GC方法进行对比. 对于n个样本的时间序列数据{X，Y}，首先对d维输入变量进行相空间重构，得到b维输入特征. 然后建立一个包含b维特征的VAR模型，其计算复杂度介于O(bn)到O(b2n2)之间[10]. 设P为设定的模型阶数的选择范围，则确定模型阶数，并创建VAR模型的计算复杂度介于O(Pb2n2)到O(Pbn)之间.
对于mBTS-CGCI，因为需要对d维输入变量进行两两计算，所以其时间复杂度介于O(Pdb2n2)到O(Pdbn)之间.
Lasso-GC模型常用最小角回归算法求解，通过建立一个回归模型就可以实现所有输入与输出的特征分析. 因此Lasso-GC的时间复杂度介于O(Pb2n2)到O(Pbn)之间[10].
对于HSIC-GL-GC，计算Gram矩阵的时间复杂度为O(bn2)，而求解HSIC-Lasso模型的时间复杂度为O(n3)[24]. 假设将原始特征分为G组，则求解HSIC-GL的时间复杂度为O(Gn3)，因此HSIC-GL-GC的计算复杂度为O(bn2+PGn3).
从上述分析可以看出，本文方法的时间复杂度低于mBTS-CGCI方法，高于Lasso-GC模型. 虽然时间复杂度高于Lasso-GC模型，但仍处于合理范围之内. 相比于Lasso-GC，本文方法增加了核映射和群组Lasso约束，扩展了模型的应用范围，并且提高了算法的分析精度，可以广泛应用于多变量系统的非线性因果关系分析.
3. 仿真实验与分析
为了验证本文方法的有效性，选取两个多变量非线性系统进行仿真. 然后利用HSIC-GL-GC分析了中国沈阳的空气质量指数和气象时间序列. 将仿真结果与mBTS-CGCI[9]、Lasso-GC[10]、KGC[13]和PMIME[17]进行了对比.
3.1 多变量非线性系统
第一组标杆数据是一个包含4维变量的5阶非线性VAR系统VAR4(5)[25]. 由式(19)产生：
(19) 式中，εi(i=1，2，3，4)表示高斯白噪声序列. 可以看到，系统中真实存在的因果关系为X1→X3，X2→X1，X2→X3和X4→X2，其因果影响关系如图 1所示.
![]()
图 1 VAR4(5)真实因果关系(黑色表示存在因果关系)Fig. 1. The real causality diagram of VAR4(5) (Black indicates that there is causality)图 2和图 3是利用BIC算法选择模型阶数和惩罚系数的三维图像的二维平面展示. 子图(a)、(b)、(c)、(d)分别表示目标变量为X1、X2、X3、X4的实验结果. 最小的BIC值对应着最优参数，即实验中选择的模型阶数和正则化参数. 图 4是不同方法得到的因果邻接矩阵. 从图中可以看到，Lasso-GC和mBTS-CGCI得到了虚假因果X3→X1的非零因果系数. 此外，mBTS-CGCI没有识别出真正的非线性因果关系X1→X3和X4→X2. 类似地，Lasso-GC没有识别出的因果关系有X2→X3和X4→X2. 这是因为这两种方法只适用于线性系统的因果关系分析，对于非线性系统可能产生错误因果识别. 同样的虚假因果非零系数也发生在PMIME和HSIC-GL-GC上. 但是，不同的是这两种方法的虚假因果关系邻接系数远远小于真实因果关系的邻接系数，因此可以进行有效区分.
表 1为VAR4(5)的显著性检验结果. 从表 1可以看到，KGC、PMIME和HSIC-GL-GC的真实因果均通过了显著性水平为0.05的检验，而mBTS-CGCI在X1→X3和X4→X2，以及Lasso-GC在X2→X1、X2→X3和X4→X2的检验结果均未通过显著性水平为0.5的显著性检验，而虚假因果关系也存在类似的结果.
表 1 VAR4(5)显著性检验结果Tab. 1. The significance test results of VAR4(5)方法 mBTS-CGCI Lasso-GC KGC PMIME HSIC-GL-GC X1→X2 2 95 0 37 0 X1→X3 15 96 98 98 100 X1→X4 2 4 5 1 3 X2→X1 98 2 100 100 99 X2→X3 95 4 96 99 96 X2→X4 0 0 1 0 0 X3→X1 2 16 3 0 1 X3→X2 1 2 0 62 3 X3→X4 0 1 2 3 0 X4→X1 0 2 0 4 0 X4→X2 45 68 97 98 96 X4→X3 2 0 0 25 3 3.2 Henon非线性耦合系统
第二组仿真数据来自于Henon系统. Henon耦合系统是一个具有复杂关联的多变量混沌系统. 三变量非线性Henon耦合系统[26]的方程如下：
(20) 式中，c表示耦合强度. 在实验中，分别设置c=0和c= 0.5进行仿真实验. 当c≠0时，存在的因果关系为X1→X2、X2→X3，其真实因果关系如图 5所示.
![]()
图 5 Henon系统真实因果关系图(黑色表示存在因果关系)Fig. 5. The real causality diagram of Henon system (Black indicates that there is causality)图 6和图 7分别是当耦合强度c=0.5时，本文方法采用BIC选取模型阶数和惩罚系数的三维图像的二维平面展示. 选取BIC最小值对应的p和λ作为模型的阶数和群组Lasso的惩罚系数. 图 8是对耦合强度c=0.5的Henon系统因果关系分析的邻接矩阵图. 从图中可以看出，Lasso-GC对虚假因果X2→X1和X3→X1产生了非零因果系数，mBTS-CGCI则存在着X2→X1的非零系数. 这是因为这两种方法适用于多变量的线性因果关系分析，无法准确识别非线性系统的因果关系，并且其对于真实的因果关系的识别效果不明显，识别能力较弱. 而KGC、PMIME和本文方法都能够正确识别系统中存在的因果关系. 相比之下，本文方法计算得到的因果关系邻接矩阵数值更大，对因果关系的识别效果更佳. 虽然本文方法得到了X3→X2的非零因果系数，但其数值远远小于真实因果X1→X2和X2→X3的因果关系系数，因此能够进行有效区分识别. 当耦合强度c取其它非零值时，得到了类似的因果关系分析结果.
![]()
图 6 耦合系数为0.5的Henon系统模型阶数选择Fig. 6. The model order selection of Henon system with coupled strength 0.5![]()
图 7 耦合系数为0.5的Henon系统惩罚系数选择Fig. 7. The penalty coefficients selection of Henon system with coupled strength 0.5![]()
图 8 耦合强度为0.5的Henon系统因果邻接矩阵Fig. 8. The causal adjacent matrix of Henon system with coupled strength 0.5表 2是耦合强度分别为0和0.5时的显著性检验结果. 从表 2可看出，mBTS-CGCI对虚假因果X1→X2和X1→X3具有统计显著性，其它方法的显著性检验结果都验证了当c=0时Henon系统变量间的非因果性. 当耦合强度为0.5时，系统存在着X1→X2和X2→X3的因果关系，根据检验结果可以看出，mBTS-CGCI和Lasso-GC都识别出了虚假因果关系. 相比之下，本文方法能够正确识别变量间的因果关系，通过了显著性检验，证明了本方法的有效性. 同时，相比于PMIME方法，本文方法在计算难度及计算成本具有更加明显的优势，尤其是当数据维度增加或规模扩大时，PMIME涉及边缘概率密度的计算，将面临着极大的困难. 相比于KGC方法，本文方法避免了对于核参数和滞后阶数的人为设定，从而减少了分析方法的人为偶然因素影响. 综上可以看出，本文方法具有明显的优越性.
表 2 耦合强度为0和0.5的Henon系统因果关系显著性检验结果Tab. 2. The significance test results of Henon system with coupled strengths 0 and 0.5方法 mBTS-CGCI Lasso-GC KGC PMIME HSIC-GL-GC c=0 c=0.5 c=0 c=0.5 c=0 c=0.5 c=0 c=0.5 c=0 c=0.5 X1→X2 100 100 50 51 22 99 0 100 2 100 X1→X3 100 100 0 0 0 100 0 0 9 23 X2→X1 0 4 65 100 0 0 0 0 21 4 X2→X3 0 100 0 100 0 98 0 100 4 100 X3→X1 0 0 14 100 0 0 0 0 31 0 X3→X2 0 1 0 100 13 0 0 0 6 10 3.3 沈阳AQI和气象时间序列
本节将对沈阳市的AQI(air quality index)和气象时间序列进行因果关系分析，并建立预测模型对因果分析结果进行验证. 数据来源于UCI数据库，采样时间为2014年3月24日至2015年2月22日，采样间隔为8 h，共计11维，每一维变量包含886组数据，包括6维空气污染数据和5维气象数据，其变量含义如表 3所示. 利用因果分析方法找出PM2.5浓度变化的主要影响因素，剔除无关和冗余变量，保留相关变量，并将该相关变量作为预测模型的输入进行建模预测，可以有效降低预测模型的规模，提高预测精度.
表 3 沈阳AQI和气象时间序列变量表Tab. 3. AQI and meteorological time series variables of Shenyang编号 物理含义 1 PM2.5 2 PM10 3 O3 4 NO2 5 SO2 6 CO 7 温度 8 露点 9 气压 10 湿度 11 风速 首先，采用ADF检验和差分平稳化方法对数据进行预处理. 然后以PM2.5为分析目标，建立HSIC-GL模型，其模型阶数和惩罚系数均采用BIC确定，选择的结果如图 9所示. 最后，求解HSIC-GL-GC，根据模型的回归系数和显著性检验结果来确定变量间的因果关系，并建立预测模型，选择相关变量作为预测模型的输入子集.
几种因果分析方法的结果如表 4所示. 从表 4可以看出本文方法识别出对PM2.5具有因果关系的变量分别是PM 10、O3、SO2、CO和风速. 沈阳市位于辽河平原中部，是我国重要的重工业基地. 同时，随着汽车的普及，当地的空气易受到工业、汽车尾气排放及燃煤等污染. 工业及汽车尾气排放主要是CO、NOx、SO2和烟尘微粒等. PM2.5是空气中直径小于2.5 μm的一类微粒的统称，其产生途径主要包括直接排放以和二次产生. 二次产生是指大气中的气态物质，如SO2、NOx和VOCs等经过一系列复杂的物理化学反应产生PM2.5的过程. 在沈阳地区，由于工业发展和汽车的普及，空气质量主要受到工业废气及汽车尾气排放的影响. PM2.5往往由空气中的其它气态排放物经过一系列物理化学变化形成，其浓度变化更易受到如SO2、NOx等的影响. 此外，风在污染物的扩散和累积中起着至关重要的作用. 一般来说，风速越大，越有利于PM2.5的扩散，PM2.5的浓度越低. 反之，风速越小，PM2.5浓度越高. 所以PM2.5的浓度很大程度上也会受到风速的影响. 综上所述，可以看到HSI-GLasso-GC的分析结果与沈阳当地PM2.5浓度变化的影响关系是一致的，进一步验证了本文方法的有效性.
表 4 PM2.5预测结果Tab. 4. The prediction results of PM2.5方法 因变量(编号) RMSE MAPE SMAPE - 全部 41.822 8 0.111 1 0.106 1 mBTS-CGCI 4、9、10、11 35.861 8 0.093 0 0.089 3 Lasso-GC 2、3、7、8、10 45.114 1 0.083 6 0.080 0 KGC 2、3、4、9、10 31.858 0 0.074 6 0.072 3 PMIME 11 34.456 5 0.109 3 0.104 8 HSIC-GL-GC 2、3、5、6、11 26.725 3 0.065 3 0.063 8 在进行因果关系分析之后，选择出与目标变量具有因果关系的原因变量作为模型的输入进行建模预测，对分析结果进行进一步的验证. 本文建立回声状态网络(echo state networks，ESN)预测模型[27]，来分别分析每一种方法所选择因变量的建模预测效果，并进行30次实验取平均值以消除偶然因素的影响. 采用均方根误差(RMSE)，平均绝对百分误差(MAPE)和对称平均绝对百分比误差(SMAPE)三个指标来定量评价预测精度，定义如下：
(21) 
(22) 
(23) 式中，yi和
图 10和图 11是PM2.5的预测效果图. 可以看到本文方法能够更加准确跟踪PM2.5的变化趋势，拟合效果更好. 表 4是不同方法预测PM2.5的精度比较. 可以看到，本文方法在3个预测指标都取得了最小值，表明其预测精度最高，具有明显的优越性. 由于实际系统各个变量之间存在着复杂的非线性关系，基于线性模型的mBTS-CGCI和Lasso-GC在面对复杂的非线性系统时，无法正确地排除所有无关及冗余变量的影响，难以提高预测精度. 而KGC、PMIME和本文的HSIC-GL-GC方法能够识别复杂系统的非线性关系，但面对类似的更为复杂的多变量系统时，PMIME方法的计算难度及计算成本往往会成倍增加，这会大大限制其应用范围，难以推广使用. 相比之下，本文方法能够有效应用于多变量系统的非线性因果关系分析，并取得良好的分析结果. 综上，可以看出本文所提方法能够有效地识别多变量复杂系统的非线性因果关系，为预测模型选择合适的输入变量，提高预测精度，具有良好的现实应用价值.
4. 结论
针对多变量系统的非线性因果关系识别问题，本文提出了HSIC-GL-GC方法，可以有效识别多变量系统的非线性因果关系. 该方法具有以下特点：
1) 利用HSIC对原始时间序列进行非线性映射，将Granger因果进行了非线性扩展，克服了传统Granger因果模型难于应用于非线性系统的不足.
2) 利用群组Lasso进行约束，通过建立一个模型就可以实现多变量因果关系分析. 另外，群组Lasso可以有效地对组派生变量进行选择，获得稀疏的变量选择结果，避免选择出无关和冗杂变量.
3) 利用贝叶斯信息准则对HSIC-GL的模型阶数和惩罚系数进行选择，避免了人为设定的干扰，可以有效消除人为因素产生的偶然误差.
但是，该方法仍存在不能直接应用于非平稳时间序列的缺陷. 当原始序列非平稳时，平稳化的过程可能会使原始序列的实际意义发生变化，从而导致因果分析方法的解释意义发生改变. 未来，将研究因果关系分析的非平稳模型，致力于采用时变模型将HSIC-GL-GC进行非平稳扩展等研究.
-
![]()
图 1 VAR4(5)真实因果关系(黑色表示存在因果关系)
Figure 1. The real causality diagram of VAR4(5) (Black indicates that there is causality)
![]()
图 5 Henon系统真实因果关系图(黑色表示存在因果关系)
Figure 5. The real causality diagram of Henon system (Black indicates that there is causality)
![]()
图 6 耦合系数为0.5的Henon系统模型阶数选择
Figure 6. The model order selection of Henon system with coupled strength 0.5
![]()
图 7 耦合系数为0.5的Henon系统惩罚系数选择
Figure 7. The penalty coefficients selection of Henon system with coupled strength 0.5
![]()
图 8 耦合强度为0.5的Henon系统因果邻接矩阵
Figure 8. The causal adjacent matrix of Henon system with coupled strength 0.5
![]()
图 10 不进行因果关系分析的PM2.5预测图
Figure 10. Prediction graph of PM2.5 without causality analysis
算法1 HSIC-GL-GC 输入：时间序列{X，Y}，其中X=[X1，X2，…，Xd]T∈Rd×n，Y=[y1，y2，…，yn]∈R1×n；
输出：因果关系邻接矩阵.
1. 平稳性检验：
2. for p=1：pmax(pmax是设定的最大滞后阶数)
3. 利用HSIC将X和Y映射到RKHS，分别得到对应的Gram矩阵K和L；
4. for λ=λmin：λmax
5. 根据式(6)构建回归模型；
6. 根据计算BIC；
7. end for
8. 根据计算得到的最小BIC值确定模型阶数和惩罚系数，建立模型并计算α；
9. end for
10. for k=1：d
11. if pvalue(k)>0.95，即通过显著性检验
12. 存在Xk→Y的因果关系；
13. else if
14. 不存在Xk→Y的因果关系；
15. end if
16. end for
17. 输出：因果关系分析结果.
下载: 导出CSV
表 1 VAR4(5)显著性检验结果
Table 1 The significance test results of VAR4(5)
方法 mBTS-CGCI Lasso-GC KGC PMIME HSIC-GL-GC X1→X2 2 95 0 37 0 X1→X3 15 96 98 98 100 X1→X4 2 4 5 1 3 X2→X1 98 2 100 100 99 X2→X3 95 4 96 99 96 X2→X4 0 0 1 0 0 X3→X1 2 16 3 0 1 X3→X2 1 2 0 62 3 X3→X4 0 1 2 3 0 X4→X1 0 2 0 4 0 X4→X2 45 68 97 98 96 X4→X3 2 0 0 25 3
下载: 导出CSV
表 2 耦合强度为0和0.5的Henon系统因果关系显著性检验结果
Table 2 The significance test results of Henon system with coupled strengths 0 and 0.5
方法 mBTS-CGCI Lasso-GC KGC PMIME HSIC-GL-GC c=0 c=0.5 c=0 c=0.5 c=0 c=0.5 c=0 c=0.5 c=0 c=0.5 X1→X2 100 100 50 51 22 99 0 100 2 100 X1→X3 100 100 0 0 0 100 0 0 9 23 X2→X1 0 4 65 100 0 0 0 0 21 4 X2→X3 0 100 0 100 0 98 0 100 4 100 X3→X1 0 0 14 100 0 0 0 0 31 0 X3→X2 0 1 0 100 13 0 0 0 6 10
下载: 导出CSV
表 3 沈阳AQI和气象时间序列变量表
Table 3 AQI and meteorological time series variables of Shenyang
编号 物理含义 1 PM2.5 2 PM10 3 O3 4 NO2 5 SO2 6 CO 7 温度 8 露点 9 气压 10 湿度 11 风速
下载: 导出CSV
表 4 PM2.5预测结果
Table 4 The prediction results of PM2.5
方法 因变量(编号) RMSE MAPE SMAPE - 全部 41.822 8 0.111 1 0.106 1 mBTS-CGCI 4、9、10、11 35.861 8 0.093 0 0.089 3 Lasso-GC 2、3、7、8、10 45.114 1 0.083 6 0.080 0 KGC 2、3、4、9、10 31.858 0 0.074 6 0.072 3 PMIME 11 34.456 5 0.109 3 0.104 8 HSIC-GL-GC 2、3、5、6、11 26.725 3 0.065 3 0.063 8
下载: 导出CSV
-
[1] 谢婷玉, 徐德刚, 阳春华, 等. 基于重要点双重评价的时间序列趋势提取[J]. 信息与控制, 2018, 47(6): 730-736. doi: 10.13976/j.cnki.xk.2018.7553 Xie T Y, Xu D G, Yang C H, et al. Trend feature extraction method for time series based on double evaluation factors of important points[J]. Information and Control, 2018, 47(6): 730-736. doi: 10.13976/j.cnki.xk.2018.7553
[2] Han M, Ren W, Xu M, et al. Nonuniform state space reconstruction for multivariate chaotic time series[J]. IEEE Transactions on Cybernetics, 2018, 49(5): 1885-1895. http://ieeexplore.ieee.org/document/8334805
[3] 李琨, 韩莹, 黄海礁. 基于IBH-LSSVM的混沌时间序列预测及其在抽油井动液面短期预测中的应用[J]. 信息与控制, 2016, 45(2): 241-247. doi: 10.13976/j.cnki.xk.2016.0241 Li K, Han Y, Huang H J. Chaotic time series prediction based on IBH-LSSVM and its application to short-term prediction of dynamic fluid level in oil wells[J]. Information and Control, 2016, 45(2): 241-247. doi: 10.13976/j.cnki.xk.2016.0241
[4] 韩敏, 张瑞全, 许美玲. 一种基于改进灰色关联分析的变量选择算法[J]. 控制与决策, 2017, 32(9): 1647-1652. https://www.cnki.com.cn/Article/CJFDTOTAL-KZYC201709015.htm Han M, Zhang R Q, Xu M L. A variable selection algorithm based on improved grey relational analysis[J]. Control and Decision, 2017, 32(9): 1647-1652. https://www.cnki.com.cn/Article/CJFDTOTAL-KZYC201709015.htm
[5] Granger C W J. Investigating causal relations by econometric models and cross-spectral methods[J]. Econometrica: Journal of the Econometric Society, 1969, 37(3): 424-438. doi: 10.2307/1912791
[6] Bressler S L, Seth A K. Wiener-Granger causality: A well established methodology[J]. Neuroimage, 2011, 58(2): 323-329. doi: 10.1016/j.neuroimage.2010.02.059
[7] Barnett L, Seth A K. The MVGC multivariate Granger causality toolbox: A new approach to Granger-causal inference[J]. Journal of Neuroscience Methods, 2014, 223: 50-68. doi: 10.1016/j.jneumeth.2013.10.018
[8] Geweke J. Measurement of linear dependence and feedback between multiple time series[J]. Journal of the American Statistical Association, 1982, 77(378): 304-313. doi: 10.1080/01621459.1982.10477803
[9] Siggiridou E, Kugiumtzis D. Granger causality in multivariate time series using a time-ordered restricted vector autoregressive model[J]. IEEE Transactions on Signal Processing, 2015, 64(7): 1759-1773. doi: 10.1109/TSP.2015.2500893
[10] Arnold A, Liu Y, Abe N. Temporal causal modeling with graphical granger methods[C]//13th ACM SIGKDD international conference on Knowledge discovery and data mining. New York, USA: ACM, 2007: 66-75.
[11] Bolstad A, Van Veen B D, Nowak R. Causal network inference via group sparse regularization[J]. IEEE Transactions on Signal Processing, 2011, 59(6): 2628-2641. doi: 10.1109/TSP.2011.2129515
[12] Ancona N, Marinazzo D, Stramaglia S. Radial basis function approach to nonlinear Granger causality of time series[J]. Physical Review E, 2004. DOI:10.1103/PhysRevE. 70.056221.
[13] Marinazz D, Pellicoro M, Stramaglia S. Kernel method for nonlinear Granger causality[J]. Physical Review Letters, 2008. DOI: 10.1103/PhysReVLett.100.144103.
[14] Hu M, Liang H. A copula approach to assessing Granger causality[J]. NeuroImage, 2014, 100: 125-134. doi: 10.1016/j.neuroimage.2014.06.013
[15] Montalto A, Stramaglia S, Faes L, et al. Neural networks with non-uniform embedding and explicit validation phase to assess Granger causality[J]. Neural Networks, 2015, 71: 159-171. doi: 10.1016/j.neunet.2015.08.003
[16] Faes L, Nollo G, Porta A. Information-based detection of nonlinear Granger causality in multivariate processes via a nonuniform embedding technique[J]. Physical Review E, 2011. DOI: 10.1103/PhysRevE.83.051112.
[17] Kugiumtzis D. Direct-coupling information measure from nonuniform embedding[J]. Physical Review E, 2013. DOI: 10.1103/PhysRevE.87.062918.
[18] Yamada M, Jitkrittum W, Sigal L, et al. High-dimensional feature selection by feature-wise kernelized Lasso[J]. Neural Computation, 2014, 26(1): 185-207. doi: 10.1162/NECO_a_00537
[19] Poignard B, Yamada M. Sparse Hilbert-Schmidt independence criterion regression[C]//International Conference on Artificial Intelligence and Statistics. Piscataway, USA: IEEE, 2020: 538-548.
[20] Gretton A, Bousquet O, Smola A, et al. Measuring statistical dependence with Hilbert-Schmidt norms[C]//International Conference on Algorithmic Learning Theory. Berlin, Germany: Springer, 2005: 63-77.
[21] Yuan M, Lin Y. Model selection and estimation in regression with grouped variables[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2006, 68(1): 49-67. doi: 10.1111/j.1467-9868.2005.00532.x
[22] Paparoditis E, Politis D N. The asymptotic size and power of the augmented Dickey-Fuller test for a unit root[J]. Econometric Reviews, 2018, 37(9): 955-973. doi: 10.1080/00927872.2016.1178887
[23] Barber R F, Drton M. High-dimensional Ising model selection with Bayesian information criteria[J]. Electronic Journal of Statistics, 2015, 9(1): 567-607. http://arxiv.org/abs/1403.3374
[24] He D, Rish I, Parida L. Transductive HSIC Lasso[C]//SIAM International Conference on Data Mining. Piscataway, USA: IEEE, 2014: 154-162.
[25] Jia Z Y, Lin Y, Liu Y X, et al. Refined nonuniform embedding for coupling detection in multivariate time series[J]. Physical Review E, 2020. DOI: 10.1103/PhysRevE.92.062113.
[26] Papana A, Kyrtsou C, Kugiumtzis D, et al. Simulation study of direct causality measures in multivariate time series[J]. Entropy, 2013, 15(7): 2635-2661. http://www.oalib.com/paper/3095864
[27] 李晓华, 李军. 基于ESN网络的连续搅拌反应釜(CSTR)辨识[J]. 信息与控制, 2014, 43(2): 223-228. doi: 10.3724/SP.J.1219.2014.00223 Li X H, Li J. Identification of continuous stirred tank reactor based on echo state network[J]. Information and Control, 2014, 43(2): 223-228. doi: 10.3724/SP.J.1219.2014.00223
-
期刊类型引用(3)
1. 傅晨，王开军. 决定系数与相关系数辅助的LASSO回归. 福建师范大学学报(自然科学版). 2024(02): 57-63+89 . 
百度学术
2. 李思宇，谭飞. 改进差分算法的仿人控制器设计与参数优化. 组合机床与自动化加工技术. 2023(02): 105-108+112 . 百度学术
3. 乔百友，郝元卿，唐忠，汪锐. 一种NMI结合HSIC0的台风移动轨迹多因素相关分析方法. 东北大学学报(自然科学版). 2023(09): 1234-1244 . 百度学术
其他类型引用(1)
计量
- 文章访问数: 200
- HTML全文浏览量: 47
- PDF下载量: 273
- 被引次数: 4
