王浩瑩, 徐賽男, 王亞杰, 馬維軍, 袁超鳳

(黑龍江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 哈爾濱, 150080)

0 引 言

高維矩陣值時間序列數(shù)據(jù)廣泛存在于金融經(jīng)濟(jì)學(xué)、氣象水文和圖像處理等許多領(lǐng)域。這類數(shù)據(jù)復(fù)雜度高，特點難以發(fā)現(xiàn)，因此，對高維矩陣值時間序列數(shù)據(jù)進(jìn)行分析并尋找規(guī)律具有十分重要的意義。

近年來，針對高維矩陣值時間序列數(shù)據(jù)，眾多學(xué)者進(jìn)行了相關(guān)研究[1-4]。2019年，Wang等提出了一個矩陣型因子模型。該模型在保證了數(shù)據(jù)的矩陣結(jié)構(gòu)的同時，實現(xiàn)了顯著的降維[1]。2019年，Chen等基于上述模型通過對行和列施加線性約束，并將先驗知識融入到了模型中，從而增強(qiáng)了觀測矩陣中潛在因子的可解釋性[5]。2020年，Chen等給出了上述矩陣型因子的另一個應(yīng)用研究，即用在動態(tài)國際貿(mào)易數(shù)據(jù)上。由于國際貿(mào)易數(shù)據(jù)可以被整理為一個動態(tài)網(wǎng)絡(luò)數(shù)據(jù)，也可視為矩陣型時序數(shù)據(jù)[6]。2019年，Yu等針對上述矩陣型因子模型給出了一個投影估計方法，并證明了其投影估計有較高的收斂階數(shù)[7]。2020年，Chen等進(jìn)一步考慮了上述矩陣型因子模型的統(tǒng)計推斷問題[8]。

針對上述矩陣型因子模型，給出了模型可識別的充分條件，并從擬似然的角度給出了一個新的參數(shù)估計方法。最后，將上述模型應(yīng)用于我國31個省(包括4個直轄市)之間的人口流動網(wǎng)絡(luò)數(shù)據(jù)，研究其網(wǎng)絡(luò)結(jié)構(gòu)的動態(tài)變化規(guī)律。

1 矩陣型因子模型

1.1 模型介紹

Xt=RFtC′+Et

(1)

vec(Xt)=(C?R)vec(Ft)+vec(Et)

(2)

基于上述模型便可將高維觀測矩陣Xt在時間上的相依性用更低維(k1×k2)的因子矩陣來驅(qū)動，達(dá)到降維的目的。

1.2 可識別性研究

Σ=(C?R)ΣF(C′?R′)+σ2Ip?Iq

(3)

定理1設(shè)θ={R,C,ΣF,σ2}為所有參數(shù)，在模型假設(shè)(IC1)和(IC2)的條件下，模型(1)的協(xié)方差陣結(jié)構(gòu)(2)式滿足可識別，即由Σ(θ)=Σ(θ*)，可得θ=θ*。

1.3 矩陣型因子的擬似然估計

對于傳統(tǒng)的向量型因子模型,2015年，Ng等利用矩陣分解得到了參數(shù)的擬似然估計[11]，2016年，Bai等利用兩步擬似然法求出了估計，并證明了擬似然估計在一定條件下滿足相合性，且具有理想的收斂速度[12]。針對矩陣型因子模型(1)，將通過擬似然的方法(Quasi-maximum likelihood estimation, Q-MLE)，對其參數(shù)進(jìn)行估計。

(4)

則E(MXX)=Σ。因此，構(gòu)建如下的擬似然函數(shù)：

L(θ)=-log|Σ|-tr(MXXΣ-1)

(5)

由Σ的特殊結(jié)構(gòu)(如式(3))，可以求出其逆矩陣及行列式的顯式表達(dá)式(見定理2)，從而可對上述擬似然進(jìn)行化簡。

定理2假設(shè)Σ=(C?R)ΣF(C′?R′)+σ2Ip?Iq，且滿足R′R=pIk1；C′C=qIk2，則其逆和行列式為：

由定理2，擬似然函數(shù)可化簡為：

(6)

接下來便通過以下三步分別給出R,C,{ΣF,σ2}的估計。

(1)給定C、 {ΣF,σ2}及R′R=pIk1, 求R的估計。

通過簡單的矩陣運(yùn)算，式(6)可化簡為：

(7)

表1 更新R為

(2)給定R、{ΣF,σ2}及C′C=qIk2，求C的估計。

方法與求R的估計類似，化簡擬似然函數(shù)式(6)為：

(8)

(3)利用EM算法，估計{ΣF,σ2}。

設(shè)θ1={ΣF,σ2}，給定R和C，可以通過EM過程進(jìn)行估計。令完全的擬似然函數(shù)為：

(9)

進(jìn)一步在正態(tài)假設(shè)下可得潛在因子的條件分布為：

vec(Ft)|vec(Xt)～N(μ,Δ)

