南京大学学报（自然科学）, 2023, 59(4): 550-560 doi: 10.13232/j.cnki.jnju.2023.04.002

基于分位数因子模型的高维时间序列因果关系分析

梁慧玲1,2, 刘慧,1,2, 刘力维1,2, 赵佳3, 阮怀军3

1.山东财经大学计算机科学与技术学院，济南，250014

2.山东省数字媒体技术重点实验室，山东财经大学，济南，250014

3.山东省农业科学院信息技术研究所，济南，250000

Causal relationship analysis of high⁃dimensional time series based on quantile factor model

Liang Huiling1,2, Liu Hui,1,2, Liu Liwei1,2, Zhao Jia3, Ruan Huaijun3

1.College of Computer Science and Technology，Shandong University of Finance and Economics，Ji'nan，250014，China

2.Key Laboratory of Digital Media Technology of Shandong Province，Shandong University of Finance and Economics, Ji'nan，250014，China

3.Institute of Information Technology，Shandong Academy of Agricultural Sciences，Ji'nan，250000，China

通讯作者: E⁃mail：liuh_lh@sdufe.edu.cn

收稿日期: 2023-06-13  

基金资助: 国家自然科学基金.  62072274
山东省科技成果转移转化项目.  2021LYXZ021
山东省泰山学者特聘专家计划.  tstp20221137

Received: 2023-06-13  

摘要

从观察数据中发现变量之间的因果关系是许多科学研究领域的关键问题，传统Granger因果模型受到维度灾难的影响，难以准确地在高维时间序列中发现因果关系.提出一种基于分位数因子模型的Granger因果分析新方法QFM⁃CGC用于高维时间序列因果关系的判定.首先，QFM⁃CGC采用赤池信息量准则进行模型选择，避免人为干预设置滞后阶数的操作；然后，对向量自回归（Vector Autoregressive，VAR）模型中的条件变量建立分位数因子模型进行降维，减少VAR模型中的待估计系数，对降维后的VAR模型重新进行条件Granger因果分析；最后，使用蒙特卡洛模拟评估不同方法识别底层系统与观测时间序列的连通性结构的能力.在不同维度变量的线性仿真系统和两组现实数据集上与基准方法和经典方法进行了比较，实验结果验证了该方法的有效性.

关键词： 高维时间序列 ; 分位数因子模型 ; 条件Granger因果分析 ; 数据挖掘

Abstract

Finding the causal relationship between variables from observed data is a key issue in many scientific research fields. Because the traditional Granger causality model is affected by the curse of dimension，it is difficult to accurately find causality in high⁃dimensional time series. In this paper,we propose a new Granger causality analysis method based on quantile factor model，QFM⁃CGC algorithm，which is used to find causality relationship in high⁃dimensional time series. Firstly，QFM⁃CGC uses Akaike information criterion to select models，which avoids setting the lag order by human intervention. Then，the quantile factor model is established to reduce the dimensionality of the conditional variables in a vector autoregressive (VAR) model，thus reducing the number of coefficients that need to be estimated. The reduced⁃dimensional VAR model is used for a conditional Granger causality analysis. Finally，Monte Carlo simulation is applied to evaluate the performance of different methods to identify the connectivity structure between the underlying system and the observation time series. Experiments compare the proposed method with benchmark and classical methods on a linear simulation system with variables in different dimensions and two sets of real data，confirming its effectiveness.

Keywords： high⁃dimensional time series ; quantile factor model ; conditional Granger causality analysis ; data mining

