CN111367959A - 一种零时滞非线性扩展Granger因果分析方法 - Google Patents

Abstract

一种零时滞非线性扩展Granger因果分析方法，属于数据挖掘技术领域。传统的Granger因果模型基于滞后项的VAR模型，但忽略零时滞会极大地改变滞后项的模型系数，从而导致错误的因果识别。其次，传统的Granger因果模型只能应用于线性系统的因果识别，对于非线性系统可能产生错误的因果识别。基于上述分析，本发明首先对传统的VAR模型进行扩展，设定变量的滞后阶数，采用高斯核函数将原始数据进行非线性映射，然后建立包含零滞后项的结构VAR模型，最后并根据结构VAR模型的残差进行Granger因果关系识别，实现如污染及气象等非线性系统的零时滞因果关系分析。本发明能够克服传统Granger因果模型的不足，实现对非线性系统的扩展Granger因果关系分析。

Description

一种零时滞非线性扩展Granger因果分析方法

技术领域

本发明属于数据挖掘技术领域，涉及一种零时滞非线性扩展Granger因果分析方法，旨 在挖掘气象及环境等非线性系统中各时间序列变量间的(零)时滞因果影响关系。

背景技术

多元时间序列是指按照时间顺序排列的一系列数字序列的集合，广泛存在于工业、医学、 金融和气象等多个领域。例如，大气污染研究领域。过去几十年，我国由于工业化的快速发 展及汽车的推广普及，城市空气往往会受到燃煤、工业废气、废料等排放物及汽车尾气的污 染，导致空气污染严重。尤其是我国北部和东部许多地区都出现了不同程度的雾霾天气。雾 霾主要由CO、SO₂、NO_x和可吸入颗粒物(主要是直径小于2.5微米的颗粒物，称为PM2.5) 组成，其中PM2.5是加重雾霾的主要原因。PM2.5等空气污染物不仅会对经济发展产生负面 影响，造成交通拥堵等社会问题，还会危害人类健康，严重时可能会损害人体心血管和呼吸 系统，增加罹患癌症的风险。因此，针对PM2.5等大气污染物进行因果关系分析，挖掘污染 物(AQI)及气象系统潜在的影响关系，并以此做出污染治理的针对性防治措施以及对污染 情况的预报，从而为大气污染治理提供理论支持，具有重要的现实意义。

向量自回归模型(VAR)模型可以用于描述动力学系统变量间的相互作用关系，即给定 一个长度为N的M维时间序列X(t)＝{X₁(t),X₂(t),…,X_M(t)},t＝1,2…,N，它们的动力学相 互作用关系可以由p阶VAR模型描述，具体为：

式中，

是VAR模型系数，p为模型阶数，k为滞后阶数，t为采样时间。A_ji(k)描 述了X_j(t)与X_i(t-k),(i,j＝1,2,…,M,k＝1,2,…,p)的依赖关系；E＝[E₁,E₂…,E_M]代表模型 残差。

因果关系是一个系统(因)对另一个系统(果)的作用关系。Granger因果模型是一种常 用的基于VAR模型的因果关系分析方法。对于两个时间序列X和Y，如果X的历史信息的加 入有助于减少仅由Y的历史信息预测Y未来时刻的预测误差，则可以说存在X→Y的Granger 因果关系。Granger因果模型一经提出，就受到了国内外学者的广泛关注，提出了大量改进模 型。但现有的绝大部分Granger因果关系分析模型侧重于多变量、非线性和时变等特性研究， 如条件Granger因果模型、核Granger因果模型和Lasso-Granger因果模型等。然而，基于传 统VAR模型的Granger因果方法依赖于模型系数来进行因果解释，但这是基于下述假设：样 本值的时间分辨率高于因果影响的时间尺度。当不满足上述假设时，瞬时效应的存在可能会 导致传统VAR模型对因果关系产生错误的解释。

等人在论文“

A,Zhang K,Shimizu S,et al.Estimation of a structural vectorautoregression model using non-gaussianity[J]. Journal of Machine LearningResearch,2010,11(May):1709-1731.”中表示如果忽略瞬时影响， 会对基于VAR模型的因果关系分析结果产生不利影响，并进行了理论证明，表明模型中零滞 后相关性的遗漏会极大地改变滞后项的模型系数，从而导致错误的因果识别。

