廖文輝,林 睿,何志鋒,黃嘉瑩

(1.廣東金融學(xué)院,廣東 廣州 510520;2.江西財(cái)經(jīng)大學(xué),江西 南昌 330032)

0 引言

近年來,顆粒物污染(PM2.5)給城市居民的生活帶來了許多負(fù)面影響,不僅影響了氣象能見度,造成交通阻塞,還會(huì)對(duì)人的健康產(chǎn)生很大影響,例如呼吸系統(tǒng)疾病、心腦血管疾病、精神健康問題、肺癌和過早死亡等都與空氣質(zhì)量有直接或間接的關(guān)系.尤其像京津冀、長(zhǎng)三角、珠三角以及四川盆地等城市群地區(qū),經(jīng)濟(jì)發(fā)達(dá),人口密集,也是空氣污染最為嚴(yán)重的區(qū)域.控制PM2.5濃度因素可以分為污染物排放和氣象要素兩個(gè)方面[1],控制PM2.5濃度的過程包括了排放、傳輸、轉(zhuǎn)換以及沉降[2～5].除排放之外,其余過程都受到天氣條件影響[6],因此構(gòu)建PM2.5濃度統(tǒng)計(jì)模型主要從這兩個(gè)方面入手.Zhu[7]等將滯后一天和滯后兩天的空氣污染物與氣象因素作為PM2.5濃度的影響因素,其中空氣污染物包括PM10、一氧化碳(CO)、二氧化氮(NO2)、臭氧(O3)和二氧化硫(SO2);氣象因素包括氣壓、溫度、相對(duì)濕度、風(fēng)速和日照時(shí)間等.柏玲[8]等研究表明,NO2、CO是影響PM2.5濃度的兩項(xiàng)主要空氣污染物,而降水量和相對(duì)濕度則是影響PM2.5濃度的兩個(gè)重要?dú)庀笠蛩?梅波[9]等研究表明,空氣污染物CO、NO2、SO2對(duì)PM2.5濃度的直接影響為正,對(duì)O3的影響為負(fù),氣象因素中氣溫、氣壓、降雨量以及前一期的降雨均為負(fù)效應(yīng),而濕度為正效應(yīng).

線性回歸是確定兩種或兩種以上變量間相互依賴定量關(guān)系的一種統(tǒng)計(jì)分析方法,應(yīng)用范圍十分廣泛.當(dāng)回歸分析中包括兩個(gè)或兩個(gè)以上的因變量和自變量,且因變量和自變量之間是線性關(guān)系時(shí),即為多元線性回歸模型(Multivariate Linear Model,MLM).在MLM中,最常用的估計(jì)方法是最小二乘法(OLS)估計(jì).現(xiàn)實(shí)中的PM2.5濃度數(shù)據(jù)有很多異常值,比如:監(jiān)控過程中機(jī)器發(fā)生故障、周圍環(huán)境的異常改變等都會(huì)引起異常值的出現(xiàn),而最小二乘法得到的參數(shù)估計(jì)非常容易受到異常值影響,從而影響PM2.5濃度的預(yù)測(cè)準(zhǔn)確性.穩(wěn)健回歸是統(tǒng)計(jì)學(xué)穩(wěn)健估計(jì)中的一種方法,其主要思路是將對(duì)異常值十分敏感的經(jīng)典最小二乘回歸中的目標(biāo)函數(shù)進(jìn)行修改,將穩(wěn)健方法應(yīng)用于回歸模型,以擬合大部分?jǐn)?shù)據(jù)存在的結(jié)構(gòu),同時(shí)可識(shí)別出潛在可能的離群點(diǎn)、強(qiáng)影響點(diǎn)或與模型假設(shè)相偏離的結(jié)構(gòu).當(dāng)誤差服從正態(tài)分布時(shí),其估計(jì)幾乎和最小二乘估計(jì)一樣好,而不滿足最小二乘估計(jì)條件時(shí),其結(jié)果優(yōu)于最小二乘估計(jì).

本文引入常用穩(wěn)健回歸中的S估計(jì)、M估計(jì)、MM估計(jì)、MCD估計(jì)和傳統(tǒng)的OLS估計(jì)建立模型,對(duì)PM2.5濃度進(jìn)行預(yù)測(cè)并比較其效果,為設(shè)計(jì)空氣污染統(tǒng)計(jì)預(yù)測(cè)模型提供選擇和參考.

1 MLM穩(wěn)健估計(jì)

1.1 MLM模型介紹

設(shè)有p個(gè)自變量x1,x2,…,xp,q個(gè)因變量y1,y2,…,yq的n次觀測(cè)樣本(x1t,x2t,…,xpt,y1t,y2t,…,yqt)t=1,2,…,n.矩陣表達(dá)式為:

于是多元線性回歸模型可表示為:

其中B為q×p維系數(shù)矩陣,β0為q維截距向量,e(t)為殘差項(xiàng).

