稀疏组Lasso Granger因果图模型及应用

高 伟,杨海忠,杨 露

(西安财经大学 统计学院,陕西 西安 710100)

一、引 言

时间序列分析的主要任务之一是发现序列间的因果联系。Granger因果性概念度量一个时间序列是否有助于解释另一个时间序列未来取值的变化,自提出以来,在多维时间序列分析中发挥了重要作用,并广泛应用于经济学、生物学和机器学习等领域检验序列间的因果关系[1-2]。用顶点表示多维时间序列的分量序列,顶点之间的有向边表示序列间的Granger因果关系,就将整个时间序列系统映射为一个Granger因果图。借助图模型结果的可视化和可解释的优点,Granger因果图广泛应用于建模和分析多维时间序列分量间复杂的因果关系[3-6]。

Granger因果图结构学习的一个关键问题是检验序列间的因果关系。传统的Granger因果关系检验基于VAR模型的方法,钟意等通过金融压力指数和宏观经济变量的Granger因果关系检验研究金融不稳定对主要宏观经济变量的溢出效应[2]。白仲林等运用SVAR-IV方法识别预估计产业指导目录政策冲击对宏观经济系统变量的动态因果效应[5]。更多详细介绍可以参考任伟杰和韩敏的综述文章[6]。

然而,随着数据采集和存储技术的发展,时间序列数据的维度和规模迅速增加,传统的基于VAR模型参数估计的方法涉及过多未知参数,给估计和推断带来困难。在很多情况下,一个序列只和系统中少数几个序列存在直接联系,这时VAR模型就存在过参数化问题。Lasso类方法通过对模型系数的稀疏性约束,在参数估计的同时实现变量选择,在高维VAR模型的参数估计和预测中得到广泛研究[7-9]。基于VAR模型的稀疏化约束,将Granger因果关系检验问题转化为变量选择问题,Arnold等提出基于Lasso的因果性概念,称为图Lasso Granger方法(Graphical Lasso Granger,简记为GLG),应用Lasso估计的回归系数辨识Ganger因果关系的存在性,并在多个领域得到应用[10-11]。Lasso方法对VAR模型系数矩阵的所有参数施加同等约束,因此存在变量过选择问题。组Lasso方法考虑到变量的分组结构,提高了变量选择的效率。将来自同一序列的不同滞后变量分为一组,Lozano等提出图组Lasso Granger方法研究多元Granger因果关系,利用变量的分组结构大大缩减了计算量[12]。基于组Lasso的Granger因果检验方法的统计性质和应用也得到进一步研究[13-14]。

组Lasso方法考虑到滞后变量按照其所在的序列进行分组,使得因果关系检验更有效。但现有的研究存在两个方面的问题:其一,组Lasso Granger方法将同一序列的不同滞后变量分为一组,约束其系数同时都为0或都不为0,提高检验效率的同时忽视了序列影响的时间滞后信息[12-14]。时间序列间的Granger因果关系包含的另一个重要信息是确定因果影响的滞后阶数,传统的Granger检验方法可以检验因果影响的具体滞后阶数[1-2]。Lasso Granger方法估计系数矩阵的所有参数,包含了因果影响的滞后信息,但存在计算复杂和过选择问题。其二,对于Lasso类方法在Granger图模型结构学习的模拟研究主要是少数几个时间序列,并且结构固定,对于该方法在各种维数和相依结构下的有效性和大样本性质没有研究[10-14]。

针对以上问题,本文主要研究内容如下:第一,提出稀疏组Lasso Granger因果图模型,应用稀疏组Lasso方法检验多维时间序列Granger因果关系的存在性和因果影响的滞后阶数。Simon等提出稀疏组Lasso算法,是组lasso方法与lasso方法的结合,既能实现组内稀疏又能实现组间稀疏[15]。通过约束同一序列的滞后变量同时都为0或都不为0检验Granger因果关系的存在性,同时对组内不同滞后变量的约束可以实现组内变量选择,进一步确定Granger影响的具体滞后变量。郭金良等用稀疏组lasso算法对获取的Granger因果特征值进行特征筛选,但并没有进行滞后期检验[16]。第二,在不同样本量和不同维数下对各种因果相依结构的时间序列数据进行模拟研究,验证了提出的方法对Granger因果关系检验和滞后变量确定的有效性。第三,国内时间序列分析实证研究中,分析变量之间的Granger因果关系一直是研究的重点内容。目前,Granger因果关系和图模型的研究主要应用于金融市场[2,4-5]。由于宏观经济变量数据多为年度数据,样本量较少,传统方法对于多维宏观经济变量的因果分析缺少有效的处理手段,主要集中于2个或3个变量的分析[17-18]。但宏观经济系统的变量之间存在直接或间接的联系,遗漏相关变量容易产生虚假因果现象。本文将稀疏组Lasso方法应用到8个宏观经济变量,建立稀疏组Lasso Granger因果图,分析变量间的直接因果关系和因果影响的具体滞后阶数,对于宏观经济形势和政策制定的研究具有十分重要的意义。