在VAR模型中，滞后时间表示观测样本值采样间隔的整数倍，零时滞表示同一时刻的观 测样本值。如果由于采集手段的不成熟，所能采集到的观测样本时间间隔较长，相邻两个样 本之间则会呈现出明显的长时滞现象。然而，实际系统变量间可能具有短时甚至零时滞相互 作用关系。原有的基于非零时滞的VAR模型可能产生错误的结果。综上所述，针对时间序列 数据，瞬时影响关系取决于时间序列的采样间隔。所以研究时间序列数据的零时滞相关关系 具有重要的研究意义。此外，传统的GC模型基于线性VAR模型，面对非线性系统可能识别 出虚假因果关系。而实际系统往往具有复杂的非线性结构，因此传统GC模型对实际系统进 行因果关系分析时，往往面临着错误因果关系识别的风险。

本发明针对非线性系统的(零)时滞因果关系分析问题，利用核函数将扩展Granger因 果模型推广到非线性系统，提出一种基于结构VAR模型的非线性扩展Granger因果分析方法， 用于气象污染等非线性系统变量间的瞬时因果影响分析。

发明内容

本发明要解决的技术问题是，针对传统GC方法没有考虑零时滞对于滞后项的模型系数 的及因果识别准确性的影响，对传统GC模型进行非线性扩展，提出一种基于扩展Granger 因果模型的非线性因果模型，实现对实际的非线性系统进行零时滞因果关系识别。本方法的 目的在于挖掘非线性系统包含(零)时滞的因果关系。

本发明面向AQI及气象等实际非线性系统，提出了挖掘系统中其他变量与主要污染物 PM2.5之间的Granger因果关系方法。先将采集得到的时序数据进行预处理，对缺失值和异 常值进行填补。然后将数据进行平稳化及去平稳化处理，以满足Granger因果建模的需要。 之后对数据进行归一化处理，以消除不同变量量纲对于分析结果的影响。最后建立包含零时 滞项的结构VAR模型，得到不同特征之间的Granger因果指数，从而定性、定量地分析系统 中其他变量与PM2.5的Granger因果关系，实现对PM2.5影响因素的分析研究。

为了实现上述目的，本发明采用的技术方案如下：

一种零时滞非线性扩展Granger因果分析方法，所述的分析方法对传统的VAR模型进行 了扩展，首先设定变量的滞后阶数，然后采用高斯核函数对时序数据进行非线性扩展，最后 构建包含零时滞项的多变量结构VAR模型，并根据结构VAR模型的残差进行Granger因果 关系识别，实现如污染及气象等非线性系统的零时滞因果关系分析。对于一个时间序列X(t)， 采用如下公式：

式中，p表示设定的VAR模型阶数，A是模型系数，Φ(·)表示进行非线性映射的高斯核函数， E为模型的拟合残差，k表示滞后阶数，t为采样时刻。

传统结构VAR模型只包含了X(t)非零滞后项，本发明充分考虑零滞后项对于VAR模型 系数的影响，以提高Granger因果模型的因果识别准确性。同时，利用高斯核函数对原始线 性系统进行非线性扩展，克服其可能对于非线性系统产生错误识别的缺点。通过公式(1)不 光可以有效识别因果关系的方向性特征，还能进行定量识别因果关系。从而达到变量选择的 效果，得到一个可以有效降低预测精度、精简预测模型的变量子集。实现对污染物及气象等 系统的时间序列的有效分析。

给定一个长度为N的M维AQI时间序列X(t)＝{X₁(t),X₂(t),…,X_M(t)},t＝1,2…,N， 其中，X₁(t),X₂(t),…,X_M(t)表示在t时刻采集的M维AQI时间序列，t表示数据的采样时刻。通过构建基于非线性结构VAR模型分析是否存在X_j→X_i的非线性零时滞因果关系，实 现对AQI及气象系统的因果关系分析，具体步骤如下：

步骤1：获取AQI及气象系统的观测站采集数据，并对多维AQI及气象系统数据集进行 预处理，对缺失值进行三次样条插值法进行缺失值补充和离群点数据替换，得到平滑的时间 序列数据，为后续步骤做准备。

步骤2：采用单位根检验方法对时间序列数据进行平稳化检验及差分平稳化处理，具体 是：如果时间序列具有单位根，表明该时间序列是非平稳时间序列，采用差分平稳化方法进 行平稳化处理；如果没有单位根，则为平稳时间序列，不需进行平稳化处理。

