魏岳嵩（淮北師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院，安徽淮北235000；西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西西安710072）

非線性結(jié)構(gòu)向量自回歸因果圖的廣義似然比辨識(shí)方法

魏岳嵩
（淮北師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院，安徽淮北235000；西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西西安710072）

利用圖模型方法研究非線性結(jié)構(gòu)向量自回歸模型的因果性問(wèn)題.構(gòu)建了非線性結(jié)構(gòu)向量自回歸因果圖模型，提出圖模型因果性的廣義似然比辨識(shí)方法.構(gòu)造同期因果關(guān)系和滯后因果關(guān)系的廣義似然比統(tǒng)計(jì)量，使用bootstrap方法來(lái)確定檢驗(yàn)統(tǒng)計(jì)量的原分布，模擬研究論述了方法的有效性.

非線性結(jié)構(gòu)向量自回歸模型；因果圖；條件獨(dú)立；廣義似然比；bootstrap方法

O211.6

1　引言

檢驗(yàn)和識(shí)別變量間的因果關(guān)系是時(shí)間序列分析的重要問(wèn)題之一.近年來(lái)，圖模型方法已經(jīng)成為分析和處理時(shí)間序列問(wèn)題的有力工具［1-3］.利用搜索算法特別是PC算法分析變量間的因果關(guān)系已經(jīng)成為當(dāng)前常用的方法［4-6］.然而利用此類(lèi)方法分析因果關(guān)系存在兩個(gè)主要問(wèn)題，即要求時(shí)間序列模型是線性的且噪聲項(xiàng)服從高斯分布.雖然線性及高斯性假設(shè)能使問(wèn)題簡(jiǎn)單化，但研究表明非線性及非高斯性可能會(huì)更貼合實(shí)際的數(shù)據(jù)生成過(guò)程.

對(duì)于非線性時(shí)間序列的圖模型問(wèn)題，雖然Chu和G lym our［7］以及W ei和T ian［8］等分別利用可加模型回歸法和信息論的方法獲得一些有益的結(jié)論，但當(dāng)前，建立有效的非線性時(shí)間序列因果關(guān)系的圖模型方法仍是尚未解決的問(wèn)題.本文主要討論如何利用圖模型方法來(lái)研究非線性結(jié)構(gòu)向量自回歸模型的因果性問(wèn)題.

2　非線性結(jié)構(gòu)向量自回歸因果圖模型

設(shè)V是一非空有限集，圖G=（V，E）是一有序集，其中V中的元素稱(chēng)為圖的頂點(diǎn)，而E=｛（a，b）|a，b∈V｝中的元素稱(chēng)為邊.若有序?qū)Γ╝，b）∈E而（b，a）/∈E，則稱(chēng)頂點(diǎn)a和b之間存在由a到b的有向邊，記為a→b，此時(shí)稱(chēng)a為b的父親（parent），a的父親集表示為pa（a）.長(zhǎng)度為n的從a到b的路徑是指由不同頂點(diǎn)組成的從a到b的序列｛a=i0，···，in=b｝，滿足對(duì)所有的k= 1，···，n，都有（ik-1，ik）∈E.如果對(duì)于所有的k=1，···，n，都有（ik-1，ik）∈E但（ik，ik-1）/∈E，則稱(chēng)該路徑為有向路徑.若在圖G中存在由a到b的有向路徑，則稱(chēng)b是a的后代.a的所有后代組成的集合稱(chēng)為a的后代集，記為de（a）.

設(shè)n維時(shí)間序列｛X（t）=（X1，t，···，Xn，t）T，-∞＜t＜∞｝由p階非線性結(jié)構(gòu)向量自回歸模型（Non linear Structu ral Vector Au toregressive M odels，簡(jiǎn)記為NSVAR）生成，即

其中εit，i=1，2···n是相互獨(dú)立的噪聲項(xiàng)，fj，i（·）和gk，i，l（·）都是一元光滑函數(shù)（可以是非線性的），且不存在序列j1，j2，···，jm（m≤n）使得fj1j2，fj2j3，···，fjm-1jm，fjmj1都不為零，存在k和i使得gk，i，p（·）/=0.

