一种多维时间序列的定量因果关系判定方法与流程
本发明属于大数据技术处理技术领域，具体的说是一种多维时间序列的定量因果关系判定方法。
背景技术：
因果分析是科学研究的核心问题，是相当多科学研究的直接目的。在致j.s.switzer的一封信中，爱因斯坦曾说过：“西方科学的发展决定于两个伟大成就：一是希腊哲学家发明形式逻辑系统，二是文艺复兴时期发现能通过系统实验推断事件间的因果关系”。现在假设我们针对两个动力事件做实验得到两条时间序列，那么一个最基本的问题来了：我们能不能就凭这两条序列，严格、如实地辨别孰是因孰是果呢？或者，如果它们之间存在“先有鸡还是先有蛋”那样的纠结，我们能否定量地刻画这其间的循环因果呢？
这是一个古老的、著名的跨学科难题，包括诺贝尔奖获得者格兰杰(clivegranger)在内的科学家们已经努力了快半个世纪。事实上，在当今的新兴学科--大数据科学--中它被列为最大的挑战之一(o’neil&schutt,2013,p.274)。近十余年来，liang&kleeman(2005)、liang(2008)、liang(2014)、liang(2016)等通过各种努力，已证明这种因果关系可以严格地由原始物理规律导出，而不必如传统方法那样以半经验的形式出现。传统的因果分析(如granger因果检验以及传递熵分析)在很多情况下验证不了这样一个被称为“零因果准则”的观测事实：如果一个事件的演化不依赖于另一个事件，则后者不是前者的因，但在这个新的体系中，这个事实作为一个已被证明的定理出现。
这种因果分析方法大致可以简述如下，对于一个m维随机动力系统
其中x＝(x1,x2,...,xm)是状态矢量，f是m维矢量场，
其中
以上理论被证明对任意动力系统都严格成立，这与现有的格兰杰因果检验或传递熵分析等经验性的因果检验形成鲜明的对比。但是(ii)式非常复杂，导致它虽然在经典问题中精确成立，但对于实际问题中几乎不能求解，极大地约束了其应用价值。为了解决这个问题，liang(2014)提出一种估计方法，并对二维的问题得到一个相当简洁的公式。该式已经在众多现行的半经验因果分析(如granger因果检验、传递熵分析)解决不了的问题中得到验证，并成功地用到越来越多的现实问题中。例如，仅凭从yahoo上下载的两条ibm与ge公司的股票价格时间序列，就发现70年代有一段时间里由ibm至ge有一个很强的、近乎单向的因果关系，并由此发掘了一个尘封了几十年的关于“七个侏儒”与一个“巨人”争夺巨型机市场的故事。在另一个例子中，stipsetal.(2016)应用此式发现二氧化碳与全球变暖的有着明确的、几乎单向的因果关系。对于最近一百多年来说，co2确实导致了全球变暖，但在一千年以上的古气候尺度上，这个因果关系可能完全颠倒过来，是全球变暖导致了co2浓度的升高。
liang(2014)对于二维liang-kleeman信息流[即m＝2时的(ii)式]的估计是非常成功的，但其成功只限于二维的情况(即m＝2的情况)，对多维情况则无能为力。
技术实现要素：
本发明要解决的技术问题是提供一种多维时间序列的定量因果关系判定方法，该方法能够定量的得到多个事件间的因果关系。
为解决上述技术问题，本发明采用的技术方案为：
一种多维时间序列的定量因果关系判定方法，其特征是：能够定量地估算出多维时间序列间的因果关系的大小与方向，具体实施方案如下：
步骤1，采集m条时间序列，每条序列对应一个事件，m条时间序列分别为x1,x2,…,xm；
步骤2，对每条时间序列进行预处理，得到具有相同时间间隔的多条时间序列；
步骤3，当步骤2中的某些时间序列不平稳时，需对这些序列进行进一步处理，使之变得近似平稳；当步骤2中的时间序列平稳时，进行步骤4；
步骤4，计算各个时间序列之间定量、定向的因果关系；具体计算方法如下：
给定m条平稳、等距的时间序列x1,x2,…,xm，记由xj到xi(i,j＝1,2,..,m,i≠j)的信息流为tj→i，则xj和xi之间的因果关系可由根据tj→i与ti→j来进行判断，tj→i可由下式估算：
其中，
理论上当
步骤5，对步骤4中的计算得到的各两两时间序列之间的因果关系数据进行整合，得到完整的因果关系网络图。
在步骤3中得到的数据平稳的时间序列中选取时间点tn，以时间点tn为中心取时间窗口[tn-l/2,tn+l/2]，在时间长度为l的时间窗口条件下执行步骤4和步骤5，得到时间点tn处的因果关系网络图；
按照时间序列中的时间顺序顺次取不同位置的时间点，重复执行步骤4和步骤5,得到随时间变化的因果关系网络图。
所述的步骤4中所得的
将xi和xj排列至第一、二的位置，此时
其中，x＝(x1,x2,...,xm)是状态矢量，在这里即表现为m条时间序列，
其中
记n时间点的状态矢量x＝(x1,x2,...,xm)为x(n)，n＝1,2,...,n,假设已知x(n)＝(x1,n,x2,n,...,xm,n)的情况下，下一步的状态矢量x(n+1)＝(x1,n+1,x2,n+1,...,xm,n+1)的概率密度，即转移概率密度为ρ，则fisher信息矩阵乘上时间总步数n为
可以证明上式为：
此矩阵中a1j,f1,b1用其估计值
对所得的矩阵(ni)求逆，则所得逆阵的第四个对角分量
给定置信水平α，查标准正态分布表得置信区间系数zα,进一步求得
当该置信区间不包括零在内，则x2是x1的因；反之则因果性不显著。对于任意序列xi与xj之间的因果关系的显著性检验，则只需把xi与xj放置到第一、二的位置，重复以上步骤即可。
当置信水平设置为90％时，置信区间系数zα＝1.65，此时置信区间为：
当置信水平设置为95％时，置信区间系数z＝1.96，此时置信区间为：
当置信水平设置为99％时，置信区间系数z＝2.56，此时置信区间为：
当xi和xj之间的因果关系不显著时，不能简单断定它们之间不存在因果关系，在资料、数据允许的情况下，应扩大数据样本，重复执行步骤1到步骤4。
该种多维时间序列的定量因果关系判定方法能够产生的有益效果为：该方法能够通过简单的计算即可推断多个事件之间的因果关系，而且这种因果关系以定量的形式输出，从而解决了因果关系推断这一难题。一般来说科学研究与工程应用中的因果推断依赖于大量的物理、化学、生物等实验或者计算机仿真模拟，耗费巨大而结果却不见得就正确，而相较之下本方法的计算成本几乎可以忽略不计。该方法在诸如金融分析、神经科学研究、气候变化研究、天气预报、地震预测、雾霾溯源等等各个方面都将有着广泛的应用。
附图说明
图1为本发明一种多维时间序列的定量因果关系判定方法的工作流程图。
图2为利用本发明方法通过股票价格推断公司之间随时间变化的因果关系示意图。
图3为利用本发明方法推断污染物源头时所用的数据采集的地理位置图。
图4为利用本发明方法推断污染物源头时各地的不同污染物浓度对应的时间序列。
图5为利用本发明方法推断污染物源头的工作流程图。
图6为利用本发明方法进行雾霾源头溯源时污染物浓度的空间分布图。
图7为利用本发明方法进行雾霾源头溯源时的工作流程图。
具体实施方式
以下结合说明书附图和具体优选的实施例对本发明作进一步描述。
一种多维时间序列的定量因果关系判定方法，其特征是：能够定量地估算出多维时间序列间的因果关系的大小与方向，具体实施方案如下：
步骤1，采集m条时间序列，每条序列对应一个事件，m条时间序列分别为x1,x2,…,xm；
步骤2，对每条时间序列进行预处理，得到具有相同时间间隔的多条时间序列；
步骤3，当步骤2中的某些时间序列不平稳时，需对这些序列进行进一步处理，使之变得近似平稳；当步骤2中的时间序列平稳时，进行步骤4；
步骤4，计算各个时间序列之间定量、定向的因果关系；具体计算方法如下：
给定m条平稳、等距的时间序列x1,x2,…,xm，记由xj到xi(i,j＝1,2,..,m,i≠j)的信息流为tj→i，则xj和xi之间的因果关系可由根据tj→i与ti→j来进行判断，tj→i可由下式估算：
其中，
理论上当
步骤5，对步骤4中的计算得到的各两两时间序列之间的因果关系数据进行整合，得到完整的因果关系网络图。
本实施例中，步骤1中采集得到的m条时间序列数量大于或等于两条，任意两条时间序列之间因果关系可能存在的情况为：相互不具有因果关系、单向具有因果关系、双向具有因果关系。
本实施例中，步骤2中对每条时间序列进行预处理的方式包括但不限于补充缺测点数据、将不规则时间分布数据映射到规则时间点上等等，最终得到数据映射精准的、具有相同时间长度的多条时间序列。
本实施例中，步骤3中当数据不平稳时对时间序列进行的预处理方式包括但不限于对序列做去除线性趋势处理、或者变换到另一序列，如将股票价格序列换算成收益率序列等。预处理后的时间序列可通过多种经典平稳判据均进行判定，本实施例中也可根据目测大致达到平稳即满足要求。
本实施例中，步骤4中计算各个时间序列之间定量、定向的因果关系，xj和xi之间的因果关系可由根据tj→i与ti→j来进行判断：如果tj→i≠0，则xj是xi的因，若为零则不是；如果ti→j≠0，则xi是xj的因，若为零则不是因。这里总共可有四种可能：(1)xj是xi的因，xi也是xj的因；(2)xj是xi的因，xi不是xj的因；(3)xj不是xi的因，xi是xj的因；(4)xj与xi之间没有因果关系。再计算令i,j＝1,2,...,m,从而得到所有时间序列之间的信息流与因果关系。
本实施例中，对步骤4中得到的每一个信息流估计数据tj→i进行显著性检验，用以判定它是否显著地异于零、抑或不显著地异于零，当显著性检验结果显著地异于零时，则证明两事件之间由因果关系，当显著性检验结果不显著地异于零时，则证明两事件之间的因果关系不显著，所述的显著性检验的具体流程如下：
将xi和xj排列至第一、二的位置，此时
其中，x＝(x1,x2,...,xm)是状态矢量，在这里即表现为m条时间序列，
其中
记n时间点的状态矢量x＝(x1,x2,...,xm)为x(n)，n＝1,2,...,n,假设已知x(n)＝(x1,n,x2,n,...,xm,n)的情况下，下一步的状态矢量x(n+1)＝(x1,n+1,x2,n+1,...,xm,n+1)的概率密度，即转移概率密度为ρ，则fisher信息矩阵乘上时间总步数n为
可以证明上式为：
此矩阵中a1j,f1,b1用其估计值
对所得的矩阵(ni)求逆，则所得逆阵的第四个对角分量
给定置信水平α，查标准正态分布表得置信区间系数zα,进一步求
当该置信区间不包括零在内，则x2是x1的因；反之则因果性不显著。对于任意序列xi与xj之间的因果关系的显著性检验，则只需把xi与xj放置到第一、二的位置，重复以上步骤即可。
进一步，当置信水平设置为90％时，置信区间系数zα＝1.65，此时置信区间为：
当置信水平设置为95％时，置信区间系数z＝1.96，此时置信区间为：
当置信水平设置为99％时，置信区间系数z＝2.56，此时置信区间为：
本实施例中，当步骤4中的计算结果出现xi和xj之间的因果关系不显著时，不能简单断定它们之间不存在因果关系，在资料、数据允许的情况下，应扩大数据样本，重复执行步骤1到步骤4。
进一步，扩大数据样本包括但不限于增加序列长度、增大数据时间分辨率、或者同时增加序列长度和增大数据时间分辨率。
本实施例中，步骤5中可以应用于任选的长度为l的时间窗口，也就是说可以得到该窗口内多条时间序列之间的因果关系网，进一步可以得到随时间变化的因果关系网络，实现方式如下：
在步骤3中得到的平稳时间序列中选取时间点tn，以时间点tn为中心取时间窗口[tn-l/2,tn+l/2]，在时间长度为l的时间窗口条件下执行步骤4和步骤5，得到时间点tn处的因果关系网络图；
按照时间序列中的时间顺序顺次取不同位置的时间点，重复执行步骤4和步骤5,得到随时间变化的因果关系网络图。具体实现步骤如下：
确定估计因果关系的时间窗口长度为l，以及起始时刻t0＝l/2+1，时间序列总长度为tn，计算[t0-l/2,t0+l/2]内的各个变量之间的因果性关系并进行显著性检验，以t0+l/2是否大于或等于tn作为判据，逐次计算在t0＝t0+1时的因果性关系数据，最终得到随时间变化的因果网络图。
进一步举出具体示例对上述技术方案进行验证：
具体实施例一：
从公司股票价格推断两公司之间的关系：
获取亚马逊(amazon)、谷歌(google)、沃尔玛(walmart)股票价格的时间序列，时间从google2004年8月19日上市开始，到2014年12月26日截止，去除交易日、节假日，总计有2607个数据点，由于价格序列不平稳，对价格数据做预处理后取其对数日收益，形成时间序列，预处理的方式具体为：如果p(n)表示第n天的收盘价，则r(n)＝lnp(n+1)-lnp(n)为对数日收益，再以1000交易日为滑动时间窗口长度求得各时间点上两者之间的因果关系，可获得如图2所示的结果，从图中可以看到在90％的概率上由amazon到google的信息流不显著，即amazon从来不是google的因，但google在2010年秋天之前是amazon的因，信息流显著不为零。其原因是因为amazon此前依赖google这个搜索引擎，为买卖提供可能，原理为网上购物依赖搜索引擎。但是2010年秋天之后，google由于违反中国法律被查封，amazon为了不失去中国这个巨大市场而及时弃用google，所以从那以后，google不再是amazon的因。
具体实施例二：
通过不同地区位置的污染物浓度推测污染物源头：
江苏经济发达，工农业排放导致其沿岸海水严重富营养化，使得各种藻类爆发，尤其是盐城外海污染问题特别严重。由于海洋污染原因复杂，某地的污染不见得是当地的排放之故，不由分说随意关停工厂可能导致巨大的经济损失。为了解江苏外海污染的原因，这里用本发明的方法对某些选定点的污染物指标序列做因果分析，计算过程如图5所示。第一步，如图3所示，选取三个点a:(120.5e,34.8n)、b:(121.5e,33.5n)、c:(122.5e,31.8n)，分别对应连云港、盐城、上海外海，取其近二十年污染物浓度数据，时间跨度为1998年初至2016年底，每一步输出间距为15天，总共456个时间点，所得的营养盐(no3、sio3、po4)与叶绿素的时间序列如见图4中所示。下面是a、b、c三点间各种营养盐与叶绿素的因果关系用步骤1至步骤5计算的信息流及其90％置信区间：
硝酸盐(no3)
tb→a＝-0.038621±0.040181
ta→b＝0.204099±0.046894
tc→a＝-0.049355±0.039681
ta→c＝0.335230±0.048638
tc→b＝0.124917±0.051953
tb→c＝0.019727±0.054564
硅酸盐(sio3)
tb→a＝0.092092±0.030884
ta→b＝-0.043766±0.019159
tc→a＝0.026947±0.024741
ta→c＝-0.030499±0.017957
tc→b＝0.258933±0.042229
tb→c＝0.059791±0.049407
磷酸盐(po4)
tb→a＝0.110934±0.055773
ta→b＝0.208629±0.073974
tc→a＝0.027468±0.053460
ta→c＝0.322688±0.070432
tc→b＝0.108395±0.066720
tb→c＝0.077843±0.066274
叶绿素-a(chl-a)
tb→a＝-0.016813±0.044458
ta→b＝0.120196±0.049673
tc→a＝0.176085±0.051851
ta→c＝0.207379±0.053923
tc→b＝0.186485±0.054840
tb→c＝0.041304±0.051044
对于这上述几项指标来说，可以看出，北边、南边的海域都是b地外海的因，所以b地外海的污染不只是b地海域自己的结果，它还来自南北的邻居。同样，就硅酸盐、磷酸盐而言，b地外海确实是a地外海、c地的因，因为信息流tb→a，tb→c显著地异于零。
但是，非常意外的是，对硝酸盐、叶绿素来说，tb→a，tb→c并不显著地异于零，也就是说，就这个数据样本来看，b地外海似乎不是上海与a地外海的因，因为信息流不显著为零。如果要讨论的是源自硝酸盐与叶绿素的污染(如藻华)，b地外海虽然污染严重，虽然长期以来认为苏北的污染通过江苏沿岸流会污染长江口，但它不是上海与a地外海的因。
具体实施例三：
通过不同地区位置大气中的pm2.5浓度进行雾霾溯源：
如图6、图7所示，通过获取多地区pm2.5浓度监测站的某个时间段内的大气中的pm2.5浓度数据，形成相应的时间序列，可通过计算得出多地点之间的雾霾因果网络图，进而找到一个或多个主源头。
由上述多个具体实施例和计算原理可以进一步说明本方法可以通过简单的计算即可推断多个事件之间的因果关系，而且这种因果关系以定量的形式输出，从而解决了因果关系推断这一难题。一般来说科学研究与工程应用中的因果推断依赖于大量的物理、化学、生物等实验或者计算机仿真模拟，耗费巨大而结果却不见得就正确，而相较之下本方法的计算成本几乎可以忽略不计。该方法在金融、气候、环境、神经科学等众多科学领域都将有着广泛的应用。
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。
- 还没有人留言评论。精彩留言会获得点赞！