步骤3：对时间序列数据进行归一化处理，消除多维变量不同度量单位对分析结果的影 响。

步骤4：设定时间序列数据的最大滞后阶数p，对时间序列数据进行滞后项扩展。

步骤5：选取合适的高斯核函数，对原始数据进行非线性扩展，然后建立如下式(3)(4)所 示的非线性结构VAR模型。

式中，公式(3)是非限制性结构VAR模型，公式(4)表示不包含X_j的限制性结构VAR模型。p 表示设定的模型阶数；X^-j＝[X₁…X_j-1X_j+1…X_M]^T表示不包含第j维变量X_j的时间序列； A_i(k),

分别表示VAR模型的模型系数；E_i,

表示模型残差；Φ(·)表示对数据进行非 线性映射的核函数。

选择高斯核函数对原始数据进行非线性扩展，其核函数如下式所示：

式中，

表示X_i的中心，

表示X_i的协方差，X_i表示原始时间序列的的第i维变量X_i。

步骤6：对公式(3)、公式(4)进行变形求解，获取因果关系方向并计算模型残差。具体为：

步骤6.1：获取时间序列间因果关系的方向；

采用非高斯数据的两两似然比，评估与残差相关的序列对之间的零滞后相互作用的方向。 对于两个标准化的数据X_i,X_j，定义两个回归模型X_j＝AX_i+e,(方向为X_i→X_j)和 X_i＝AX_j+d,(方向为X_j→X_i)，在n次试验中计算两个模型出现的概率比值，如下式所示：

如果R＞0则表示因果方向为X_i→X_j，反之R＜0则表示因果方向为X_j→X_i。

步骤6.2：根据最小二乘方法求解VAR模型，计算模型残差：对于传统的AR模型，借助最小二乘方法实现对模型残差E(t)的估计。即对于公式(1)，得：

所以，由公式(3)、公式(4)可得，模型残差的估计值为：

式中，

和

为根据最小二乘法求解得到的模型系数估计值。

步骤7：然后，根据Granger因果关系的定义，计算得到非线性扩展Granger因果关系指 数F_j→i，如下式所示：

式中，var表示变量方差，

和

分别为公式(3)、公式(4)的模型残差估计值。

步骤8：针对一个变量，获取不同变量对该变量的因果关系后，根据因果关系系数选择 与目标变量相关的变量子集，建立回声状态网络进行建模预测，实现对目标变量的分析预测 研究，验证因果关系分析的分析结果。

与现有技术相比，本发明具有以下明显优势：

本发明着眼于非线性系统的零时滞和非零时滞Granger因果关系分析。首先，本发明采 用高斯核函数将原始数据进行非线性映射，在再生核希尔伯特空间中进行因果关系分析。将 传统的Granger因果模型进行非线性扩展，克服了传统方法只能应用于线性系统的缺点。然 后，本发明通过建立结构VAR模型，将因果关系分析结果与时间序列的采样间隔联系起来， 解决了传统Granger因果方法分析时由于测量值的时间分辨率低于因果影响的时间尺度而导 致的VAR模型的错误因果识别问题。即充分考虑了零时滞在因果关系分析中的重要影响关 系。

附图说明

图1为包含瞬时因果关系的方向示意图(虚线表示瞬时因果关系)；

图2为本发明的流程图；

图3为PMIME对PM2.5的预测误差曲线图。

图4为mBTS-CGCI对PM2.5的预测误差曲线图。

图5为本发明方法对PM2.5的预测误差曲线图。

具体实施方式

以下结合具体实例，并结合仿真图，对本发明做进一步详细说明。

本发明所使用的硬件设备包括PC机一台。

图1为本发明提供的一种零时滞非线性扩展Granger因果分析方法，具体包括以下步骤：

一种零时滞非线性扩展Granger因果分析方法，所述的分析方法首先采用高斯核函数对 数据进行非线性扩展，并采用结构VAR模型分析样本集特征之间的(零)时滞Granger因果 关系。例如，对于一组标准数据集，如下式所示：