NSVAR模型可以利用因果關(guān)系加以描述，即在模型（1）中，若fj，i（·）/=0，則稱(chēng)Xj，t是Xi，t的直接原因，并稱(chēng)Xj，t和Xi，t之間存在同期因果關(guān)系；若gk，i，l（·）/=0，則稱(chēng)Xk，t-l是Xi，t的直接原因，并稱(chēng)Xk，t-l和Xi，t之間存在滯后因果關(guān)系.此外，由于同期變量間不存在反饋關(guān)系，因此該模型能用有向非循環(huán)圖加以解釋.

定義2.1（非線性結(jié)構(gòu)向量自回歸因果圖）對(duì)于p階NSVAR模型（1），記

稱(chēng)有向非循環(huán)圖G=（V，E）為非線性結(jié)構(gòu)向量自回歸因果圖（Non linear Structu ral Vector Autoregressive Model Causal G raphs，簡(jiǎn)記為NSVARCG），如果圖G滿足：

定義2.2（因果M arkov條件）對(duì)于p階NSVAR模型（1），記V=｛X1t，···，Xnt，X1，t-1，···，Xn，t-1，···，X1，t-p，···，Xn，t-p｝，有向非循環(huán)圖G=（V，E）是NSVAR模型（1）的因果圖，則圖G滿足因果M arkov條件，即對(duì)V的任意子集A總有A⊥V（de（A）∪pa（A））|pa（A），其中符號(hào)⊥表示獨(dú)立，pa（A）及de（A）分別表示不包含集合A中元素的那些A中元素的父親集和后代集合：pa（A）=∪i∈Apa（i）A，de（A）=∪i∈Ade（i）A.

3　同期因果關(guān)系的廣義似然比辨識(shí)法

定理3.1設(shè)圖G是模型（1）的因果圖，Xt=（X1t，···，Xnt）T，Xit和Xjt是Xt的任意兩元素，誘導(dǎo)得到的G的子圖，則Xit和Xjt在G中被Xct和Xld-分離當(dāng)且僅當(dāng)Xit和Xjt在G1中被Xctd-分離.

證設(shè)Xit和Xjt在G中被Xct和Xld-分離，則若π是G中Xit和Xjt之間的一條僅包含Xt中點(diǎn)的路徑，π一定是阻塞路徑，因此G1中Xit和Xjt之間的任意路徑都是阻塞路徑，即Xit和Xjt在G1中被Xctd-分離.反之，當(dāng)Xit和Xjt在G1中被Xctd-分離，則如果π是G中一條Xit和Xjt之間的相對(duì)于Xct的d-連通路徑，π一定包含Xl中的點(diǎn).由于不存在由Xt中點(diǎn)指向Xl中點(diǎn)的有向邊，因此π中屬于Xl的點(diǎn)都為非對(duì)沖點(diǎn)，從而該路徑相對(duì)于Xct和Xl是阻塞的，故Xit和Xjt在G中被Xct和Xld-分離.

注定理說(shuō)明當(dāng)前變量之間的因果關(guān)系不受滯后變量的影響，因此在研究當(dāng)前變量之間的因果關(guān)系時(shí)，分離集只需考慮當(dāng)前變量集合的子集，而無(wú)需考慮滯后變量的影響，由此縮小了分離集的研究范圍，降低模型結(jié)構(gòu)辨識(shí)的難度.

設(shè)在NSVARCG圖中，Xit和Xjt是任意兩個(gè)同期變量，由于同期變量之間不存在反饋關(guān)系，若二者存在因果關(guān)系則必為Xit→Xjt或者Xjt→Xit，分別對(duì)應(yīng)于兩個(gè)假設(shè)：

因此，可以將確定Xit和Xjt之間的因果性方向轉(zhuǎn)化為上述的假設(shè)檢驗(yàn)問(wèn)題.這是一個(gè)非參數(shù)模型對(duì)非參數(shù)模型的檢驗(yàn)問(wèn)題.這里采用Fan等［9］提出的廣義似然比（Generalized likelihood ratio，簡(jiǎn)記為GLR）方法處理上述的假設(shè)檢驗(yàn)問(wèn)題.為此，首先考慮如下兩個(gè)簡(jiǎn)單的假設(shè)檢驗(yàn)問(wèn)題：

