王 慧, 吳茜茜

(合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230601)

0 引 言

空間自回歸(spatial autoregression,SAR)模型是空間計量模型的一個特例,常被用于分析具有空間相關(guān)性的數(shù)據(jù)。由于其獨特的博弈結(jié)構(gòu)并且可以被解釋為反應(yīng)函數(shù),被廣泛地應(yīng)用于空間計量經(jīng)濟學(xué)和社交網(wǎng)絡(luò)建模中[1]。建立模型的過程離不開對模型中未知參數(shù)的估計,而針對SAR模型中參數(shù)的估計,文獻(xiàn)[2]提出編碼法;文獻(xiàn)[3]利用最小二乘估計方法研究SAR模型;文獻(xiàn)[4]概括了SAR模型的似然函數(shù),并給出相應(yīng)的極大似然估計,但SAR模型的極大似然估計量的收斂速度取決于空間權(quán)重矩陣的特征[5];文獻(xiàn)[6]提出SAR模型的有限信息的極大似然估計量,但由于空間相關(guān)性,得到的是非一致估計。隨著廣義矩估計方法的提出,很多學(xué)者將其應(yīng)用到空間計量模型中[7-8],但在異方差存在的情況下,該估計是非一致的[9]。

隨著貝葉斯理論的發(fā)展,越來越多的學(xué)者開始重視貝葉斯統(tǒng)計推斷這一有力工具。文獻(xiàn)[10-12]研究了SAR模型和受限因變量SAR模型的貝葉斯方法,突破了SAR模型經(jīng)典方法的限制;文獻(xiàn)[13]的研究結(jié)果表明,當(dāng)SAR模型中存在負(fù)的空間相關(guān)性時,通過馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法取得的貝葉斯估計量比穩(wěn)健的廣義矩估計量更有效。實際上,貝葉斯統(tǒng)計推斷主要依靠后驗分布,為此需要計算似然函數(shù)并給出合理的先驗[14]。當(dāng)模型參數(shù)空間維數(shù)增加或先驗結(jié)構(gòu)復(fù)雜時,似然函數(shù)越來越復(fù)雜而無法獲取其解析表達(dá)式,無法計算后驗密度,此時不能直接運用貝葉斯方法。近似貝葉斯計算(approximate Bayesian computation,ABC)最早由文獻(xiàn)[15]提出并應(yīng)用于種群基因的研究。與經(jīng)典貝葉斯方法相比,ABC方法可以更合理地利用先驗,作為一種“似然自由”的方法解決實際問題。目前ABC方法在生物及其他領(lǐng)域應(yīng)用廣泛,但在經(jīng)濟領(lǐng)域的應(yīng)用偏少,本文將ABC方法引入SAR模型的推斷分析。文獻(xiàn)[16]將ABC方法應(yīng)用于量子核磁共振模型推斷;文獻(xiàn)[17]將ABC方法應(yīng)用于含噪聲數(shù)據(jù)的隨機伽馬過程;文獻(xiàn)[18]通過ABC方法和穩(wěn)態(tài)信號模擬,從轉(zhuǎn)錄組和蛋白質(zhì)組數(shù)據(jù)中的逆向工程指導(dǎo)基因調(diào)控網(wǎng)絡(luò)。

本文首先介紹貝葉斯推斷的思想、常見的馬爾科夫鏈蒙特卡洛、近似貝葉斯計算方法、空間自回歸模型;其次提出適應(yīng)估計該模型參數(shù)的2種算法;最后基于長三角地區(qū)26個城市2007-2016年間的空間面板數(shù)據(jù),建立相應(yīng)的SAR模型,根據(jù)相應(yīng)算法對該模型進(jìn)行估計,并給出結(jié)果和討論。

1 貝葉斯推斷

1.1 貝葉斯定理

貝葉斯定理是貝葉斯推斷的核心,最早由英國數(shù)學(xué)家Thomas Bayes提出[14]。模型分析中,通常利用已知數(shù)據(jù)(記為Y)與模型形式估計未知參數(shù)(記為θ)。

根據(jù)貝葉斯公式,有

(1)

其中:π(θ)為參數(shù)θ的先驗信息;L(Y|θ)為似然函數(shù);π(θ|Y)為參數(shù)θ的后驗分布。

1.2 馬爾可夫鏈蒙特卡洛(MCMC)方法

MCMC方法是貝葉斯統(tǒng)計推斷中一種重要的抽樣方法。它通過構(gòu)建一條馬爾可夫鏈?zhǔn)蛊淦椒€(wěn)分布為后驗分布,由此產(chǎn)生的該鏈上的樣本點{θ(1),θ(2),…,θ(n)}可以看作近似服從后驗分布π(θ|Y)。該方法的關(guān)鍵是如何構(gòu)造“合適”的馬爾科夫鏈,常見的算法有M-H(Metropolis-Hastings)抽樣和Gibbs抽樣。