(10)

其中μ=ΣF(C′?R′)Σ-1vec(Xt)，Δ=ΣF-ΣF(C′?R′)Σ-1(C?R)ΣF, 則Q函數(shù)為:

(11)

最大化Q函數(shù)有:

(12)

(13)

綜上三步，采用如下的交替最大化算法來估計參數(shù)，如表2所示。

表2 交替最大化算法

2 模擬研究

將通過數(shù)值模擬來驗證上述擬似然估計的有效性，并與文獻(xiàn)[1]給出的估計方法進(jìn)行比較。

2.1 數(shù)據(jù)生成

根據(jù)模型(1)生成觀測數(shù)據(jù)Xt，其步驟如下：

(4) 觀測數(shù)據(jù)Xt通過式(1)生成。

分別考慮了k1=k2=1、k1=k2=2及p,q,T=20, 50, 80, 100的情況，每種情形下重復(fù)500次。

2.2 模擬結(jié)果

表3呈現(xiàn)了k1=k2=1時參數(shù)的估計結(jié)果。其中，σ2=1，Φ=0.55，則因子的方差ΣF=1.433 7。由表3可以看出，隨著p、q和T的增加，兩種參數(shù)的估計效果逐漸變好，且采用的估計方法始終優(yōu)于文獻(xiàn)[1]中的方法。此外，文獻(xiàn)[1]中的方法僅能給出載荷矩陣的估計，而本文的方法不僅能給出載荷的估計，還能給出與方差有關(guān)的參數(shù)的估計。

表3 k1=k2=1的模擬結(jié)果

表4 k1=k2=2的模擬結(jié)果

3 人口流動網(wǎng)絡(luò)數(shù)據(jù)實證分析

人口的遷徙流動一直是涉及地理學(xué)、管理學(xué)和人口學(xué)等學(xué)科的熱門研究課題。目前，人口的流動也是國內(nèi)新冠疫情擴(kuò)散的主要原因[15-17]。本節(jié)將對由百度遷徙大數(shù)據(jù)平臺得到的我國各省之間的人口流動網(wǎng)絡(luò)數(shù)據(jù)進(jìn)行分析，研究其網(wǎng)絡(luò)結(jié)構(gòu)的動態(tài)變化規(guī)律。

3.1 數(shù)據(jù)介紹與整理

本文數(shù)據(jù)來源于百度遷徙大數(shù)據(jù)平臺(https://qianxi.baidu.com/)，該平臺利用百度地圖定位服務(wù)LBS(Location based service)，通過手機(jī)定位信息記錄用戶行為軌跡，記錄出人口流動過程中的起止城市，路徑流強(qiáng)度。研究數(shù)據(jù)包括從2020年10月1日到2021年8月29日(共330日)的逐日的31個省級行政區(qū)(包括直轄市和自治區(qū)，不含港澳臺)之間的省際流動數(shù)據(jù)。初始數(shù)據(jù)包含關(guān)于每個省遷入、遷出人口的規(guī)模指數(shù)，和其遷入、遷出來源省份的比率。時間跨度包含2021年元旦前、春運(yùn)和暑假。由于百度遷徙網(wǎng)站缺失7月8日、8月12日和8月25日的數(shù)據(jù)，故這三日忽略不計。

3.2 結(jié)果與解釋

首先，對于每個窗口，利用文獻(xiàn)[5]中的方法選出因子的個數(shù)，即k1,k2,結(jié)果如表5所示。同時，也列出了當(dāng)k1=k2=3時的累計方差貢獻(xiàn)率。對于每個窗口，當(dāng)因子數(shù)k1=k2=3時，模型至少解釋了總方差的80%以上。故為了便于解釋和比較，對所有時間窗口均采用k1=k2=3來估計載荷矩陣R和C。

表5 因子Ft的維度的估計和k1=k2=3時的累計方差貢獻(xiàn)率

圖1 因子個數(shù)為k1=k2=3時的載荷熱圖

觀察圖1中的(R1)，除1月以廣東為主導(dǎo)外，其他窗口都是長三角(長江三角洲城市群)城市所在省，例如上海、江蘇、浙江、安徽，占據(jù)了第1個遷出樞紐的主導(dǎo)地位。觀察(C1)，除1月以外，長江三角洲城市群也占據(jù)了第1個遷入樞紐的主導(dǎo)地位。這與江浙滬皖三省一市經(jīng)濟(jì)實力最強(qiáng)，城市聚集程度高密切相關(guān)，從而使得人口的短期流動最為活躍。1月對應(yīng)的窗口實際包含12月、1月、2月，包含春運(yùn)期間，是以泛珠三角主要省份(湖南、廣西、江西)為第1個遷入樞紐，以廣東為第1個遷出樞紐。這與他們的地理位置及密切的經(jīng)濟(jì)往來相關(guān)。

圖1中(R2, C2)可分為兩部分，11～12月、4～7月分別是以北京為主導(dǎo)地位的京津冀(北京、天津、河北)地區(qū)為第2個遷出和遷入樞紐。1～3月則是以泛珠三角省份(湖南、廣西、江西)為主的第2個遷出樞紐，同時也是以廣東為主的第2個遷入樞紐。

由圖1中的R3，從11月、12月、1月、2～3月、4～7月數(shù)據(jù)可知，第3個遷出樞紐依次由河北，變?yōu)閺V東，變?yōu)殚L三角、變?yōu)閺V東、變?yōu)楹颖?。由圖1中的C3，在上述時間段中，第3個遷入樞紐依次由河北，變?yōu)榉褐槿恰⒆優(yōu)殚L三角、變?yōu)榉褐槿?，變?yōu)楹颖薄?/p>