和

其中

將上式關(guān)于參數(shù)σ2最大化可得備擇假設(shè)下的似然：

同理可得到在原假設(shè)下的對(duì)數(shù)似然為：

其中

由此可定義GLR統(tǒng)計(jì)量為：

從而對(duì)于同期變量因果關(guān)系定向的假設(shè)檢驗(yàn)問(wèn)題H0?H1，構(gòu)造GLR統(tǒng)計(jì)量：

為了避免復(fù)雜的漸近方差的計(jì)算，這里利用條件bootstrap方法來(lái)計(jì)算GLR統(tǒng)計(jì)量Tij的p值，具體步驟如下：

步驟2計(jì)算GLR檢驗(yàn)統(tǒng)計(jì)量Tl，Tij及備擇模型下的殘差項(xiàng)

步驟3對(duì)于原始數(shù)據(jù)Xk，用中心化的殘差序列｛來(lái)生成其中是的平均值.

步驟4用這組條件bootstrap樣本來(lái)構(gòu)造GLR統(tǒng)計(jì)量和；

步驟6利用經(jīng)驗(yàn)分布

4　滯后因果關(guān)系的廣義似然比辨識(shí)法

采取和上節(jié)相同的方法，構(gòu)造GLR統(tǒng)計(jì)量：

其中

這里同樣采用條件條件bootstrap方法來(lái)計(jì)算統(tǒng)計(jì)量T的p值.檢驗(yàn)步驟如下：

步驟1由局部多項(xiàng)式估計(jì)法擬合原假設(shè)H0和備擇假設(shè)H1下的方程模型；

步驟2計(jì)算GLR檢驗(yàn)統(tǒng)計(jì)量T及備擇模型下的殘差項(xiàng)ε?ki，k=1，2，···，m；

步驟3對(duì)于原始數(shù)據(jù)Xk，用中心化的殘差序列來(lái)生成構(gòu)造條件bootsrap樣其中f?和是在步驟1中原假設(shè)下所估計(jì)的回歸函數(shù)；

步驟4用這組條件bootstrap樣本來(lái)構(gòu)造GLR統(tǒng)計(jì)量；

步驟6利用經(jīng)驗(yàn)分布

作為T(mén)在原假設(shè)下分布的近似分布.對(duì)選擇的顯著性水平α，構(gòu)造拒絕域（-∞）和（C1-，+∞），其中Cα表示分布B（x）的α分位點(diǎn)；

5　模擬分析

為了驗(yàn)證所給方法的有效性，對(duì)以下三個(gè)模型進(jìn)行數(shù)值模擬，三個(gè)模型分別選取為線性結(jié)構(gòu)向量自回歸模型、非線性可加模型和非線性結(jié)構(gòu)向量自回歸模型.同時(shí)為了加以比較，也分別利用PC算法、文［7］提出的可加自回歸方法（簡(jiǎn)記為ARM方法）以及文［8］提出的信息論方法（簡(jiǎn)記為ITM方法）對(duì)所給模型做模擬分析.

模型1：

模型2：

模型3：

其中εit～N（0，1），i=1，2，3是相互獨(dú)立的噪聲項(xiàng).

對(duì)于每一個(gè)因果模型，生成樣本容量為1000的模擬序列，并分別利用ITM方法，GLR方法，PC算法和ARM方法辨識(shí)原模型的因果結(jié)構(gòu).為了研究方便，所有的滯后階數(shù)選取為2，其中在應(yīng)用ITM方法時(shí)，帶寬選取∈=1.0（帶寬選擇具體參見(jiàn)文［8］）.在應(yīng)用GLR方法時(shí)，分別采取局部多項(xiàng)式估計(jì)法和局部線性估計(jì)法構(gòu)造GLR統(tǒng)計(jì)量.由于高維空間中局部數(shù)據(jù)的稀疏性，因此局部多項(xiàng)式估計(jì)時(shí)采用局部二次多項(xiàng)式擬合方法確定GRL統(tǒng)計(jì)量，條件bootstrap方法中B=99，顯著性水平α=0.1，實(shí)驗(yàn)運(yùn)行1000次.表1至表9給出了所得模擬結(jié)果，其中同期因果關(guān)系表中符號(hào)（以（Xt，Yt）為例）??表示Xt和Yt之間不存在邊，?-?表示Xt和Yt之間存在邊但無(wú)法定向，?→?和?←?分別表示因果關(guān)系Xt→Yt和Xt←Yt.