其分量e(t)之間相互獨(dú)立,t=1,2,…,n,每個(gè)e(t)的均值向量為0,協(xié)方差矩陣為∑e.假設(shè)(X,Y)服從聯(lián)合正態(tài)分布,μ表示(X,Y)的聯(lián)合分布的期望向量,∑表示(X,Y)的聯(lián)合分布的協(xié)方差矩陣.顯然,可以利用樣本估計(jì)多元正態(tài)總體參數(shù)μ和∑,再通過分塊矩陣獲得線性模型(3)的參數(shù)β0和B的估計(jì).

1.2 M估計(jì)

M估計(jì)(Maximum Likelihood Type Estimates,M-Estimator)由Huber[10]于1964年提出,也稱作廣義最大似然估計(jì),是最常用的穩(wěn)健估計(jì)方法.許多其它穩(wěn)健估計(jì)方法都是在M估計(jì)的基礎(chǔ)上發(fā)展起來的.

M估計(jì)應(yīng)用在回歸方程中,準(zhǔn)則如下:

其中ri(t)是殘差項(xiàng),ρ(t)為目標(biāo)函數(shù),也稱為殘差損失函數(shù),滿足下述假設(shè)條件:

(ⅰ)ρ(t)是光滑正值實(shí)函數(shù);

(ⅱ)t≥0時(shí),ρ(t)單調(diào)不減;

(ⅲ)ρ(0)=0,maxρ(t)=1.

1.3 S估計(jì)

S估計(jì)追求樣本馬氏距離的某M-尺度統(tǒng)計(jì)量最小,也是單變量尺度參數(shù)M估計(jì)的推廣,故稱為最小化殘差尺度S估計(jì)[11].對(duì)于樣本數(shù)為n的集合,記各樣本馬氏距離為目標(biāo)函數(shù)為ρ(t),也稱為殘差損失函數(shù),這里取光滑且極限為1的正值函數(shù).通過求解得為樣本M-尺度統(tǒng)計(jì)量,這里δ取0.5.顯然樣本馬氏距離與參數(shù)估計(jì)量有關(guān),則多元總體參數(shù)的S估計(jì)定義為

S估計(jì)通常沒有顯性表達(dá)式,一般通過數(shù)值方法近似求解,需要一個(gè)穩(wěn)健且方便計(jì)算但效率一般的估計(jì)量作為迭代的初始值.

1.4 MM估計(jì)

S估計(jì)存在低維情況下估計(jì)效率不高的問題,這主要是因?yàn)榭ǚ椒植荚趐較小時(shí)非對(duì)稱所致.MM估計(jì)則能有效提高低維情況下的估計(jì)效率,其計(jì)算過程如下:首先尋找一個(gè)有偏高穩(wěn)健的初始估計(jì),求初始的樣本馬氏距離與S估計(jì)類似,定義樣本M-尺度統(tǒng)計(jì)量對(duì)于每一個(gè)集合定義估計(jì)量為使得成立的解,那么MM估計(jì)為滿足的解,其中c用于調(diào)整有效性[12,13].

1.5 MCD估計(jì)

MCD估計(jì),也稱最小化協(xié)差陣行列式估計(jì),由Rousseeuw[14]提出.常用樣本協(xié)差陣行列式最小的h個(gè)樣品構(gòu)造的子集估計(jì)均值向量和協(xié)差陣.考慮一個(gè)n行p列的矩陣Xn×p,從中隨機(jī)抽取p+1個(gè)樣本組成子集J,分別計(jì)算T0=ave(J),S0=cov(J),如果S0的行列式 det(S0)=0,則子集J隨機(jī)增加一個(gè)數(shù)據(jù)再計(jì)算S0,直至 det(S0)≠0;計(jì)算n個(gè)樣本數(shù)據(jù)到初始中心T0的馬氏距離選出這n個(gè)距離中最小的h個(gè);然后再通過這h個(gè)樣本計(jì)算樣本均值T1和協(xié)方差矩陣S1,如此迭代稱C-步迭代,可以證明 det(Si)≤det(Si+1),直到當(dāng) det(Sm)≈det(Sm-1)時(shí)停止迭代,這時(shí)再通過Sm進(jìn)行加權(quán)計(jì)算,求出穩(wěn)健的協(xié)方差矩陣估計(jì)量.h值選擇確保穩(wěn)健方法有最高崩潰點(diǎn).Rousseeuw[15]還提出一種快速M(fèi)CD估計(jì)算法,目前在各統(tǒng)計(jì)軟件包中均有使用.

2 實(shí)驗(yàn)分析

根據(jù)世界衛(wèi)生組織數(shù)據(jù)顯示,每年因空氣污染造成的死亡人數(shù)達(dá)700萬.因此,如果能開發(fā)優(yōu)秀的統(tǒng)計(jì)模型對(duì)空氣污染進(jìn)行預(yù)測(cè)掌控,就可以降低或避免其對(duì)居民生活的影響[16].空氣污染主要受氣象條件和污染排放的影響,而這兩方面的實(shí)時(shí)監(jiān)控都存在不少困難,特別是污染排放.本文以廣州萬頃沙環(huán)境監(jiān)測(cè)站2014年9月至2016年12月的空氣污染物濃度數(shù)據(jù)為例(數(shù)據(jù)來源:粵港澳空氣質(zhì)量監(jiān)測(cè)網(wǎng)絡(luò)),引入常用穩(wěn)健方法中的S估計(jì)、M估計(jì)、MM估計(jì)、MCD估計(jì)和傳統(tǒng)的OLS估計(jì)建立回歸模型,對(duì)PM2.5濃度進(jìn)行預(yù)測(cè),比較其效果的差別.