式中，t表示采样时刻，X₁(t)，X₂(t)，X₃(t)表示3维时间序列数据，δ是零时滞模型系数， X_i(t-1),X_i(t-2)分别表示在t-1时刻和t-2时刻的时间序列数据，a_ij(1),a_ij(2)(i,j＝1,2,3)分 别表示在t-1时刻和t-2时刻变量X_j→X_i(i,j＝1,2,3)的VAR模型系数，e_i为模型的拟合残 差。可以看到，式(11)所示的3维变量间可能存在零时滞影响关系。即当δ≠0表示存在零时 滞因果关系，反之当δ＝0时表示不存在零时滞因果关系。其示意图如图2所示，实线表示非 零时滞因果关系，虚线表示零时滞因果关系，数字表示延迟时间。

给定一个长度为N的M维AQI及气象时间序列 X(t)＝{X₁(t),X₂(t),…,X_M(t)},t＝1,2…,N。

该方法包括以下步骤：

步骤1：获取北京地区2014年3月24日—2015年2月22日共计886组大气AQI及气 象系统的观测站采集数据，数据共11维，每8小时采集一次，其各变量含义编号如表1所示。 对多维AQI及气象系统数据集对缺失值及异常值进行预处理。

表1北京AQI及气象数据编号对照表

步骤2：采用单位根检验方法对时间序列数据进行平稳化检验及差分平稳化处理。

步骤3：对时间序列数据进行归一化处理，消除多维变量不同度量单位对分析结果的影 响。

步骤4：设定时间序列数据的滞后阶数p＝8对时间序列数据进行滞后项扩展。

步骤5：选用核宽为σ＝0.12的高斯核函数进行非线性映射，然后建立如式(3)(4)所示的 非线性结构VAR模型。

步骤6：对式(3)(4)进行变形求解，获取因果关系方向，求解因果模型的模型残差。具体 为：

步骤6.1：获取序列间因果关系的方向，利用Hyvarinen和Smith提出的非高斯数据的两 两似然比，评估与残差相关的序列对之间的零滞后相互作用的方向。对于两个标准化的数据 X_i,X_j，定义两个回归模型X_j＝AX_i+e,(方向为X_i→X_j)和X_i＝AX_j+d,(方向为X_j→X_i)，在n次试验中计算两个模型出现的概率比值，如下式所示：

如果R＞0则表示因果方向为X_i→X_j，反之R＜0则表示因果方向为X_j→X_i。

定性因果关系如表2所示：当表中的因果指数大于0表示表中第一列对应变量对PM2.5 存在因果关系，反之如果小于0则表示PM2.5对对应的第一列变量存在因果关系。可以看到， 在进行定性分析时，北京当地对PM2.5存在因果关系的变量分别是PM10、CO和气温等。

表2针对PM2.5因果指数

步骤6.2：根据最小二乘方法求解VAR模型，求解模型残差。

步骤7：根据Granger因果关系的定义求得因果系数，得到如表2所示的AQI及气象系统 变量间的非线性扩展Granger因果关系指数，对因果关系进行定量分析，并将结果与基于混 合嵌入的偏互信息(PMIME)和基于后向选择的条件Granger因果方法(mBTS-CGCI)的分 析结果进行比较。其中，mBTS-CGCI和本文方法均满足显著性为0.05的显著性检验。从表 中可以看出，本文方法得出PM10、气温和风速对PM2.5存在显著性水平为0.05的显著性因 果关系。PMIME则识别到了PM10对PM2.5的因果关系，mBTS-CGCI的结果为PM10和CO 对PM2.5存在因果关系。

步骤8：根据因果系数选择与目标变量相关的变量子集，建立回声状态网络进行建模预测， 对因果关系分析的结果进行检验。三种方法均设定回声状态网络的储备池维数为125，稀疏 度为0.051，谱半径为0.99，进行30次实验取平均值，以消除偶然误差的影响，预测指标有 均方根误差(RMSE)和对称平均绝对百分比误差(SMAPE)。其定义如下：

式中，x_i和

分别表示真实值和预测值，n表示样本数。这两个指标都是越小表示预测效果 越好。

预测结果如图3、4和5所示。从图中可以看出，采用本发明的方法选择的变量子集对PM2.5 的预测效果图更能准确跟踪PM2.5的变化趋势。30次实验的误差如表3所示，可以看到，本 发明方法的两个预测指标值分别为RMSE＝9.8571，SMAPE＝0.0357均取得了最小值，进一步 验证了本发明方法能够更加准确预测，提高预测精度的优势。