表1　G LR算法所得模型同期因果關(guān)系比率

表2　ITM算法所得模型同期因果關(guān)系比率

表3　PC算法所得模型同期因果關(guān)系p值

表4　ARM算法所得模型同期因果關(guān)系比率

表5　ITM算法所得模型滯后因果關(guān)系比率

表6　GLR算法所得模型滯后因果關(guān)系比率（局部多項(xiàng)式估計(jì)）

表7　GLR算法所得模型滯后因果關(guān)系比率（局部線性估計(jì)）

表8　PC算法所得模型滯后因果關(guān)系比率

表9　ARM算法所得模型滯后因果關(guān)系比率

表1-表4給出了不同算法所得到的三個(gè)因果模型的同期因果關(guān)系.模擬結(jié)果顯示，除了模型1之外，ITM算法對(duì)于模型2和模型3都能得到理想的結(jié)果.對(duì)于模型1，ITM算法雖然能正確的捕捉到所有同期變量之間的關(guān)系，但無(wú)法確定因果方向，造成這一結(jié)果的主要原因在于ITM算法主要利用條件互信息來(lái)給出變量間的因果關(guān)系，而僅僅由條件I（Xt，Zt|Yt）≤I（Xt，Zt）是無(wú)法區(qū)分結(jié)構(gòu)Xt→Yt→Zt和Zt→Yt→Xt，事實(shí)上此時(shí)得到的是等價(jià)結(jié)構(gòu)關(guān)系.表5則給出的是由ITM算法所得到的三個(gè)模型的滯后因果關(guān)系.ITM算法能正確的辨識(shí)原模型中幾乎所有的滯后因果關(guān)系，只是在模型1中，關(guān)于Xt-1→Yt和Yt-1→Zt的模擬結(jié)果稍顯偏大.

表1，表6和表7的模擬結(jié)果顯示，對(duì)于三個(gè)模型，GLR方法都能正確的確定同期變量之間的因果關(guān)系.對(duì)于三個(gè)模型中所蘊(yùn)含的滯后因果關(guān)系，相比于ITM方法，GLR方法雖然能正確的確定原模型中幾乎所有的滯后因果關(guān)系，但也產(chǎn)生過(guò)多的因果關(guān)系.如基于局部線性估計(jì)的GLR方法所得模型2的結(jié)果以及基于局部多項(xiàng)式估計(jì)的GLR方法所得模型1的結(jié)果.從所有的結(jié)果可以看出，當(dāng)采用局部多項(xiàng)式估計(jì)時(shí)，對(duì)于非線性模型2和模型3，GLR方法能夠得到更為滿意的結(jié)果.然而對(duì)于線性模型1，GLR方法得到的結(jié)果較差，可能原因在于構(gòu)造GLR統(tǒng)計(jì)量時(shí)，采用的是局部二次多項(xiàng)式擬合原模型，從而對(duì)于線性模型1可能會(huì)造成更大的偏差.采用局部線性估計(jì)法得到的GLR方法所得結(jié)果情況則恰恰相反.這說(shuō)明在使用GLR方法辨識(shí)變量間的因果關(guān)系時(shí)，模擬的結(jié)果也依賴于構(gòu)造的GLR統(tǒng)計(jì)量時(shí)所使用的估計(jì)法.因此，如果在試驗(yàn)前能夠獲取原模型一定的先驗(yàn)知識(shí)（如線性模型或是非線性模型），采用更為合理的估計(jì)方法，相應(yīng)的會(huì)提高GLR方法所得結(jié)果的精確度.

