王 慧, 吳茜茜

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

0 引 言

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

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

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

1 貝葉斯推斷

1.1 貝葉斯定理

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

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

(1)

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

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

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

1.2.1M-H抽樣

選取建議分布q(·),通?？梢圆捎脤ΨQ式(q(x|y)=q(y|x),q(x|y)=q(|x-y|)),獨(dú)立式(q(x|y)=q(x))或隨機(jī)游走的形式(y=x+η)。具體步驟如下:

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

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

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

(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)置,并且可以計(jì)算參數(shù)θi的后驗(yàn)條件密度π(θ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)其中某個(gè)后驗(yàn)條件分布難以直接采樣時(shí),可以利用M-H抽樣針對該分布進(jìn)行抽樣,由此構(gòu)成Metropolis within Gibbs算法。

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

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

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

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

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

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

對于比較小的d0,最終獲得的模擬樣本分布π(θ|S(Y,Y′)≤d0)可用于近似后驗(yàn)分布π(θ|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表示隨機(jī)誤差項(xiàng),服從MVN(0,σ2In)。模型(2)的待估參數(shù)向量為θ=(β,σ2,ρ)T,且假設(shè)獨(dú)立性存在,似然函數(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、ρ相互獨(dú)立,先驗(yàn)分別為π(β)～N(μ、Σ)、π(σ2)～I(xiàn)G(α,γ)、π(ρ)～U(a,b)[12]。

由獨(dú)立性及似然的表達(dá)可以得出后驗(yàn)的表達(dá)式,可利用M-H抽樣進(jìn)行采樣,但由于待估參數(shù)維度高且參數(shù)類型不同,采用Gibbs抽樣更加合適。根據(jù)Gibbs抽樣,在后驗(yàn)條件分布抽樣過程中,第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的條件后驗(yàn)分別為正態(tài)分布、逆伽瑪分布,可以直接抽樣。但是ρ的條件后驗(yàn)形式比較復(fù)雜,這里選擇應(yīng)用獨(dú)立的M-H抽樣對ρ的條件后驗(yàn)抽樣,由此得到的Metropolis within Gibbs算法的步驟如下:

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

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

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

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

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

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

2.2 ABC方法

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

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

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

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

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

(8)

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

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

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

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

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

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

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

3.1 模型假設(shè)

本文使用長三角26個(gè)城市2006-2017年的面板數(shù)據(jù)進(jìn)行實(shí)際論證。數(shù)據(jù)分別來自《中國城市統(tǒng)計(jì)年鑒》及江蘇、浙江、安徽省統(tǒng)計(jì)年鑒,對于缺失1年的樣本采用均值法補(bǔ)充,對于缺失2年及以上的數(shù)據(jù)采用線性插值法補(bǔ)充。假定社會投入和產(chǎn)出符合Cobb-Douglas生產(chǎn)函數(shù),通過計(jì)算這26個(gè)城市的地區(qū)生產(chǎn)總值(per capita regional gross domestic product, PGDP)的Moran’s I指數(shù)及其統(tǒng)計(jì)量,結(jié)果表明,利用該數(shù)據(jù)建立SAR模型是可行的[18]。通過Pearson相關(guān)性檢驗(yàn)研究變量之間的相關(guān)性,對于存在多重共線性的變量保留其中一個(gè),保證估計(jì)值的唯一性。本文選擇解釋變量為服務(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)為隨機(jī)誤差項(xiàng)向量;wjk為空間權(quán)重矩陣W中的元素，即

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

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

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

利用26個(gè)城市9年的數(shù)據(jù),分別應(yīng)用Metropolis within Gibbs算法和ABC算法對模型(9)進(jìn)行參數(shù)估計(jì),并利用估計(jì)結(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算法的估計(jì)過程中,需要通過調(diào)整建議分布的方差(c-2)控制接受概率(確保接受概率在20%～30%)。通過對每M個(gè)樣本取平均值降低抽樣樣本的自相關(guān)性,從而使每批次的均值近似獨(dú)立,獲得更有效的估計(jì)參數(shù),最終需要的樣本容量設(shè)定為N=600(其中前100個(gè)樣本作為預(yù)燒期),M=50,c=0.1。

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

首先選取采用Metropolis within Gibbs算法估計(jì)的參數(shù)β1的結(jié)果作為展示,如圖1所示。從圖1a可以看出,截取前N=100個(gè)樣本作為預(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個(gè)參數(shù)的抽樣序列圖也都能快速達(dá)到收斂趨勢,且相應(yīng)的ACF也隨Lag值增大而減小,由于篇幅有限,在此就不一一闡述。

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

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

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

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

4 結(jié) 論

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