二、稀疏组Lasso Granger因果图模型

本节给出稀疏组Lasso Granger因果图模型的定义和估计过程。

(一)模型设定

Granger因果性概念推广到d维(d>2)时间序列X(t)=(X1(t),X2(t),…,Xd(t))T,定义如下:对于序列Xi,在给定所有其他序列Xm(m=1,2,…,d,m≠j)的过去值条件下,如果加入Xj的过去值可以改进Xi未来值预测的精度,就称Xj是Xi的Granger原因。

设d维时间序列X(t)=(X1(t),X2(t),…,Xd(t))T,t∈Z来自于VAR(p)模型,

(1)

其中A(1),A(2),…,A(p)为d×d系数矩阵,

ε(t)=(ε1(t),ε2(t),…,εd(t))T是d维误差向量,与X(t)不相关。假设X(t),ε(t)都是0均值,ε(t)的协方差矩阵Cov(ε(t))是对角阵。

多元Granger因果关系通过检验式(1)中VAR模型的系数矩阵A(k)得到。

(二)稀疏组Lasso Granger因果关系检验方法

(2)

Bolstad等指出,多数情况下序列只和系统中部分序列存在直接的联系,提出了稀疏VAR模型的概念[13]。本文假设模型满足如下的稀疏条件:

假设1 定义因果关系集合S,如果序列Xj(t)是序列Xi(t)(1≤i,j≤d)的Granger原因,那么(j→i)∈S。稀疏性条件要求集合S中的元素个数|S|≤Md,M>1。

设Z(t)=(X1(t-1),…,X1(t-p),…,Xd(t-1),…,Xd(t-p))为1×pd维滞后变量集合。令Zj(t)=(Xj(t-1),…,Xj(t-p))∈R1×p表示Z(t)中序列Xj(t),j=1,2,…,d的滞后变量。按照滞后变量所属的序列,将Z(t)中的滞后变量分成d个组,Z1(t),…,Zd(t)。按照分组,式(1)的VAR(p)模型可改写为:

(3)

矩阵表示:

X(t)=BZ(t)T+ε(t)

(4)

Xi(t)=Bi·Z(t)T+εi(t)

(5)

这样系数矩阵B的估计问题分解成d个优化问题的集合,每一个子问题处理一个顶点。

给定模型中自变量集合分组的d组变量,式(5)的稀疏组Lasso估计求解如下优化问题,

(6)

稀疏组Lasso估计中,参数α控制组内稀疏性,1-α越大,组内不为0的参数越少。参数α=1时,稀疏组Lasso估计转化为组Lasso估计,组Lasso方法将每一个序列的不同期滞后变量分为一组,通过组间约束使得对同组变量系数的估计同时都为0或同时都不为0,提高了Granger检验的效率。但组Lasso约束损失了变量影响的时间滞后信息。本文提出稀疏组Lasso方法检验Granger因果联系,将组lasso方法与lasso方法相结合,既能实现组内稀疏又能实现组间稀疏,组间稀疏约束Granger因果联系的存在性,组内稀疏进一步选择导致因果影响的具体滞后变量。式(6)可通过次梯度算法求解[19]。

(7)

为方便,计算时令λ1=λnα,λ2=λn(1-α)。本文用BIC准则选择罚参数,即选择λ1,λ2,τ使式(8)达到最小,

(8)

参数λ1,λ2,τ的选择采用格点搜索法,例如:可首先设定搜索区间为(0,1)×(0,1)×(0,c),常数c>0,以0.1为间隔,计算各种组合下的BIC值;如果最优搜索结果在区间内,则停止;如果最优搜索结果在边界上,如0.1,0.9或(c-0.1),则重新选择搜索区间(0,0.1),(0.9,1)或(c-0.1,c)以0.01为间隔进行逐步搜索,以此进行下去,直到得到最优的λ1,λ2,τ值。

(三)图模型的建立

(9)

根据定义2,稀疏组Granger因果图中顶点集表示多维时间序列的分量序列,顶点间的有向边表示对应序列间的直接Granger因果关系,借助图的直观表达可以同时反映序列之间的复杂联系,有效区别直接因果关系和由中间因素影响导致的间接因果关系。滞后信息矩阵进一步给出了因果影响的具体滞后阶数。