1.2.1M-H抽樣

選取建議分布q(·),通常可以采用對稱式(q(x|y)=q(y|x),q(x|y)=q(|x-y|)),獨立式(q(x|y)=q(x))或隨機游走的形式(y=x+η)。具體步驟如下:

(1) 令t=0,由π(θ)抽樣θ(0)。

(2) 由q(·)抽樣θ*。

(3) 計算接受概率α。計算公式為：

(4) 由U(0,1)抽樣u,若u≤α,則θ(t+1)=θ*,否則θ(t+1)=θ(t)。

(5) 令t=t+1,重復(fù)上述步驟直到t=N。

1.2.2Gibbs抽樣

上述M-H抽樣針對低維問題較為有效。如果待估參數(shù)為向量θ=[θ1θ2…θn]T,其中T表示向量的轉(zhuǎn)置,并且可以計算參數(shù)θi的后驗條件密度π(θi|θ-i),其中θ-i=[θ1…θi-1θi+1…θn]T是由參數(shù)向量θ剔除了參數(shù)θi所構(gòu)成的,那么可以采用Gibbs抽樣,具體步驟如下:

(3) 令t=t+1,重復(fù)上述步驟直到t=N。特別地,當(dāng)其中某個后驗條件分布難以直接采樣時,可以利用M-H抽樣針對該分布進(jìn)行抽樣,由此構(gòu)成Metropolis within Gibbs算法。

1.3 近似貝葉斯計算(ABC方法)

ABC方法避免了(1)式中似然函數(shù)L(Y|θ)的計算,其核心思想是從先驗分布中直接抽樣,通過距離函數(shù)等條件,篩選參數(shù)集來近似估計參數(shù)后驗分布[15]。首先需要定義距離函數(shù)S(·)和容差閾值d0,ABC方法(通常又稱為拒絕算法)的算法步驟如下:

(1) 由π(θ)抽樣θ*。

(2) 基于模型和參數(shù)θ*計算模擬值Y′。

(3) 通過S(·)計算Y與Y′之間的距離d。

(4) 若d≤d0,則接受θ*;反之,則拒絕θ*。

對于比較小的d0,最終獲得的模擬樣本分布π(θ|S(Y,Y′)≤d0)可用于近似后驗分布π(θ|Y)[15]。

1.4 SAR模型

SAR模型的一般形式為:

Y=ρWY+Xβ+ε,ε～MVN(0,σ2In)

(2)

其中:Y為(n×1)觀測值;X為(n×k)解釋變量矩陣;W為空間權(quán)重矩陣;β為(k×1)系數(shù);ρ為空間自回歸系數(shù);In為n階單位矩陣;ε=[ε1ε2…εn]T表示隨機誤差項,服從MVN(0,σ2In)。模型(2)的待估參數(shù)向量為θ=(β,σ2,ρ)T,且假設(shè)獨立性存在,似然函數(shù)為:

(3)

其中:f(·)=exp(-(2σ2)-1(A(·)Y-Xβ)T·(A(·)Y-Xβ));A(ρ)=In-ρW。

2 算 法

2.1 Metropolis within Gibbs算法

假設(shè)參數(shù)β、σ2、ρ相互獨立,先驗分別為π(β)～N(μ、Σ)、π(σ2)～I(xiàn)G(α,γ)、π(ρ)～U(a,b)[12]。

由獨立性及似然的表達(dá)可以得出后驗的表達(dá)式,可利用M-H抽樣進(jìn)行采樣,但由于待估參數(shù)維度高且參數(shù)類型不同,采用Gibbs抽樣更加合適。根據(jù)Gibbs抽樣,在后驗條件分布抽樣過程中,第t+1步迭代為:

β(t+1)～π(β|ρ(t),(σ2)(t),Y,W)

(4)

(σ2)(t+1)～π(σ2|ρ(t),β(t+1),Y,W)

(5)

ρ(t+1)～π(ρ|ρ(t),β(t+1),(σ2)(t+1),Y,W)

(6)

其中:β、σ2的條件后驗分別為正態(tài)分布、逆伽瑪分布,可以直接抽樣。但是ρ的條件后驗形式比較復(fù)雜,這里選擇應(yīng)用獨立的M-H抽樣對ρ的條件后驗抽樣,由此得到的Metropolis within Gibbs算法的步驟如下:

(1) 令t=0,由先驗分布π(β)、π(σ2)、π(ρ)分別抽樣θ(0)=[β(0)(σ2)(0)ρ(0)]T。

(2) 由(4)式、(5)式分別抽樣β(t+1)、(σ2)(t+1)。