表3各种方法预测结果

以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围 的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做 出若干变形和改进，这些均属于本发明的保护范围。

Claims (2)

1.一种零时滞非线性扩展Granger因果分析方法，其特征在于，所述的分析方法对传统的VAR模型进行了扩展，首先设定变量的滞后阶数，然后采用高斯核函数对时序数据进行非线性扩展，最后构建包含零时滞项的多变量结构VAR模型，并根据结构VAR模型的残差进行Granger因果关系识别，实现如污染及气象等非线性系统的零时滞因果关系分析；

给定一个长度为N的M维AQI时间序列X(t)＝{X₁(t),X₂(t),…,X_M(t)},t＝1,2…,N，其中，X₁(t),X₂(t),…,X_M(t)表示在t时刻采集的M维AQI时间序列，t表示数据的采样时刻；通过构建基于非线性结构VAR模型分析是否存在X_j→X_i的非线性零时滞因果关系，实现对AQI及气象系统的因果关系分析；其中，时间序列X(t)，采用如下公式：

式中，p表示设定的VAR模型阶数，A是模型系数，Φ(·)表示进行非线性映射的高斯核函数，E为模型的拟合残差，k表示滞后阶数，t为采样时刻；

所述因果分析方法具体步骤如下：

步骤1：获取AQI及气象系统的观测站采集数据，并对多维AQI及气象系统数据集进行预处理，对缺失值进行三次样条插值法进行缺失值补充和离群点数据替换，得到平滑的时间序列数据；

步骤2：采用单位根检验方法对时间序列数据进行平稳化检验及差分平稳化处理，具体：如果时间序列具有单位根，表明该时间序列是非平稳时间序列，采用差分平稳化方法进行平稳化处理；如果没有单位根，则为平稳时间序列；

步骤3：对时间序列数据进行归一化处理，消除多维变量不同度量单位对分析结果的影响；

步骤4：设定时间序列数据的最大滞后阶数p，对时间序列数据进行滞后项扩展；

步骤5：选取合适的高斯核函数，对原始数据进行非线性扩展，并建立如下式(3)、(4)所示的非线性结构VAR模型；

式中，公式(3)是非限制性结构VAR模型，公式(4)表示不包含X_j的限制性结构VAR模型；p表示设定的模型阶数；X^-j＝[X₁…X_j-1X_j+1…X_M]^T表示不包含第j维变量X_j的时间序列；A_i(k),

分别表示VAR模型的模型系数；E_i,

表示模型残差；Φ(·)表示对数据进行非线性映射的核函数；

步骤6：对公式(3)、公式(4)进行变形求解，获取因果关系方向并计算模型残差；具体为：

步骤6.1：获取时间序列间因果关系的方向；

采用非高斯数据的两两似然比，评估与残差相关的序列对之间的零滞后相互作用的方向；对于两个标准化的数据X_i,X_j，定义两个回归模型X_j＝AX_i+e,(方向为X_i→X_j)和X_i＝AX_j+d,(方向为X_j→X_i)，在n次试验中计算两个模型出现的概率比值，如下式所示：

如果R＞0则表示因果方向为X_i→X_j，反之R＜0则表示因果方向为X_j→X_i；

步骤6.2：根据最小二乘方法求解VAR模型，计算模型残差：对于传统的AR模型，借助最小二乘方法实现对模型残差E(t)的估计；即对于公式(1)，得：

所以，由公式(3)、公式(4)可得，模型残差的估计值为：

式中，

和

为根据最小二乘法求解得到的模型系数估计值；

步骤7：然后，根据Granger因果关系的定义，计算得到非线性扩展Granger因果关系指数F_j→i，如下式所示：

式中，var表示变量方差，

和

分别为公式(3)、公式(4)的模型残差估计值；

步骤8：针对一个变量，获取不同变量对该变量的因果关系后，根据因果关系系数选择与目标变量相关的变量子集，建立回声状态网络进行建模预测，实现对目标变量的分析预测研究，验证因果关系分析的分析结果。

2.根据权利要求1所述的一种零时滞非线性扩展Granger因果分析方法，其特征在于，所述的步骤5中，选择高斯核函数对原始数据进行非线性扩展，其核函数如下式所示：

式中，

表示X_i的中心，

表示X_i的协方差，X_i表示原始时间序列的的第i维变量X_i。

Images