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

非線性結構向量自回歸因果圖的廣義似然比辨識方法

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

利用圖模型方法研究非線性結構向量自回歸模型的因果性問題.構建了非線性結構向量自回歸因果圖模型，提出圖模型因果性的廣義似然比辨識方法.構造同期因果關系和滯后因果關系的廣義似然比統(tǒng)計量，使用bootstrap方法來確定檢驗統(tǒng)計量的原分布，模擬研究論述了方法的有效性.

非線性結構向量自回歸模型；因果圖；條件獨立；廣義似然比；bootstrap方法

O211.6

1　引言

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

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

2　非線性結構向量自回歸因果圖模型

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

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

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

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

定義2.1（非線性結構向量自回歸因果圖）對于p階NSVAR模型（1），記

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

定義2.2（因果M arkov條件）對于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條件，即對V的任意子集A總有A⊥V（de（A）∪pa（A））|pa（A），其中符號⊥表示獨立，pa（A）及de（A）分別表示不包含集合A中元素的那些A中元素的父親集和后代集合：pa（A）=∪i∈Apa（i）A，de（A）=∪i∈Ade（i）A.

3　同期因果關系的廣義似然比辨識法

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

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

注定理說明當前變量之間的因果關系不受滯后變量的影響，因此在研究當前變量之間的因果關系時，分離集只需考慮當前變量集合的子集，而無需考慮滯后變量的影響，由此縮小了分離集的研究范圍，降低模型結構辨識的難度.

設在NSVARCG圖中，Xit和Xjt是任意兩個同期變量，由于同期變量之間不存在反饋關系，若二者存在因果關系則必為Xit→Xjt或者Xjt→Xit，分別對應于兩個假設：

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

和

其中

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

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

其中

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

從而對于同期變量因果關系定向的假設檢驗問題H0?H1，構造GLR統(tǒng)計量：

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

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

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

步驟4用這組條件bootstrap樣本來構造GLR統(tǒng)計量和；

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

4　滯后因果關系的廣義似然比辨識法

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

其中

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

步驟1由局部多項式估計法擬合原假設H0和備擇假設H1下的方程模型；

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

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

步驟4用這組條件bootstrap樣本來構造GLR統(tǒng)計量；

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

作為T在原假設下分布的近似分布.對選擇的顯著性水平α，構造拒絕域（-∞）和（C1-，+∞），其中Cα表示分布B（x）的α分位點；

5　模擬分析

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

模型1：

模型2：

模型3：

其中εit～N（0，1），i=1，2，3是相互獨立的噪聲項.

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

表1　G LR算法所得模型同期因果關系比率

表2　ITM算法所得模型同期因果關系比率

表3　PC算法所得模型同期因果關系p值

表4　ARM算法所得模型同期因果關系比率

表5　ITM算法所得模型滯后因果關系比率

表6　GLR算法所得模型滯后因果關系比率（局部多項式估計）

表7　GLR算法所得模型滯后因果關系比率（局部線性估計）

表8　PC算法所得模型滯后因果關系比率

表9　ARM算法所得模型滯后因果關系比率

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

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

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

表4和表9給出的是ARM方法所得結果.從模擬結果可以發(fā)現(xiàn)，由ARM算法得到的模型1和模型2的因果結構類似于ITM算法得到的結果.然而ARM算法得到模型3的結構和原模型由較大出入，它錯誤的去除了原模型中的許多因果關系，如Yt→Xt和Yt→Zt，造成這一現(xiàn)象的原因是因為模型3中同期變量間的因果關系是非線性的，違反了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

國家自然科學基金（60375003；61201323）；安徽省高校自然科學研究項目（K J2015A 035）