另外,本文的图模型仅考虑序列Xj和Xi(j≠i)间的Granger因果关系,不考虑自身的因果关系,即Xi(t-k),k=1,2,…,p对Xi的影响,其在图中是一个自环。本文构造的滞后变量集Z(t)仍包含滞后变量Xi(t-k),k=1,2,…,p,主要为了研究在给定所有滞后变量的条件下,序列Xj对Xi的直接因果关系。

由数据建立稀疏组Lasso Granger因果图模型的算法步骤如下:

步骤1:输入样本数据Xi(t),i=1,2,…,d,t=1,2,…,n,滞后阶数p,参数λ1,λ2,τ,构造样本矩阵Z(t)=(Z1(t),Z2(t),…,Zd(t)),Zj(t)=(Xj(t-1),Xj(t-2),…,Xj(t-p)),j=1,2,…,d;

三、模拟比较分析

本节用数值模拟方法验证稀疏组Lasso方法的有效性,并与组Lasso Granger方法比较,验证其检验Granger因果影响具体滞后变量的优越性,最后验证了稀疏组Lasso方法的稳健性。本节模拟中的两个模型,模型1中|S|=5,模型2中|S|=d-1,均满足稀疏性条件假设1。

对于Granger因果图模型结构学习,用F1评分作为评价结果的优劣的度量[12]。F1评分是精确率(Precision)和召回率(Recall)的调和均值,定义如下:

(10)

其中,P表示精确率,是正确辨识的边数和辨识的总边数的比值,R表示召回率,是正确辨识的边数占模型真实的边数的比重。F1评分在精确率和召回率之间平衡自然选择的评估指标,只有精确率和召回率都大时,才能得到比较大的F1值。

(一)模型结构和参数固定情形的有效性验证

本小节利用结构和参数固定,Granger因果关系有不同的滞后期数的模型,考察不同样本量下,稀疏组Lasso和组Lasso方法在存在因果影响的序列不同滞后期上的非0比率,验证稀疏组Lasso方法辨识因果影响具体滞后信息的能力。

模型1:设5维时间序列数据产生如下VAR(3)模型:

(11)

其中εi(t)～N(0,0.12),i=1,2,3,4,5是独立的高斯噪声。真实的因果图结构见图1。

图1 模型1的Granger因果图

分别产生n=100,200,500个样本观测值,用组Lasso和稀疏组Lasso方法检验Granger因果关系,建立Granger因果图。用第2节的BIC准则选择罚参数λ,λ1,λ2,τ,两种方法的门限参数τ=0.05。另外,在非0的Granger因果关系Xj→Xi上,进一步判断影响的滞后信息,即计算1 000次模拟中,滞后变量Xj(t-1),Xj(t-2),Xj(t-3)的系数非0的比率,作为评价两种方法在检验因果影响具体滞后值的性能度量。结果见表1。

表1给出了两种方法1 000次模拟的平均F1评分,值都大于0.8,且都随样本量增加而显著提高,表明两种方法都以高的正确率检验出Granger因果关系。相同样本量下,稀疏组Lasso方法要优于组Lasso方法,样本量较小时优势更明显,表明了稀疏组Lasso方法在Granger因果图模型结构辨识上的优良性。

接下来比较两种方法对于Granger因果关系滞后期的检验效果。由表1可知,两种方法都能以接近或等于1的比率检测出正确的滞后期,但组Lasso方法在其余两个滞后期上不为0的比例也相当高,且随着样本量的增加,在无影响的滞后期上非0比率降低不明显,只有在样本量n=500时,X1(t-3)→X4(t),X5(t-1)→X4(t),X4(t-2)→X5(t),X4(t-3)→X5(t),X4(t-3)→X4(t),X5(t-3)→X5(t)的比率小于0.1。而稀疏组Lasso方法在其余两个滞后期上不为0的比率小得多,并且随着样本量的增加,在其余两个滞后期上为0的比率显著降低。在样本量n=100时,基本在0.1左右;样本量n=200时,除X1(t-1)→X4(t)上比率为0.081外,其他值都小于0.05;样本量n=500时,无影响的滞后值非0比率除X1(t-1)→X4(t)和X4(t-2)→X4(t)分别为0.023和0.013外,其他值都小于0.003。验证了稀疏组Lasso方法通过施加组内约束,可以辨识因果影响的具体滞后信息。

表1 模型1的1 000次模拟结果比较

(二)高维随机滞后影响模型的比较