2.1 變量選擇

所用數(shù)據(jù)集的樣本容量n=852;自變量個(gè)數(shù)p=12(包括六個(gè)氣象指標(biāo):風(fēng)速(Spd)、位勢(shì)高度(Hgt)、溫度(Temp)、露點(diǎn)溫度(Dewpt)、海平面氣壓(Slp)、相對(duì)濕度(RH);六個(gè)滯后一天的污染物濃度:二氧化硫(Pre_SO2)、二氧化氮(Pre_NO2)、一氧化碳(Pre_CO)、臭氧(Pre_O3)、顆粒物2.5 mm(Pre_PM2.5)、顆粒物10 mm(Pre_PM10));目標(biāo)變量為當(dāng)前的顆粒物2.5 mm(PM2.5)濃度.目前,氣象預(yù)報(bào)有比較成熟的體系,可以較準(zhǔn)確地預(yù)測(cè)未來三天的氣象條件.因此,基于當(dāng)前氣象條件和歷史污染物濃度作為自變量來構(gòu)建MLM預(yù)測(cè)當(dāng)前空氣污染物濃度是穩(wěn)妥做法.根據(jù)廣義可加模型(GAM)分別對(duì)12個(gè)解釋變量進(jìn)行篩選.首先建立單因素的GAM模型,選擇Identity連接函數(shù),使用樣條平滑函數(shù)進(jìn)行擬合,根據(jù)得出的方差解釋率進(jìn)行初步的變量篩選.結(jié)果見表1.從計(jì)算的結(jié)果分析可得,所有變量的方差解釋度都超過了10%,因此都可以選進(jìn)回歸方程.

表1 PM2.5濃度單因素GAM模型擬合結(jié)果

2.2 穩(wěn)健回歸模型構(gòu)建與預(yù)測(cè)

為了選擇更好的穩(wěn)健估計(jì)方法,進(jìn)行100次隨機(jī)模擬,每次從樣本中抽取30%樣品作為測(cè)試集,其余樣品作為訓(xùn)練集.采用均方預(yù)測(cè)誤差根(mean squared prediction error roots)

自變量包括當(dāng)天的氣象指標(biāo)和前一天的污染指標(biāo)共12個(gè)指標(biāo),因變量為PM2.5.比較方法包括S估計(jì)、M估計(jì)、MM估計(jì)、MCD估計(jì)和傳統(tǒng)的OLS估計(jì).為增加結(jié)果的參考價(jià)值,在樣本中加入一定比例的異常值作比較研究,異常值比例a取0,0.1,0.2和0.4這四種情況,從均值和標(biāo)準(zhǔn)差都為100的正態(tài)分布中取值.比較結(jié)果見表2.

表2 不同污染下各方法的MSPER和Pcor

從表2可以看出,在實(shí)證數(shù)據(jù)集中沒有污染的情況下,S估計(jì)、M估計(jì)、MM估計(jì)和MCD估計(jì)的預(yù)測(cè)結(jié)果都很理想,與傳統(tǒng)最小二乘法較接近.當(dāng)數(shù)據(jù)出現(xiàn)異常值時(shí),傳統(tǒng)方法OLS估計(jì)的預(yù)測(cè)效率馬上下降,測(cè)試集的預(yù)測(cè)誤差明顯變大,預(yù)測(cè)值相關(guān)系數(shù)變小.而穩(wěn)健方法的預(yù)測(cè)效果具有較好的穩(wěn)定性,MCD估計(jì)在各種污染比例下的預(yù)測(cè)效率均最優(yōu),尤其在污染比例小于0.4時(shí),隨著污染比例上升,MCD估計(jì)預(yù)測(cè)效率超過了S估計(jì)、M估計(jì)、MM估計(jì).考慮到污染比例不能超過0.5,上述結(jié)果表明MCD估計(jì)在各種污染比例上均有較好的預(yù)測(cè)效率.

3 結(jié)束語

本文采用多種穩(wěn)健估計(jì)方法,對(duì)空氣污染物濃度進(jìn)行預(yù)測(cè),并比較了這些方法的預(yù)測(cè)效率.對(duì)空氣污染物濃度預(yù)報(bào)的實(shí)證分析結(jié)果表明,在不同污染率的情況下,穩(wěn)健方法的預(yù)測(cè)效果具有較好的穩(wěn)定性,其中MCD估計(jì)在各種污染比例下的預(yù)測(cè)效率均為最優(yōu),能更有效地消除或減弱離群點(diǎn)對(duì)回歸系數(shù)估值的影響,具有更強(qiáng)的穩(wěn)健性,且當(dāng)污染率愈大時(shí),這種效果愈明顯.本文實(shí)證結(jié)果可以為開發(fā)空氣污染統(tǒng)計(jì)預(yù)測(cè)模型提供更多的選擇和思路.