PDF (1113KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

梁慧玲, 刘慧, 刘力维, 赵佳, 阮怀军. 基于分位数因子模型的高维时间序列因果关系分析. 南京大学学报（自然科学）[J], 2023, 59(4): 550-560 doi:10.13232/j.cnki.jnju.2023.04.002

Liang Huiling, Liu Hui, Liu Liwei, Zhao Jia, Ruan Huaijun. Causal relationship analysis of high⁃dimensional time series based on quantile factor model. Journal of nanjing University[J], 2023, 59(4): 550-560 doi:10.13232/j.cnki.jnju.2023.04.002

时间序列是指属于同一统计指标的数值按其时间发生的先后顺序排列而形成的一组随机变量，可以分成一元时间序列和多元时间序列，其中多元时间序列是将多个一元时间序列组合形成的时间序列.1969年，英国计量经济学家Granger［1］首次提出一种经济学上的统计学假设检验方法，通常称为Granger因果关系检验，广泛应用于神经科学［2-3］、计量经济学［4-5］等其他研究领域［6］.根据Granger最初提出的概念，如果Y的预测模型中包括X时，Y的预测得到了改善，则变量X会Granger导致变量Y.

多元时间序列在进行传统Granger因果关系分析时忽略了其他变量的存在，因此双变量因果关系测度在估计真实因果关系时不准确［7］.针对这些问题，1982年Geweke［8］提出条件Granger因果分析（Conditional Granger Causality Analyse，CGCA）.多元时间序列的因果关系分析依赖于从一个观测变量到另一个观测变量的直接因果关系的估计，并考虑其他观察变量的存在.相关影响由条件Granger因果关系指数（Conditional Granger Causality Index，CGCI）量化，由线性向量自回归（Vector Autoregressive，VAR）模型推导.

随着信息时代的发展，数据的维度不断增加，对于数据的研究也已转变到高维空间［9-10］，然而，目前大部分因果分析研究仍然集中于二元或多元时间序列，对于高维和超高维的时间序列因果分析缺少有效的处理方法.通常低维主时间序列能提供动态特征的良好表示，易于解释和可视化.因此，对于高维和超高维的时间序列因果分析需要使用降维方法来限制VAR模型.

主成分分析（Principal Component Analysis，PCA）是目前最主要的降维方法之一［11］，它将原始的高维数据投影到一个较低维的子空间上，使原始高维数据可以由一组低维变量表示［12］.但是，PCA在降维时没有捕获隐藏的因子，值得注意的是相关因子可能会改变时间序列的分布特征（矩或分位数），而不是其均值.Chen et al［13］提出分位数因子模型（Quantile Factor Models，QFM）及其估计程序，简称分位数因子分析（Quantile Factor Analysis，QFA）.与其他因子模型不同，QFM还捕获移动可观测分布的其他相关部分的未观测因子.QFM的一个重要优点是它能同时提取决定QFM因素结构的所有均值和额外（非均值）因子，而PCA只能提取平均因子，所以QFA克服了PCA没有捕捉隐藏因子的能力的问题.为此，Chen et al［13］通过蒙特卡洛模拟［14］说明存在异常值时，使用QFA有明显优势［13］.

由于传统Granger因果分析方法无法准确区分直接因果关系和间接因果关系，且高维时间序列在VAR模型中容易受到维度灾难的影响，难以在高维时间序列中准确地发现因果关系，本文将分位数因子模型与条件Granger因果关系分析相结合，提出QFM⁃CGC算法来处理高维数据的因果分析.QFM⁃CGC算法将降维技术运用在降低VAR模型条件项的计算中，减少VAR模型中待估计系数，再对降维后的VAR模型重新进行条件Granger因果分析，避免传统方法受到的维度灾难的影响.在线性仿真系统和现实宏观经济数据集上与四种基准方法和经典方法进行对比实验，结果表明，在两个不同维度的仿真实验中，本文提出的方法识别正确因果关系的概率平均提高6%和3.46%.另外，在现实宏观经济数据的实验中发现，从均方根误差、平均绝对百分误差和对称平均绝对百分比误差三个指标来看，本文提出的方法的因果关系序列的预测效果相较于最优的对比方法，分别降低了10.47%，5.18%和2.09%.最后，本文方法的北京空气质量指标（Air Quality Index，AQI）及气象时间序列的预测效果相较于最优的对比方法，三个误差指标分别降低了11.87%，14.09%和9.79%.

1 相关理论

1.1　条件Granger因果关系指数

为了解决传统的Granger因果模型在多变量系统中容易生成虚假的因果关系的问题，Geweke［8］引入条件变量，提出了条件Granger因果分析方法.

假设Xt=X1,t,X2,t,…,XK,tt=1,…,N是长度为N的K维平稳时间序列.从驱动变量Xi到响应变量Xj的CGCI的定义涉及Xj的两个VAR模型，也称动态回归模型［15］.第一个模型是无限制模型［16］（U⁃模型），表示为：

Xj,t=∑k=1Kajk,1Xk,t-1+⋯+ajk,pXk,t-p+uj,t

其中，p是模型的阶数，ajk,lk=1,⋯,K,l=1,⋯,p是U⁃模型的系数.

第二个模型是从U⁃模型排除Xi的滞后导出的受限模型（R⁃模型），表示为：

Xj,t=∑k=1,k≠iKbjk,1Xk,t-1+⋯+bjk,pXk,t-p+ej,t

其中，bjk,lk=1,…,K但k≠i且l=1,2,⋯,p是R⁃模型的系数.uj,t和ej,t是均值为零，方差分别为σU2和σR2的白噪声.用普通最小二乘法拟合U⁃模型和R⁃模型，得到残差σ^U2和σ^R2的估计值.条件Granger因果关系由CGCI量化，定义为R⁃模型和U⁃模型的误差方差之比的对数［17］，如式（3）所示：

CGCIXi→Xj=lnσ^R2σ^U2

显然，当Xi不改进Xj的预测时，即U⁃模型和R⁃模型给出的拟合误差方差大致相同，CGCI处于零水平；当Xi改进Xj预测时，CGCI获得更大的正值，表明Xi Granger导致Xj.

1.2　分位数因子模型

分位数因子模型（QFM）是代表高维面板数据的一类新的因子模型，旨在提取大型面板数据集分布均值处的公共因子［13］.设Xit是一个由N个观测单位组成的面板，且每个变量都有T个观测值.那么，在τ∈0,1处Xit有如下的因子结构：

QXitτftτ=λi'τftτs.t. i=1,2,⋯,N且t=1,2,⋯,T

其中，公共因子ftτ是rτ×1维不可观测因子向量，λτi是rτ×1维因子载荷向量且rτ≪

N，ftτ和λτi在每个τ上可能不同.为了简化符号，下文中抑制ftτ，λτi和rτ对τ的依赖性，分别改写为ft，λi和r，因此，式（4）中的因子结构可以通过最小化目标函数获得，如式（5）所示：

LNTθ=1NT∑i=1N∑t=1TρτXit-λi'ft

其中，θ=(λ1',⋯,λN',f1',⋯,fT')'为待估参数，ρτu=τ-1u≤0为分位数损失函数.对因子与因子载荷进行以下规范化：

1T∑t=1Tftft'=Ir1N∑i=1Nλiλi'为对角元素非增的对角

Chen et al［13］给出了当样本矩阵Y的维度趋于无穷时因子模型估计量的渐近性质，提出迭代分位数回归（Iterative Quantile Regression，IQR）算法，可以有效地找到目标函数的平稳点.令：

Λ=(λ1,λ2,⋯,λN)'
F=(f1,f2,⋯,fT)'

并定义以下平均值：

Mi,Tλ,F=1T∑i=1TρτXit-λ'ft
MN,tΛ,f=1N∑i=1TρτXit-λi'f

IQR的迭代过程如下：

（1）随机选择起始参数F(0)；

（2）给定Fl-1，对i=1,…,N，计算λil-1=argminλΜi,Tλ,Fl-1；

给定Λl-1，对t=1,…,T，计算ftl-1=argminfΜt,NΛl-1,f；

（3）对于l=1,…,L，迭代第二步直到ΜNTθL接近ΜNTθL-1，其中θl=vechΛl',vechFl''；

（4）规范化Λl和Fl，使其满足（3）的规范化.

通过Chen et al［13］提出的基于秩最小化的方法来确定不可观测因子的数量.

2 基于分位数因子模型的条件Granger因果关系分析

基于分位数因子模型的条件Granger因果关系分析模型包括以下几个部分.

2.1　平稳性分析

由于Granger因果模型的时间序列具有平稳性，即该时间序列均值和方差没有系统的变化且严格消除了周期性变化，因此，在建模前需要对数据进行平稳性分析.本文采用Augmented Dickey⁃Fuller （ADF）［18］检验，确定序列中是否存在单位根，帮助判断该序列是否平稳，如式（7）所示：

ΔXt=α+βt+δXt-1+∑i=1mβiΔXt-i+εt

其中，α是一个常数；β是趋势项系数；δ是自回归系数，描述数据的平稳性；m是自回归模型的阶数；εt是误差项.若序列不存在单位根，表示时间序列是平稳的，可以直接进行Granger因果关系分析；反之，时间序列是非平稳的，需要进行差分，实现时间序列平稳化后再进行因果关系分析.

2.2　模型选择

VAR模型中，如果解释变量的最大滞后阶数p太小，残差可能存在自相关，导致参数估计不一致.虽然适当增加滞后阶数p可以解决此问题，但p过大会使待估计参数增多，严重降低自由度，最终影响模型参数估计的有效性［19］，所以VAR模型中解释变量的最大滞后阶数p的选择很重要.本文采用AIC （Akaike Information Criterion）［20］来自动选择合适的模型阶数，以消除人为选择的不确定性的干扰，如式（8）所示：

AIC=2k-2lnL

其中，L表示似然函数，k是拟合模型中参数的数量.假设模型误差服从独立正态分布，设n为观测值数目，RSS为残差平方和，则式（8）可以改写为：

AIC=2k-nlnRSS/n

其中，第一项表示模型拟合的情况，第二项表示对模型复杂度的惩罚，最终达到满足模型有效性和可靠性条件下参数个数最少的目的.

2.3　QFM⁃CGC算法描述

根据上述推导和分析过程，总结QFM⁃CGC算法如下所示.

算法

基于分位数因子模型的条件Granger因果分析

输入：时间序列X，X=X1,X2,⋯,XKT∈Rk×n

输出：因果关系连接矩阵

（1）ADF检验：

ΔXt=α+βt+δXt-1+∑i=1mβiΔXt-i+εt

（2）for i=1:K

（3） for j=1:K

（4） if i=j，结束本次循环

（5） else

（6） 根据计算得到的最小AIC确定模型阶数

（7） for p=1:pmax（pmax是时间序列X的最大滞后阶数）

（8） Z=setdiffXi,j,Xi,1,Xi,2,⋯,Xi,pmax,Xj,p

（9） 对条件变量进行分位数因子分析得到降维后的ZNEW

（10） 对Xi，Xj,p和条件变量ZNEW进行条件Granger因果分析建模

（11） if pvalue>0.9，即通过显著性检验

存在Xi→Xj,p的因果关系

（12） else 不存在Xi→Xj,p的因果关系

（13） end

（14） end

（15） end

（16）end

3 仿真实验与分析

在仿真模拟研究中，比较QFM⁃CGC，经典方法CGC［8］和基准方法PCA⁃CGC［21］，mBTS⁃CGC［22］，PMIME［23］的性能.Geweke［8］向VAR模型中引入条件变量，提出条件Granger因果模型，改善了传统方法无法判断直接因果关系和间接因果关系的缺陷.Zhou et al［21］提出PCA⁃CGC方法，将PCA与条件Granger因果模型相结合来处理高维大脑神经网络的计算，与传统方法相比，降低了计算成本.Siggiridou and Kugiumtzis［22］采用backward⁃in⁃time方法对每个变量的滞后阶数使用有监督的逐步向前选择，有效减少VAR模型阶数，并与条件Granger因果模型结合，提出mBTS⁃CGC方法.Kugiumtzis［23］将度量混合嵌入的条件互信息（Conditional Mutual Information from Mixed Embedding，MIME）拓展到多变量时间序列，形成可以检测直接耦合的部分MIME （Partial MIME，PMIME）.PMIME在由非均匀嵌入方案导出的滞后变量X，Y和Z的联合状态空间的子空间中重构一个点（向量），目的是最好地解释Y的演化，得到的混合嵌入向量只包含所有变量中最相关的成分，避免大维度会恶化估计的情况.

实验考虑的仿真模拟系统是两个不同维度变量的线性仿真系统，且多项式平稳随机.生成的多变量时间序列的平稳性要求每个时间序列的数据在时间函数的合理范围内进行经验评估，实验结果在显著性水平α=0.1下确定［24］.

将QFM⁃CGC方法运用到宏观经济时间序列并建立预测模型对因果分析结果进行验证，最后将仿真结果与CGC，PCA⁃CGC，mBTS⁃CGC和PMIME进行对比.

3.1　多变量线性时间序列

在两个不同维度的线性系统的仿真模拟时间序列上评估了因果关系测度的集合，共计19个，且两个随机系统均为假设.仿真系统如下所示.

第一组数据是一个5维变量的4阶线性VAR系统VAR5（5）［25］.由式（10）产生：

X1,t=0.3X1,t-1+0.4X2,t-3+ϵ1,tX2,t=0.4X2,t-1+0.4X5,t-3+ϵ2,tX3,t=0.4X3,t-2-0.4X1,t-1+ϵ3,tX4,t=0.4X4,t-1+0.2X4,t-3-0.4X2,t-1+ϵ4,tX5,t=0.4X3,t-1+0.4X5,t-2-0.4X4,t-1+ϵ5,t

其中，ϵi i=1,…,5表示高斯白噪声序列.时间序列长度N=500.仿真系统VAR5（5）中真实存在的因果关系为X1→X3，X2→X1，X2→X4，X3→X5，X4→X5和X5→X2，共计六个.其因果关系如图1所示.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   VAR5(5)真实因果关系(黑色表示存在因果关系)

Fig.1   The real causality diagram of VAR5(5) (Black indicates causality)


利用AIC算法选择最佳模型阶数，图2a~e分别代表目标变量为X1,X2,X3,X4,X5的AIC算法的实验结果，最小的AIC对应最优的延迟阶数.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   VAR5(5)模型阶数选择

Fig.2   The order selection of VAR5(5)


表1为VAR5（5）在100次蒙特卡洛实验中因果关系的p>0.9的频率，其中，选择QFM⁃CGC方法的分位数为0.5.由表可见，CGC和PCA⁃CGC不能完全正确识别式（10）中的因果关系.其中，CGC识别正确因果关系X3→X5的频率仅为3%，错误识别直接因果关系X1→X2的概率高达98%，而PCA⁃CGC除了X1→X3，其他识别正确因果关系的概率均未超过50%.虽然PMIME，mBTS⁃CGC和QFM⁃CGC都能正确识别所有因果关系，但仅有CGC和PMIME受到虚假的因果关系影响，尤其是PMIME受到比CGC更多的虚假的因果关系X1→X2，X2→X3，X3→X1，X3→X2，X4→X1，X4→X2，X4→X3，X5→X1和X5→X3的干扰，假阳性更高.和mBTS⁃CGC和QFM⁃CGC相比，综合来看，本文方法QFM⁃CGC识别正确因果关系的概率更高.

表1   VAR5(5)的100次蒙特卡洛实验中因果关系结果的频率

Table 1  The frequency of causality results in 100 Monte Carlo implementations of VAR5(5)

方法CGCPMIMEPCA⁃CGCmBTS⁃CGCQFM⁃CGC
X1→X298%1%000
X1→X399%95%100%94%99%
X2→X198%99%40%72%99%
X2→X304%000
X2→X4100%98%50%98%100%
X3→X108%000
X3→X201%000
X3→X53%86%43%96%100%
X4→X104%000
X4→X201%000
X4→X302%000
X4→X5100%97%26%100%99%
X5→X102%000
X5→X2100%99%18%100%99%
X5→X301%000

新窗口打开| 下载CSV


第二组数据是一个10维变量的4阶线性VAR系统VAR10（5）.由式（11）产生：

X1,t=0.4X1,t-3-0.4X1,t-4+0.4X5,t-3+ϵ1,tX2,t=0.4X2,t-1-0.3X2,t-5+0.5X1,t-1+ϵ2,tX3,t=0.4X3,t-5+0.3X3,t-4-0.4X5,t-3+ϵ3,tX4,t=0.6X4,t-3+0.3X5,t-3-0.4X2,t-4+ϵ4,tX5,t=0.3X5,t-1+0.4X5,t-3-0.3X4,t-1+ϵ5,tX6,t=0.4X1,t-3+0.4X6,t-2-0.6X8,t-2+ϵ6,tX7,t=0.3X7,t-1-0.4X6,t-4+0.3X7,t-3+ϵ7,tX8,t=0.5X8,t-4-0.2X8,t-2+0.4X5,t-3+ϵ8,tX9,t=0.5X9,t-3-0.5X10,t-3+0.4X1,t-3+ϵ9,tX10,t=0.5X7,t-2-0.5X10,t-3-0.3X10,t-1+ϵ10,t

其中，ϵi i=1,…,10表示高斯白噪声序列.仿真系统VAR10（5）中真实存在的因果关系为X1→X2，X1→X6，X1→X9，X2→X4，X4→X5，X5→X1，X5→X3，X5→X4，X5→X8，X6→X7，X7→X10，X8→X6和X10→X9，共计13个.其因果关系如图3所示.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   VAR10(5)真实因果关系(黑色表示存在因果关系)

Fig.3   The real causality diagram of VAR10(5) (Black indicates causality)


图4是利用AIC算法选择最佳模型阶数，图4a~j分别代表目标变量为X1，X2，X3，X4，X5，X6，X7，X8，X9，X10的AIC算法的实验结果，最小的AIC对应最优的延迟阶数.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   VAR10(5)模型阶数选择

Fig.4   The order selection of VAR10(5)


表2为VAR10（5）的100次蒙特卡洛实验中因果关系的p>0.9的频率，选择QFM⁃CGC的分位数为0.5.与线性VAR系统VAR5（5）的实验结果相似，CGC和PCA⁃CGC无法完全正确识别式（11）中的因果关系，其中CGC识别正确因果关系X2→X4，X6→X7的频率为0，PCA⁃CGC识别正确因果关系X1→X9的概率未超过50%.本仿真系统中PMIME仍受许多虚假的因果关系干扰，识别正确因果关系X5→X4的概率仅为69%.mBTS⁃CGC和QFM⁃CGC都能正确识别所有因果关系，但mBTS⁃CGC正确识别的概率不高.虽然本文方法错误识别了虚假因果关系X6→X8，概率为8%，但综合比较，本文方法识别正确因果关系的概率更高，表现更好.

表2   VAR10(5)的100次蒙特卡洛实验中因果关系结果的频率

Table 2  The frequency of causality results in 100 Monte Carlo implementations of VAR10(5)

方法CGCPMIMEPCA⁃CGCmBTS⁃CGCQFM⁃CGC
X1→X299%99%100%99%100%
X1→X632%95%75%57%91%
X1→X935%98%40%55%100%
X2→X40100%85%99%100%
X2→X501%000
X4→X588%99%099%100%
X5→X151%85%100%57%100%
X5→X336%100%100%99%100%
X5→X495%69%100%97%100%
X5→X823%99%100%80%94%
X6→X70100%99%89%100%
X6→X800008%
X7→X10100%100%100%99%100%
X8→X6100%96%100%99%100%
X9→X801%000
X10→X801%000
X10→X9100%100%100%99%100%

新窗口打开| 下载CSV


虽然CGC理论上能正确区分直接因果关系和间接因果关系，但实验结果表明，CGC在仿真系统VAR5（5）中仍然错误识别了X1→X2因果关系，这可能是受到间接因果关系X1→X3→X5→X2的影响.由于PMIME是基于KNN算法的，该算法受维度灾难的影响，对于多维度的数据处理不准确，因为随着维度的增加，“看似相近”的两个点的距离越来越大，就会越来越“不像”，对于高度依赖距离的KNN算法其结果会影响准确率.PCA⁃CGC和QFM⁃CGC虽然方法类似，但PCA在降维时没有捕获隐藏的因素，尤其是这些因素可能改变时间序列的分布特征，造成PCA⁃CGC不能准确识别因果关系.mBTS⁃CGC对每个变量的滞后阶数使用有监督的逐步向前选择，有效地减少VAR模型阶数，但在噪声的干扰下，其中一个条件变量选择错误会引起其他因果关系的判断不准确，最终造成该方法的假阴性较高.

3.2　宏观经济时间序列

使用宏观经济时间序列进行因果分析并建立预测模型，对因果分析结果进行验证，主要目标是从高维宏观经济变量面板中找寻与实际GDP变化趋势有因果关系的时间序列.该数据集由1960年第一季度至2019年第二季度的211个美国宏观经济变量组成（N=211，T=238），其中的数据会及时更新，可以在网站（http：∥research.stlouisfed.org/econ/mccracken/）免费下载.计算之前，对每个序列进行平稳性处理，代码也可以在FRED⁃QD数据网站上获得.利用因果关系的方法找出影响宏观经济变量的主要因素，剔除无关变量，保留相关变量，并将该相关变量作为预测模型的输入进行建模预测，根据预测误差反向验证因果分析方法的有效性.

与Chen et al［13］相同，设置估计量的最大因子数k=8.使用秩最小化估计器［13］估计分位数为0.01,0.05,0.1,0.25,0.75,0.9,0.95,0.99时的因子估计数如表3所示.由表可见，QFA因子的数量在不同分位数之间存在显著差异，表明该数据集存在非标准因子结构.为了比较QFA因子和PCA因子，将QFA因子的每个元素与选择的八个PCA因子进行回归并计算这些回归中的R2，结果如表4［13］所示.很明显，当τ接近0.5时，QFA因子与PCA因子高相关，R2均在0.9以上.相比之下，τ=0.9时的第一个QFA因子（表示为F^QFA0.90）和τ=0.95，0.99时的第一个QFA因子（分别表示为F^QFA0.95和F^QFA0.99）与PCA因子的相关性较低，R2低于0.4.因此，F^QFA0.90，F^QFA0.95和F^QFA0.99包含可能有助于预测宏观经济变量的额外信息，在此应用程序中有使用QFA的空间.由表4可得，由于F^QFA0.90，F^QFA0.95和F^QFA0.99的R2分别为0.316，0.261和0.266，与其他QFA因子相比，F^QFA0.95与F^QFA0.90和F^QFA0.99有非常高的相关性，它们具有类似的捕获额外信息的能力，这些信息能够帮助预测宏观经济变量.因此，在后续分析中重点关注F^QFA0.90和F^QFA0.99的预测能力.

表3   不同分位数下的因子估计数

Table 3  Estimation of factors at different quantiles

分位数τ因子个数
0.011
0.051
0.102
0.254
0.505
0.755
0.902
0.951
0.991

新窗口打开| 下载CSV


表4   F^QFA和F^PCA的比较结果

Table 4  Comparison of F^QFA and F^PCA

分位数τF̂QFAτ的元素的个数
12345
0.010.657
0.050.733
0.100.7960.871
0.250.9520.9320.9390.890
0.500.9930.9760.9640.9450.923
0.750.9060.9450.9430.9030.882
0.900.3160.911
0.950.261
0.990.266

新窗口打开| 下载CSV


使用不同方法进行因果关系分析后，选出与目标变量具有因果关系的原因变量作为模型的输入进行预测，并对分析结果进行进一步的验证.采用CNN⁃LSTM预测模型来分析每一种方法得出的因变量进行建模的预测效果，进行30次实验，取平均值来消除偶然因素对实验结果的影响.最后，采用均方根误差（RMSE）、平均绝对百分误差（MAPE）和对称平均绝对百分比误差（SMAPE）三个指标来定量评价预测精度，三个评价指标的定义如下：

RMSE=1n∑i=1nyi-y^i2
MAPE=100%n∑i=1ny^i-yiyi
SMAPE=100%n∑i=1ny^i-yiyi+y^i2

其中，yi和y^i分别是真实值和预测值，n是样本个数.

表5是不同方法预测GDP的精度比较，表中黑体字表示最优值.图5~10展示了不同方法预测的GDP变化趋势.由表5可见，本文方法在0.90分位数条件下的RMSE，MAPE和SMAPE都是最小的，并在预测对比图中具有明显的优越性拟合效果，能更精准地追踪GDP的变化趋势.

表5   GDP预测结果

Table 5  The prediction of GDP

对比方法因变量（编号）RMSEMAPESMAPE
CGC

28，64，74，104，

116，162

2.535582.303541.06974
mBTS⁃CGC4，9，11，16，18，26，36，60，66，86，89，102，138，141，148，2032.378872.017091.06333
PCA⁃CGC6，70，71，77，141，1482.246802.065301.03592
PMIME70，137，161，1631.909921.544850.95669
QFM⁃CGC (τ=0.90)

2，4，7，10，21，

79，163，171

1.703791.468410.93114
QFM⁃CGC (τ=0.99)

2，4，7，10，21，26，

79，160，163，189

1.808521.587680.95876

新窗口打开| 下载CSV


图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   PCA⁃CGC的GDP预测图

Fig.5   The prediction of GDP by PCA⁃CGC


图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   mBTS⁃CGC的GDP预测图

Fig.6   The prediction of GDP by mBTS⁃CGC


图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   CGC的GDP预测图

Fig.7   The prediction of GDP by CGC


图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   PMIME的GDP预测图

Fig.8   The prediction of GDP by PMIME


图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   QFM⁃CGC (τ=0.90)的GDP预测图

Fig.9   The prediction of GDP by QFM⁃CGC (τ=0.90)


图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   QFM⁃CGC (τ=0.99)的GDP预测图

Fig.10   The prediction of GDP by QFM⁃CGC (τ=0.99)


从表5还可以看出，QFM⁃CGC识别出对GDP具有因果关系的变量主要与个人消费支出、私人固定投资、生产制造和消费有关.内需、投资和出口俗称拉动经济增长的“三驾马车”，尤其是消费需求是生产的目的，消费可以创造出生产的动力，并刺激投资需求，以此促进经济发展.然而，CGC未能识别出投资与GDP具有因果关系，因此使用CGC来预测GDP造成的误差最大.虽然mBTS⁃CGC可以识别出许多因变量，但其中可能包含错误的因变量（如货币存量、国库券等），这些因变量会干扰预测结果，导致预测误差较大.PCA⁃CGC和PMIME识别出制造业和非监督员工的平均每周工作小时数与GDP之间有Granger因果关系，加班和额外的工作时间可能会增加生产和服务活动，对GDP产生积极影响.然而，过度的长时间工作可能导致劳动力疲劳，影响效率或产生健康问题，可能减少GDP.尽管制造业和非监督员工的平均每周工作小时数可能与GDP相关，但它们之间不一定存在Granger因果关系，因为过高或过低都会对GDP产生负面效应.实验结果与GDP影响关系是一致的，进一步验证了本文方法的有效性.

3.3　北京AQI及气象时间序列

使用北京AQI及气象数据集进行因果分析，并建立预测模型，对因果分析的结果进行验证.该数据集选用2016年1月1日到2016年6月15日的每小时数据，共4008个样本，每个样本包括6维AQI时间序列和5维气象时间序列，详细描述见表6.

表6   北京AQI及气象时间序列编号及变量对照表

Table 6  The number and variable comparison table of Beijing AQI and meteorological time series

编号123456
变量PM2.5PM10SO2NO2COO3
编号7891011
变量气温气压露点降雨量风速

新窗口打开| 下载CSV


对每个对比模型的因果分析结果进行预测，根据预测的准确性来判断不同因果分析方法的有效性.与上节相同，设置估计量的最大因子数k=8.

使用Chen et al［13］的秩最小化估计器在分位数为0.01,0.05,0.1,0.25,0.75,0.9,0.95,0.99时进行因子估计，不同分位数条件下的因子估计数同为1.此外，将QFA因子的每个元素与选择的八个PCA因子进行回归并计算这些回归中的R2，分位数在0.01,0.05,0.1,0.25,0.75,0.9,0.95,0.99时的R2分别为0.794，0.786，0.999，0.999，0.999，0.984，0.945，0.961.τ=0.01时的QFA因子（表示为F^QFA0.01）和τ=0.05时的QFA因子（表示为F^QFA0.05）与PCA因子的相关性较低，R2低于0.8.因此，F^QFA0.01和F^QFA0.05包含可能有助于预测宏观经济变量的额外信息，在此应用程序中有使用QFA的空间.

使用不同的方法进行因果关系分析后，选出与目标变量具有因果关系的原因变量作为模型的输入进行预测，并对分析结果进行进一步的验证.表7是不同方法对NO2的因果关系分析结果并使用得到的因变量采用CNN⁃LSTM预测模型预测NO2的精度比较，表中黑体字表示最优值.

表7   NO2的预测结果

Table 7  The prediction of NO2

对比方法因变量（编号）RMSEMAPESMAPE
mBTS⁃CGC5，9，10，115.895890.496750.02104
PCA⁃CGC2，5，93.370850.567040.02329
PMIME，CGC6，7，115.372640.498210.02043

QFM⁃CGC

τ=0.01,0.05

1，2，3，5，72.464990.424110.01817

新窗口打开| 下载CSV


从表7可以发现，本文方法在低分位数0.01和0.05时的RMSE，MAPE和SMAPE都是最小的.此外，结合气象学等领域的实际背景，进一步分析因果关系分析结果的应用价值.露点是空气中水蒸气达到饱和所需的温度，而SO2是一种气体，在大气中可能对人类健康和环境造成负面影响，虽然SO2的排放会随着湿度升高而下降，但这不是因为露点与SO2之间存在因果关系，mBTS⁃CGC和PCA⁃CGC错误判断了露点与SO2之间的因果关系.PM2.5，PM10和SO2都是大气污染物，SO2在大气中与水蒸气、氧气等物质相互作用，形成硫酸盐颗粒物，并与大气中的其他颗粒物结合生成复合颗粒物.这些复合颗粒物包括PM2.5和PM10等，因此SO2与PM2.5和PM10之间可能存在一定的因果关系，然而其他方法没有捕获到因变量PM2.5和PM10.因此，通过从理论出发分析所提模型和对比模型的因果关系分析的结果，进一步证明了本文方法的有效性.

4 结论

传统的Granger因果模型由VAR模型推导，但当VAR模型中涉及大量待估计系数时容易产生维度灾难，因此，使用一般的因果分析方法在高维时间序列中难以准确地发现因果关系，目前对于高维和超高维的时间序列因果分析缺少有效的处理方法.针对高维线性时间序列因果关系的识别问题，本文提出QFM⁃CGC方法，使用分位数因子模型对条件Granger因果分析方法的条件变量进行降维，解决了VAR模型中存在大量待估计系数的问题，可以有效识别高维线性时间序列的因果关系，尤其是在数据包含异常值时，使用分位数因子模型降维具有更明显的优势.大量实验证明，将降维技术与条件Granger因果分析相结合，能准确识别直接因果关系和间接因果关系.

本文提出的方法能在平稳和线性的时间序列下进行建模，并实现较好的因果分析结果，但现实中部分数据具有非线性特征，因此，未来将对高维非线性时间序列因果关系的分析展开研究.此外，当时间序列非平稳时，处理时间序列平稳化的过程可能会使原时间序列的结果和意义发生变化，导致因果分析方法的解释意义发生变化.

参考文献

Granger C W J.

Investigating causal relations by econometric models and cross⁃spectral methods

Econometrica，1969，37(3)：424-438.

[本文引用: 1]

Wismüller A， Vosoughi M A， DSouza A，et al.

Exploring directed network connectivity in complex systems using large⁃scale augmented Granger causality

∥Proceedings of SPIE 12033，Medical Imaging 2022：Computer⁃Aided Diagnosis. San Diego，CA，USA：SPIE，2022：168-177.

[本文引用: 1]

Maradana R P， Pradhan R P， Dash S，et al.

Innovation and economic growth in European Economic Area countries：The Granger causality approach

IIMB Management Review，2019，31(3)：268-282.

[本文引用: 1]

Billio M， Getmansky M， Lo A W，et al.

Econometric measures of connectedness and systemic risk in the finance and insurance sectors

Journal of Financial Economics，2012，104(3)：535-559.

[本文引用: 1]

Chang T， Gupta R， Inglesi⁃Lotz R，et al.

Renewable energy and growth：Evidence from heterogeneous panel of G7 countries using Granger causality

Renewable and Sustainable Energy Reviews，2015(52)：1405-1412.

[本文引用: 1]

Hlinka J， Hartman D， Vejmelka M，et al.

Reliability of inference of directed climate networks using conditional mutual information

Entropy，2013，15(6)：2023-2045.

[本文引用: 1]

Blinowska K J， Kuś R， Kamiński M.

Granger causality and information flow in multivariate processes

Physical Review E，2004，70(5)：050902.

[本文引用: 1]

Geweke J.

Measurement of linear dependence and feedback between multiple time series

Journal of the American Statistical Association，1982，77(378)：304-313.

[本文引用: 4]

李松，胡晏铭，郝晓红,等.

基于维度分组降维的高维数据近似k近邻查询

计算机研究与发展，2021，58(3)：609-623.

[本文引用: 1]

Li S， Hu Y M， Hao X H，et al.

Approximate k⁃nearest neighbor query of high dimensional data based on dimension grouping and reducing

Journal of Computer Research and Development，2021，58(3)：609-623.

[本文引用: 1]

刘淑伟，陈威，赵伟，等.

基于簇内乘积量化的最近邻检索方法

计算机学报，2020，43(2)：303-314.

[本文引用: 1]

Liu S W， Chen W， Zhao W，et al.

Nearest neighbor search based on product quantization in clusters

Chinese Journal of Computers，2020，43(2)：303-314.

[本文引用: 1]

Wold S， Esbensen K， Geladi P.

Principal component analysis

Chemometrics and Intelligent Laboratory Systems，1987，2(1-3)：37-52.

[本文引用: 1]

Abdi H， Williams L J.

Principal component analysis

WIREs Computational Statistics，2010，2(4)：433-459.

[本文引用: 1]

Chen L， Dolado J J， Gonzalo J.

Quantile factor models

Econometrica，2021，89(2)：875-910.

[本文引用: 10]

Mooney C Z. Monte Carlo simulation. Thousand Oaks：Sage Publications，1997，103.

[本文引用: 1]

Pankratz A. Forecasting with dynamic regression models. Hoboken：John Wiley & Sons，2012，400.

[本文引用: 1]

Brandt P T， Williams J T.

Multiple time series models

Sage Publications，2006，120.

[本文引用: 1]

Guo S X， Seth A K， Kendrick K M，et al.

Partial Granger causality⁃eliminating exogenous inputs and latent variables

Journal of Neuroscience Methods，2008，172(1)：79-93.

[本文引用: 1]

Dickey D A， Fuller W A.

Distribution of the estimators for autoregressive time series with a unit root

Journal of the American Statistical Association，1979，74(366)：427-431.

[本文引用: 1]

Barber R F， Drton M.

High⁃dimensional Ising model selection with Bayesian information criteria

Electronic Journal of Statistics，2015，9(1)：567-607.

[本文引用: 1]

Akaike H.

A new look at the statistical model identification

IEEE Transactions on Automatic Control，1974，19(6)：716-723.

[本文引用: 1]

Zhou Z Y， Chen Y H， Ding M Z，et al.

Analyzing brain networks with PCA and conditional Granger causality

Human Brain Mapping，2009，30(7)：2197-2206.

[本文引用: 2]

Siggiridou E， Kugiumtzis D.

Granger causality in multivariate time series using a time⁃ordered restricted vector autoregressive model

IEEE Transactions on Signal Processing，2016，64(7)：1759-1773.

[本文引用: 2]

Kugiumtzis D.

Direct⁃coupling information measure from nonuniform embedding

Physical Review E，2013，87(6)：062918.

[本文引用: 2]

Quiroga R Q， Kraskov A， Kreuz T，et al.

Performance of different synchronization measures in real data：A case study on electroencephalographic signals

Physical Review E，2002，65(4)：041903.

[本文引用: 1]

Jia Z Y， Lin Y F， Liu Y X，et al.

Refined nonuniform embedding for coupling detection in multivariate time series

Physical Review E，2020，101(6)：062113.

[本文引用: 1]

/

〈 〉

版权所有 © 2021 《南京大学学报(自然科学)》

地址：江苏省南京市鼓楼区汉口路22号，南京大学学报编辑部，210093

电话：025-83592704 , 传真：025-83592704        技术支持：北京玛格泰克科技发展有限公司