本小节模拟在序列的滞后影响期数随机改变,序列的维数也不同的情形下,稀疏组Lasso和组Lasso方法检验多维时间序列Granger因果关系的结果比较。

模型2:设观测样本来自于d维时间序列(X1,2,…,Xd),服从如下VAR(2)模型:

(12)

图2 模型2的Granger因果图

其中εi(t)～N(0,0.12),i=1,2,3,4,5。在d维VAR(2)模型中,序列Xi通过随机选择的不同阶滞后变量Xi-1(t-1)或Xi-1(t-2)依赖于序列Xi-1。真实的因果图结构见图2。

本文用均匀分布U[0,1]的随机数U控制因果关系的滞后期。对不同的维数d=10,50,100和样本量n=100,200,500,分别进行1 000次模拟实验。用第2节的BIC准则选择罚参数λ,λ1,λ2,τ,门限参数分别为τ=0.08(组Lasso方法)和τ=0.1(稀疏组Lasso方法)。表2给出了1 000次模拟中F1评分和两个滞后检验正确率的平均值。本小节用两个正确率度量评价两种方法检验因果影响滞后信息的性质。正确率L指当序列Xj对序列Xi(i≠j)因果影响的滞后变量为k时,估计的Xj(t-k)系数不为0,而Xj(t-h),h≠k,k,h=1,2,…,p的系数为0的比率;正确率S表示序列Xi自身滞后对当前值影响的滞后期数辨识的正确率,即当只有Xi(t-k)对Xi(t)有影响时,估计的Xi(t-k)系数不为0,而Xi(t-h),h≠k,k,h=1,2,…,p的系数为0的比率。

表2的结果表明,对于模型2,组Lasso和稀疏组Lasso方法的F1评分值相近。相同维数下,F1得分随样本量增加而增加;相同样本量下,F1得分随维数d增加而递减;对于样本量n=500,差别变小。表明了两种方法检验Granger因果关系的有效性及其大样本性质。在进一步辨识因果影响的滞后期数上,两种方法在正确率的结果上相差很大:组Lasso方法在d=10,n=500时最大正确率在检验自身滞后值上,只有0.487 7;而稀疏组Lasso方法基本上在0.8以上,在d=100,n=100时最小,值为0.620 5,并且在检验自身滞后信息上都在0.95以上。表明组Lasso方法只控制组间稀疏性,使同一组的变量系数都为0或都不为0,损失了因果关系的滞后信息。而稀疏组Lasso在组间稀疏的同时控制组内稀疏性,可以进一步选择不为0的组内的显著滞后变量。

表2 模型2的1 000次模拟结果比较

(三)方法的稳健性分析

Bolstad等在自回归系数满足一定的约束条件下证明了组Lasso方法变量选择的一致性。但没有模拟不同系数对结果的影响[13]。本小节考虑不同的系数对稀疏组Lasso方法辨识Granger因果联系及滞后信息的影响。仍用模型1的因果关系结构,其中的系数都取相等值,记为C。样本量n=3 000,系数C=0.1和0.05,分别进行1 000次模拟实验。表3给出了1 000次模拟中滞后期检验的正确率。系数C=0.1时,样本量n=3 000,门限参数τ=0.01,罚参数(λ1,λ2)=(0.3,0.2),1 000模拟的F1得分平均值为0.981 4;系数C=0.05时,样本量n=3 000,门限参数τ=0.005,罚参数(λ1,λ2)=(0.5,0.3),1 000模拟的F1得分平均值为0.756 8。相同样本量下,在系数C=0.1时要明显优于系数C=0.05的结果。

进一步计算因果影响在具体滞后期上的拒绝率,结果见表3。

表3 稀疏组Lasso方法不同系数的检验结果

由表3可知,与系数C=0.1相比,系数C=0.05时,检测出正确的滞后期的比率有所降低,特别是对于存在密切滞后的序列X4,X1(t-1)→X4(t),X5(t-3)→X4(t)和X4(t-1)→X4(t),分别为0.62,0.579和0.639;同时在其余两个滞后期上不为0的比例也有提高。模拟结果验证了作为一种基于惩罚约束的变量选择方法,稀疏组Lasso方法的有效性对联系的强度有一定要求。总的来说,稀疏组Lasso方法在系数较小时基本上也可以得到正确的Granger因果图结构和联系的滞后信息。

四、宏观经济数据实证分析

本节将稀疏组Lasso Granger因果图方法应用到中国宏观经济序列,分析序列间的直接因果关系。

(一)数据预处理