(3) 由建議分布ρ*=c·η(η～N(0,1))抽樣ρ*。

(4) 計算接受概率。計算公式為：

由U(0,1)抽樣u,若u≤α,則ρ(t+1)=ρ*;否則ρ(t+1)=ρ(t)。

(5) 令t=t+1,重復(fù)上述步驟直到t=N。

2.2 ABC方法

根據(jù)ABC方法,需要通過模型計算距離函數(shù)S(Y,Y′)。

模型(2)中存在誤差項ε,且ε～MVN(0,σ2In),只要控制住了誤差項ε(即εi,i=1,2,…,n),就可以近似計算出模擬數(shù)值Y′。

因為εi在(±3σ)中的概率為99.74%,所以有99.74%的把握將εi控制在一個區(qū)間內(nèi),只要計算出εi的區(qū)間端點處的值,然后代入計算即可。

在SAR模型中抽取參數(shù)θ*,模擬數(shù)值Y′計算公式為:

Y′=(In-ρ*W)-1(Xβ*+ε*)

(8)

其中:ε*=max(-3σ*In, 3σ*In)。具體的算法步驟如下:

(1) 令t=0,由先驗分布π(β)、π(σ2)和π(ρ)分別抽樣β*、(σ2)*和ρ*,記為θ*。

(2) 根據(jù)(8)式計算模擬數(shù)據(jù)Y′。

(3) 計算距離d(采用切比雪夫距離)。計算公式為：

若d≤d0,則θ(t+1)=θ*;否則返回步驟1。

(4) 令t=t+1,重復(fù)上述步驟，直到t=N。

3 模型建立及數(shù)值模擬

3.1 模型假設(shè)

本文使用長三角26個城市2006-2017年的面板數(shù)據(jù)進(jìn)行實際論證。數(shù)據(jù)分別來自《中國城市統(tǒng)計年鑒》及江蘇、浙江、安徽省統(tǒng)計年鑒,對于缺失1年的樣本采用均值法補充,對于缺失2年及以上的數(shù)據(jù)采用線性插值法補充。假定社會投入和產(chǎn)出符合Cobb-Douglas生產(chǎn)函數(shù),通過計算這26個城市的地區(qū)生產(chǎn)總值(per capita regional gross domestic product, PGDP)的Moran’s I指數(shù)及其統(tǒng)計量,結(jié)果表明,利用該數(shù)據(jù)建立SAR模型是可行的[18]。通過Pearson相關(guān)性檢驗研究變量之間的相關(guān)性,對于存在多重共線性的變量保留其中一個,保證估計值的唯一性。本文選擇解釋變量為服務(wù)業(yè)集聚度(agglomeration, AGG)、固定資本投入(fixed asset input, FI)、外商直接投資(foreign direct investment, FDI)以及勞動力資本投入(Human)。

為便于研究各要素投入對生產(chǎn)地區(qū)增長的貢獻(xiàn)率,消除異方差和數(shù)據(jù)的劇烈波動,建立對數(shù)化SAR模型如下:

(9)

其中:i=1,2,3,4;j=1,2,…,26;Yjt為城市j在第t年的地區(qū)生產(chǎn)總值(PGDP);X1jt、X2jt、X3jt、X4jt分別為城市j在第t年的服務(wù)業(yè)集聚度(AGG)、固定資本投入(FI)、外商直接投資(FDI)、勞動力資本投入(Human);ρ為空間自回歸系數(shù);εjt～N(0,σ2In)為隨機誤差項向量;wjk為空間權(quán)重矩陣W中的元素，即

其中：sjk為城市j與城市k之間的距離;L為所選城市的數(shù)量。

該模型中未知參數(shù)向量為θ=[β1β2β3β4ρσ2]T。

3.2 參數(shù)估計結(jié)果

利用26個城市9年的數(shù)據(jù),分別應(yīng)用Metropolis within Gibbs算法和ABC算法對模型(9)進(jìn)行參數(shù)估計,并利用估計結(jié)果對第10年的數(shù)據(jù)進(jìn)行擬合分析。假設(shè)π(β)～N(μ,Σ),π(σ2)～I(xiàn)G(α,γ),π(ρ)～U(a,b),其中超參數(shù)的值分別取為:μ=(0.5,0.5,0,0)T,Σ=diag(0.05,0.05,0.1,0.1),α=100,γ=10,a=-1,b=1。在應(yīng)用Metropolis within Gibbs算法的估計過程中,需要通過調(diào)整建議分布的方差(c-2)控制接受概率(確保接受概率在20%～30%)。通過對每M個樣本取平均值降低抽樣樣本的自相關(guān)性,從而使每批次的均值近似獨立,獲得更有效的估計參數(shù),最終需要的樣本容量設(shè)定為N=600(其中前100個樣本作為預(yù)燒期),M=50,c=0.1。