表3和表8給出的是PC算法所得結(jié)果.模擬結(jié)果顯示PC算法僅僅對(duì)于線性時(shí)間序列模型1所得到結(jié)果較好，而對(duì)于非線性可加模型2和非線性結(jié)構(gòu)向量自回歸模型3，PC算法得到的結(jié)果明顯較差，由PC算法得到的模型2和3的因果結(jié)構(gòu)明顯缺失了較多的邊，如模型2中的Xt-1→Xt和Yt-1→Yt，模型3中的Yt-1→Yt和Zt-1→Zt等.此結(jié)論也進(jìn)一步證實(shí)了模型是線性模型與否對(duì)PC算法所得結(jié)果有著顯著的影響.

表4和表9給出的是ARM方法所得結(jié)果.從模擬結(jié)果可以發(fā)現(xiàn)，由ARM算法得到的模型1和模型2的因果結(jié)構(gòu)類(lèi)似于ITM算法得到的結(jié)果.然而ARM算法得到模型3的結(jié)構(gòu)和原模型由較大出入，它錯(cuò)誤的去除了原模型中的許多因果關(guān)系，如Yt→Xt和Yt→Zt，造成這一現(xiàn)象的原因是因?yàn)槟Ｐ?中同期變量間的因果關(guān)系是非線性的，違反了ARM算法使用的前提條件.

［1］Oxley L，Reale M，W ilson G.Constructing structural VAR models w ith conditional independence graphs［J］.M athem atics and Com pu ters in Simu lation，2009，79:2910-2916.

［2］Eich lerM.Graphicalmodelling ofmu ltivariate time series［J］.Probability Theory and Related Fields，2012，153:233-268.

［3］Gao Wei，Tian Zheng.Latent ancestral graph of structure vector autoregressivemodels［J］. Journal of System s Engineering and Electronics，2010，21:233-238.

［4］Hoover K.Automatic inference of the contem poraneous causal order of a system of equations［J］.Econom etric Theory，2005，21:69-77.

［5］M oneta A.G raphical causalm odels and VARs:an em p irical assessm ent of the real business cycles hypothesis［J］.Em pirical Econom ics，2008，35（2）:275-300.

［6］Wei Yuesong，Tian Zheng，Xiao Yanting.Learning mu ltivariate time series causal graphs based on cond itionalmu tual in form ation［J］.Jou rnal of System s Science and System s Engineering，2013，22:38-51.

［7］Chu Tian jiao，G lymour C.Search for additivenonlinear time series causalmodels［J］.Journal of M achine Learning Research，2008，9:67-91.

［8］Wei Yuesong，Tian Zheng，Xiao Yanting.Learning causal graphs of non linear structural vector autoregressivemodelsusing information theory criteria［J］.Journal of System s Science &Com p lexity，2014，27（6）:1213-1226.

［9］Fan Jianqin，Zhang Chunm ing，Zhang Jian.Generalized likelihood ratio statistics and W ilks phenomenon［J］.The Annals of Statistics，2001，29（1）:153-193.

M R Su b jec t C lassifica tion:62M 10

Generalized likelihood ratio app roach fo r iden tifying non linear structu ral vector au toregressive causal graphs

WEIYue-song
（School of M athem atical Science，Huaibei Norm al University，Huaibei 235000，China；Department of App lied M athematics，Northwest Polytechnical University，Xi’an 710129，China）

In this paper，the causal relationships among variables of nonlinear structure vector autoregressivemodel are studied using graphicalmodelmethod.The non linear structure vector autoregressive causal graph is p resented，and a generalized likelihood ratio app roach is developed to in fer the causal relationships.The generalized likelihood ratio statistics of contem poraneous and lagged are presen ted respectively，and a bootstrap m ethod is considered for determ ining the nu ll distribu tion of the test statistic.The validity of the p roposed method is con firmed by simu lations analysis.

nonlinear structu ral vector autoregressivem odels；graphicalm odels；conditional independence；generalized likelihood ratio；bootstrap method

A

1000-4424（2016）02-0143-10

2015-01-28

2016-03-13

國(guó)家自然科學(xué)基金（60375003；61201323）；安徽省高校自然科學(xué)研究項(xiàng)目（K J2015A 035）