綜合圖1來看，山東、河南的人口流動情況較長三角、珠三角和京津冀都較弱，長三角省群樞紐主要分布在第一、三個樞紐，且在所有時間段上一直都是人口流動情況明顯的省群，京津冀省群樞紐在第二、三個樞紐都占了比較明顯的主導(dǎo)位置，主導(dǎo)時間主要是11～12月、4～7月，而珠三角省省群樞紐在一、二、三個樞紐中都占有較大載荷，主導(dǎo)時間主要是12～3月。

圖2 潛在樞紐的人口流動網(wǎng)絡(luò)圖: (a) 2021年2月; (b) 2021年7月

觀察圖2，2月份、7月份的第1個遷出、遷入樞紐均為長三角省群。2月份的第2個遷出、遷入樞紐為泛珠三角省群，7月份以北京為主。2月份的第3個遷出、遷入樞紐為廣東，7月份的為河北。由實線部分可以看出所有樞紐的內(nèi)部人口流動性最強(qiáng)，比如遷出樞紐1與遷入樞紐1之間的實線最粗，說明長三角省群內(nèi)部的人更偏好于在長三角流動，這是受地理因素影響比較大的原因。

最后，為了呈現(xiàn)省份間流動情況的相似性，分別基于2月份和7月份的載荷矩陣進(jìn)行層次聚類(圖3)。本文選取歐幾里得距離和離差平方和法進(jìn)行層次聚類[19-20]。圖3(a)依次2月份R和C的結(jié)果，圖3(b)依次為7月份R和C的結(jié)果。

由圖3可以很容易地識別出長三角群體和珠三角群體。一般來說，地理位置相近、經(jīng)濟(jì)和文化聯(lián)系密切的省份通常會屬于一個群體，人口流動情況會相似。比如2月的廣東、湖南和廣西人口流動很活躍，故整個2月可標(biāo)記為“珠三角活躍”、“長三角-京津冀活躍”和“不那么活躍”。而省份之間的關(guān)系也在發(fā)生變化，2月的珠三角在人口流動網(wǎng)絡(luò)中占主導(dǎo)地位，而后京津冀在人口網(wǎng)絡(luò)中越來越活躍，7月成為了一個單獨的集群。此聚類圖和圖1、圖2所說明的情況也相呼應(yīng)，這也說明了本文的結(jié)果具有一定的可靠性。

圖3 根據(jù)遷出和遷入載荷矩陣對省份聚類: (a) 2021年2月; (b) 2021年7月

4 結(jié) 論

為了更好地理解矩陣值時間序列數(shù)據(jù)的動態(tài)性質(zhì)，本文針對文獻(xiàn)[1]中所研究的矩陣型因子，提出了可識別性條件，并基于擬似然的參數(shù)估計方法，不僅求出了載荷矩陣的估計，還求出了潛在因子矩陣的估計。同時，通過模擬驗證了所提出的估計方法在矩陣維數(shù)較大或者在時間長度較大時也是更具有優(yōu)越性的。用所提出的估計方法針對人口流動網(wǎng)絡(luò)數(shù)據(jù)進(jìn)行了分析，研究其網(wǎng)絡(luò)結(jié)構(gòu)的動態(tài)變化規(guī)律，研究結(jié)果能為中國省際人口流動的動態(tài)性研究提供一種新型的可參照方法。

本文所用的高維矩陣型時間序列模型可以更直觀地呈現(xiàn)出網(wǎng)絡(luò)型數(shù)據(jù)的動態(tài)性變化。由于并沒有對網(wǎng)絡(luò)模型提出任何分布假設(shè)，所得到的載荷矩陣和潛在因子矩陣也都是直接從數(shù)據(jù)網(wǎng)絡(luò)中得到的，而且所提出的方法可以有效地降低動態(tài)網(wǎng)絡(luò)的維數(shù)，并揭示其核心結(jié)構(gòu)，得到的載荷矩陣的估計和潛在因子的估計可以進(jìn)一步用于動態(tài)型網(wǎng)絡(luò)數(shù)據(jù)的研究和預(yù)測，也可以廣泛運(yùn)用到其他動態(tài)因子模型。