為方便對比,在ABC算法的計算中,樣本容量設(shè)定為N=500;另外,需要通過降低d0來獲得更為精確的估計,同時需要保證計算效率,因此取d0=1.7。

首先選取采用Metropolis within Gibbs算法估計的參數(shù)β1的結(jié)果作為展示,如圖1所示。從圖1a可以看出,截取前N=100個樣本作為預(yù)燒期,剩余樣本構(gòu)成的馬氏鏈已經(jīng)收斂。根據(jù)圖1b,隨著階數(shù)(Lag)值的增大,自相關(guān)函數(shù)(autocorrelation function,ACF)遞減,這說明抽樣容量是合理的,并且該ACF在2倍標(biāo)準(zhǔn)差范圍內(nèi)波動,說明其自相關(guān)系數(shù)截尾,該序列具有平穩(wěn)性特征。其他5個參數(shù)的抽樣序列圖也都能快速達(dá)到收斂趨勢,且相應(yīng)的ACF也隨Lag值增大而減小,由于篇幅有限,在此就不一一闡述。

圖1 參數(shù)β1的MCMC抽樣結(jié)果

通過匯總2類算法的結(jié)果,各參數(shù)的后驗均值、標(biāo)準(zhǔn)誤差、擬合結(jié)果見表1所列。由表1可以看出,根據(jù)Metropolis within Gibbs算法和ABC算法估計參數(shù)得出的估計結(jié)果接近。其中β1和β2后驗均值分別近似為0.6和0.5,這說明服務(wù)業(yè)集聚度和固定資產(chǎn)投資對長三角區(qū)域經(jīng)濟增長效果最顯著。β3和β4均在0附近,說明外商直接投資和勞動力資本投入并不是影響區(qū)域經(jīng)濟的主要因素?？臻g自回歸系數(shù)ρ的估計值約為0.15,說明該長三角地區(qū)存在空間自相關(guān)性,地區(qū)經(jīng)濟增長的強弱會受到相鄰地區(qū)增長的影響。2類算法估計參數(shù)的后驗標(biāo)準(zhǔn)誤差均小于0.005,表明2種算法抽取的樣本對總體來說具有代表性,具有較高的可靠度。在樣本容量相同的情況下,Metropolis within Gibbs算法的四分位距(inter-quartile range,IQR)與標(biāo)準(zhǔn)誤差均小于ABC算法,即利用前者估計參數(shù)比后者具有更高的精度,并且耗時更短。此外,通過這2種算法獲得的估計擬合第10年數(shù)據(jù)的平均誤差約為0.2。

表1 SAR模型估計結(jié)果

綜合來看,Metropolis within Gibbs算法與ABC算法相比,可以提供更為有效和較高精度的估計,這是因為后者依賴于容差閾值和先驗分布的選擇,這兩者對求解精度和效率起到關(guān)鍵作用;與此同時,Metropolis within Gibbs算法對于建議分布的選取、馬氏鏈的收斂以及樣本序列的自相關(guān)性都有較高的要求,一旦要求滿足,后驗概率的估計會相對高效。但Metropolis within Gibbs算法的實現(xiàn)需要似然函數(shù)的表達(dá)式,因此在實際應(yīng)用中,當(dāng)似然表達(dá)式未知時,是難以實現(xiàn)的。此時,在保證符合條件的樣本量足夠多的情況下,ABC算法能夠避開似然的求解,提供另一種有效的方法。

4 結(jié) 論

本文基于長三角26個城市2006-2017年的面板數(shù)據(jù)建立SAR模型,并提出適應(yīng)估計該模型參數(shù)的Metropolis within Gibbs算法和ABC算法進(jìn)行貝葉斯推斷和分析。結(jié)果表明:服務(wù)業(yè)集聚度及固定資產(chǎn)投資的參數(shù)估計對長三角地區(qū)生產(chǎn)總值具有顯著的影響;外商直接投資和勞動力資本投入并不是影響區(qū)域經(jīng)濟的重要因素,說明這兩者并不具備該地區(qū)生產(chǎn)總值的代表性,以上結(jié)論對地區(qū)發(fā)展經(jīng)濟具有參考價值。估計結(jié)果顯示,Metropolis within Gibbs算法雖然能夠有效估計參數(shù),但是依賴于具體的似然函數(shù)的解析式,而ABC算法避免了求解似然函數(shù),能夠有效估計后驗概率分布,2種算法在求解模型中各有優(yōu)勢。本文提出的SAR模型下參數(shù)估計的近似貝葉斯計算方法為空間計量模型的建立提供了一種新思路,可以將它運用于其他帶有隨機誤差項的模型中。