考虑到数据的完整性,选取8个宏观经济指标,包括消费指标:社会消费品零售总额X1;财政指标:国家财政支出X2和国家财政收入X3;对外贸易指标:出口总值X4和进口总值X5;交通运输指标:货运量X6;投资指标:房地产投资X7和固定资产投资完成额X8;选取2000年1月到2017年12月的月度观测值,共216个样本。数据来自于国家统计局http:∥data.stats.gov.cn/。

首先检验原始数据的平稳性。由于季节和趋势因素,原始序列是不平稳的。进行一阶差分运算消除长期趋势,ADF检验表明,8个变量一阶差分序列是平稳的,即都为1阶单整序列。然后对序列进行12步季节差分消除周期趋势,处理后的变量仍然用X1,X2,…,X8表示,预处理后的样本量n=203。

稀疏组Lasso检验Granger因果关系以VAR模型为基础,首先需要确定VAR模型的滞后阶数p。表4的结果表明,SC准则选择的最优滞后阶数为p=1,HQ准则选择的最优滞后阶数为p=2,本文选择VAR(2)模型。

表4 VAR模型定阶结果

图3 8个宏观经济变量的稀疏组Lasso图

(二)稀疏组Lasso Granger因果图

将稀疏组lasso方法式(6)用于VAR(2)模型,估计参数,按照本文提出的建立稀疏组Lasso Granger因果图模型的算法步骤建立图模型。在实际应用中,选择一个门限值τ,约束估计值小于τ的参数为0。使BIC准则达到最小的罚参数为(0.009,0.05),门限值τ=0.1,得到12条边,模型满足稀疏性假设和对系数的要求。估计的图见图3。

图4 8个宏观经济变量的组Lasso图

由图3可知,最大自由度的变量是X1(社会消费品零售总额),其中X3(国家财政收入),X8(固定资产投资完成额)是X1的Granger原因,X1是X7(房地产投资)和X5(进口总值)的Granger原因。表明社会消费品零售总额是反映宏观经济运行状况的重要经济统计指标。其次,X5(进口总值)的自由度较大,除X1(社会消费品零售总额)外,X3(国家财政收入),X7(房地产投资)和X8(固定资产投资完成额)是其Granger原因。表明现阶段中国对外贸易依存度不断攀升,进口的供给效应在经济发展中的作用愈发明显。

另外,X2(国家财政支出)是X3(国家财政收入),X4(出口总值)和X6(货运量)的Granger原因,表明财政投资具有生产性,能够带动经济。X7(房地产投资)是X5(进口总值)的Granger原因,表明房地产对宏观经济的发展有影响[17]。X8(固定资产投资完成额)影响X1和X5,表明固定资产投资可以积极推动中国社会经济增长[18]。X4(出口总值)和X6(货运量)互为Granger原因。图3还表明X2(国家财政支出)和X8(固定资产投资完成额)在经济运行中是主要变量。

将稀疏组Lasso方法得到的图模型(图3)与组Lasso方法得到的图模型(图4)相比较,组Lasso Granger因果图少了X2(国家财政支出)对X6(货运量)、X1(社会消费品零售总额)对X7(房地产投资)、X4(出口总值)对X6(货运量)的Granger影响;并且稀疏组Lasso方法的最小BIC值(6.283 6)小于组Lasso方法的最小BIC值(6.359 1),进一步表明了稀疏组Lasso方法的优越性。稀疏组Lasso图的滞后信息矩阵如下(其中“-”表示无联系):组Lasso方法只能检验Gragner因果关系的存在性,而不能确定因果关系是由哪一个滞后变量的影响导致的,在应用Granger因果关系时,确定因果影响的具体滞后变量对实现精确的管理和决策非常重要。而稀疏组Lasso方法可以通过控制组内稀疏性,检测因果影响的滞后阶数。

(13)

五、结 论

本文通过将多维时间序列的Granger因果关系检验问题转化为变量选择问题,提出用稀疏组Lasso方法检验Grange因果关系的存在性和因果影响的滞后期,建立稀疏组Lasso Granger因果图模型直观分析多维序列间的因果相依结构。与组Lasso方法相比,稀疏Lasso方法在对变量按照不同分量序列分组的基础上,通过组间约束检验Granger因果关系,用组内约束在不为0的组中进一步选择滞后变量,保留了Granger因果影响的滞后信息。数值模拟验证了提出的稀疏的Lasso方法在检验Granger因果关系及其滞后信息上的有效性。最后将方法应用到宏观经济变量数据,建立8个变量间的系数Granger因果图模型,分析因果关系及其影响的滞后阶数。

本文主要考虑了宏观经济数据的应用,其他如基因数据等也可以应用本文方法进行分析。非线性因果关系如动态Granger因果关系的图模型方法也是下一步的研究方向。