劉玲，宋馬林

(安徽財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽蚌埠　233030)

?

基于ARMA模型的南京市PM2.5濃度分析與預(yù)測(cè)

劉玲，宋馬林*

(安徽財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽蚌埠233030)

[摘要]針對(duì)PM2.5濃度的時(shí)間分布，運(yùn)用南京市2014年P(guān)M2.5日濃度時(shí)間序列，利用時(shí)間序列分析的Box-Jenkins建模思想，結(jié)合EViews7得到符合南京市PM2.5濃度變化的ARMA(5,3)模型，并預(yù)測(cè)出2015年1月1日至5日的PM2.5濃度，通過與實(shí)際數(shù)據(jù)相對(duì)照，發(fā)現(xiàn)模型具有較好的擬合性，適用于大氣中PM2.5濃度分析與預(yù)測(cè).

[關(guān)鍵詞]PM2.5；時(shí)間序列；EViews7；ARMA模型；預(yù)測(cè)

0引言

大氣為地球上的生命繁衍以及人類發(fā)展提供了理想的環(huán)境.它的狀態(tài)和變化，無時(shí)無刻影響著人類的生產(chǎn)、生活甚至生存.近年來，隨著經(jīng)濟(jì)的迅猛發(fā)展，大氣污染問題也越來越嚴(yán)重.空氣質(zhì)量問題已然是政府、環(huán)境保護(hù)部門和全國(guó)人民關(guān)注的熱點(diǎn)問題.據(jù)中國(guó)國(guó)家發(fā)展和改革委員會(huì)2013年7月11日在官方網(wǎng)站上公布的一份報(bào)告披露：自2013年初以來，中國(guó)發(fā)生大范圍持續(xù)霧霾天氣.據(jù)統(tǒng)計(jì)，我國(guó)約1/4的國(guó)土面積、6億人受到霧霾的影響，包括華北平原、黃淮、江淮、江漢、江南、華南北部等地區(qū)[1]，其中，PM2.5被認(rèn)為是造成霧霾天氣的“元兇”.

大氣中PM2.5的污染已經(jīng)引起了國(guó)內(nèi)外環(huán)境和大氣科學(xué)家的廣泛關(guān)注，其中以美國(guó)對(duì)于PM2.5的研究開展得最早，也最為深入.中國(guó)對(duì)于PM2.5等顆粒物的研究起步相對(duì)較晚，近年來由于空氣污染日益嚴(yán)重，國(guó)內(nèi)的研究工作越來越受到重視.前期研究工作主要是分析PM2.5的濃度特征以及成分和來源，楊凌霄(2008)通過對(duì)濟(jì)南市PM2.5的深入研究，闡述了其PM2.5的污染水平和影響因素[2].近年來對(duì)于PM2.5濃度與氣象條件關(guān)系的研究也逐漸增多，劉輝(2011)等通過對(duì)實(shí)測(cè)數(shù)據(jù)的分析，得出氣象條件是除污染源排放外影響PM2.5濃度的重要因素[3].隨著目前學(xué)術(shù)界對(duì)PM2.5研究的逐漸深入，國(guó)內(nèi)相關(guān)研究也逐步增多，近兩年尤以PM2.5的時(shí)空分布居多，張振華(2014)對(duì)北京市PM2.5的污染水平和時(shí)空分布特征進(jìn)行了研究[4]，盧鵬(2014)等運(yùn)用高斯擴(kuò)散模型研究了PM2.5的時(shí)間分布與演變擴(kuò)散[5].

南京市作為長(zhǎng)江三角洲輻射帶動(dòng)中西部地區(qū)發(fā)展的重要門戶，近幾年經(jīng)濟(jì)得到迅速發(fā)展，高能耗企業(yè)頻出，人口密集進(jìn)一步加大，城市建設(shè)與環(huán)境保護(hù)之間的矛盾日益加劇.2013年12月4日南京首次發(fā)布霾紅色預(yù)警，與此同時(shí)，南京市氣象部門統(tǒng)計(jì)數(shù)據(jù)顯示，2013年南京市有242天出現(xiàn)霾，即南京全年有五分之三的時(shí)間處于霧霾之中，這一數(shù)據(jù)達(dá)到了有氣象統(tǒng)計(jì)以來的最高值[6].據(jù)2014年南京市環(huán)境狀況公報(bào)顯示，南京市環(huán)境空氣質(zhì)量超標(biāo)天數(shù)高達(dá)175天[7]，約達(dá)全年天數(shù)的半成，再一次敲響了南京市PM2.5防治的警鐘.為此本文利用南京市2014年P(guān)M2.5日濃度時(shí)間序列，基于ARMA模型對(duì)南京市PM2.5濃度進(jìn)行分析和預(yù)測(cè)，以期為南京市當(dāng)前和未來的PM2.5防治工作提供第一手參考資料.

1ARMA模型的理論介紹及建模步驟

1.1ARMA模型的理論介紹

ARMA(Auto Regressive Moving Average Model)自回歸滑動(dòng)平均模型，是一種隨機(jī)時(shí)序模型，由美國(guó)統(tǒng)計(jì)學(xué)家Box、英國(guó)統(tǒng)計(jì)學(xué)家Jenkins一起創(chuàng)立，因此簡(jiǎn)稱為B-J方法[8].作為一種時(shí)序短期預(yù)測(cè)方法，其具有較高的精度，通常寫為ARMA(p, q)，基本公式如下：

Xt=φ1Xt-1+φ2Xt-2+…+φpXt-p+εt-θ1εt-1-θ2εt-2-…-θqεt-q

(1)

其中，{εt}是白噪聲序列，φ1,φ2,…,φp為自回歸系數(shù)，θ1,θ2,…,θq為移動(dòng)平均系數(shù)，都是模型的待估參數(shù).顯而易見，AR(p)模型和MA(q)模型均為ARMA(p,q)模型的特殊情況.當(dāng)公式(1)中q=0時(shí)，則是自回歸模型AR(p)，當(dāng)公式(1)中p=0時(shí)，則成為移動(dòng)平均模型MA(q).

以上模型針對(duì)的是平穩(wěn)序列，然而非平穩(wěn)序列，需要經(jīng)過差分變換轉(zhuǎn)化為平穩(wěn)序列后才能應(yīng)用于ARMA模型.此時(shí)需要在公式(1)中引入滯后算子B，得到公式(2)：

φ(B)Xt=θ(B)εt

(2)

1.2時(shí)間序列特性的分析工具與方法

(1)序列的平穩(wěn)性檢驗(yàn)

如果時(shí)間序列{Xt},t∈N滿足：

②對(duì)任何t∈N,EXt=μ；對(duì)任何t,s∈N,E[(Xt-μ)(Xs-μ)]=γt-s，就稱{Xt}是平穩(wěn)時(shí)間序列，簡(jiǎn)稱平穩(wěn)序列.

針對(duì)獲取的時(shí)間序列，在建立模型之前，通常需要對(duì)其進(jìn)行平穩(wěn)性檢驗(yàn)，其中ADF檢驗(yàn)法運(yùn)用最為廣泛，下面對(duì)ADF檢驗(yàn)法進(jìn)行簡(jiǎn)要說明：

首先，假定序列{Xt},t∈N服從AR(p)過程.檢驗(yàn)方程為：

▽Xt=γXt-1+ξ1▽Xt-2+…+ξp▽Xt-p+1+εt

其中，εt是白噪聲，若參數(shù)γ<0，則序列是平穩(wěn)的，而當(dāng)參數(shù)γ=0時(shí)，序列至少存在一個(gè)單位根，序列是爆炸性的，沒有實(shí)際意義.因此ADF檢驗(yàn)的原假設(shè)可以表示為：H0:γ=0H1:γ<0.

其次，構(gòu)造ADF檢驗(yàn)統(tǒng)計(jì)量：

(3)

通過公式(2)可以計(jì)算得到一個(gè)t統(tǒng)計(jì)量，然后通過查詢ADF檢驗(yàn)統(tǒng)計(jì)量的臨界值表，決定是否拒絕原假設(shè).如果計(jì)算得到的t統(tǒng)計(jì)量絕對(duì)值超過ADF檢驗(yàn)的臨界值，則拒絕原假設(shè)，認(rèn)為序列是平穩(wěn)的，反之，則認(rèn)為序列是非平穩(wěn)的.

(2)ARMA模型的識(shí)別

ARMA模型的識(shí)別其實(shí)就是對(duì)模型進(jìn)行定階，確定p、q的取值，常用的方法有序列的自相關(guān)圖和偏自相關(guān)圖、AIC準(zhǔn)則.

首先，通過對(duì)平穩(wěn)序列的自相關(guān)圖和偏自相關(guān)圖進(jìn)行分析，有如表1所示的模型識(shí)別依據(jù)[9]：

表1　ARMA模型識(shí)別依據(jù)

當(dāng)運(yùn)用自相關(guān)圖和偏自相關(guān)圖難以確定階數(shù)時(shí)，再采用AIC準(zhǔn)則進(jìn)行判斷.針對(duì)某一p、q取值，當(dāng)AIC(p, q)取最小值時(shí)，p、q為最佳的模型階數(shù)，此時(shí)模型達(dá)到最優(yōu).

1.3ARMA模型的建模步驟

Box-Jenkins提出了針對(duì)時(shí)間序列的建模思想，其建模基本步驟如下：

(1)檢驗(yàn)原始序列進(jìn)行平穩(wěn)性，若序列表現(xiàn)為不平穩(wěn)，可對(duì)序列進(jìn)行差分變換使其達(dá)到平穩(wěn)；

(2)求出該觀察值序列的樣本自相關(guān)函數(shù)和樣本偏自相關(guān)函數(shù)，結(jié)合AIC準(zhǔn)則，選擇適當(dāng)p、q擬合ARMA(p, q)模型；

(3)運(yùn)用OLS估計(jì)ARMA(p, q)模型中的未知參數(shù)；

(4)檢驗(yàn)?zāi)Ｐ偷暮侠硇裕?dāng)擬合模型無法通過檢驗(yàn)時(shí)，需返回步驟(2)，重新確定p、q；

(5)模型優(yōu)化，即時(shí)擬合的模型通過檢驗(yàn)，仍然可以返回步驟(2)重新確定p、q來模型進(jìn)一步優(yōu)化，從而選取最優(yōu)模型；

(6)根據(jù)最終確定的擬合模型，預(yù)測(cè)序列的未來趨勢(shì).

2南京市M2.5濃度的實(shí)證分析與預(yù)測(cè)

通過中國(guó)空氣質(zhì)量在線監(jiān)測(cè)平臺(tái)選取南京市2014年1月1日至12月31日PM2.5濃度數(shù)據(jù)[9]，共365個(gè)樣本.下面旨在利用ARMA模型的建模理論結(jié)合EViews7進(jìn)行南京市PM2.5濃度的實(shí)證分析與預(yù)測(cè).

2.1原始數(shù)據(jù)平穩(wěn)性檢驗(yàn)

首先，運(yùn)用EViews7畫出南京市2014年1月1日至12月31日PM2.5濃度時(shí)間序列圖，如圖1所示，通過對(duì)時(shí)序圖的分析，可以大致認(rèn)為該序列沒有明顯的趨勢(shì)，初步判斷序列是平穩(wěn)的.

圖1　PM2.5濃度時(shí)間序列圖

接下來，采用ADF檢驗(yàn)對(duì)原序列進(jìn)行平穩(wěn)性檢驗(yàn)，若數(shù)據(jù)沒有通過檢驗(yàn)，則說明原始序列不平穩(wěn).運(yùn)用EViews7繼續(xù)ADF檢驗(yàn)，結(jié)果如表2所示.

表2　序列PM2.5的ADF檢驗(yàn)結(jié)果

ADF檢驗(yàn)的原假設(shè)為序列PM2.5存在一個(gè)單位根，由表2知，ADF檢驗(yàn)的t統(tǒng)計(jì)量值為-8.32，比1%顯著性水平臨界值-3.45要小，因此在99%的置信水平下，可以拒絕原假設(shè)，認(rèn)為序列PM2.5不存在單位根，即序列平穩(wěn).

圖2　PM2.5序列自相關(guān)—偏自相關(guān)分析圖

2.2模型識(shí)別

運(yùn)用EViews7畫出原始序列的自相關(guān)圖和偏自相關(guān)圖，如圖2所示，接下來，通過分析圖2進(jìn)行模型的識(shí)別和定階.

總觀圖2，可以看出，原始序列的自相關(guān)圖和偏自相關(guān)圖都沒有明顯的截尾性，因此需要嘗試使用ARMA模型.首先，通過觀察圖2左邊的序列自相關(guān)圖，可以看出序列明顯地短期相關(guān)，延遲1階、2階、3階以及5階的相關(guān)系數(shù)均顯著不為0，延遲4階的相關(guān)系數(shù)近似為0；且序列相關(guān)系數(shù)由非零衰減為小值波動(dòng)的過程較為連續(xù)和緩慢，因此可以判定該序列的自相關(guān)系數(shù)具有拖尾性，可以考慮取q=2或q=3.再觀察圖2右側(cè)的序列偏自相關(guān)圖，可以看出偏自相關(guān)系數(shù)除了在延遲1階、2階以及5階時(shí)顯著大于2倍標(biāo)準(zhǔn)差，其他延遲階數(shù)時(shí)，基本在2倍標(biāo)準(zhǔn)差范圍內(nèi)，因此，可以考慮P取5.

綜上，本文可以對(duì)PM2.5序列建立ARMA(5,2)模型或者ARMA(5,3)模型.

2.3模型的參數(shù)估計(jì)

確定模型的階數(shù)后，需要對(duì)模型進(jìn)行參數(shù)估計(jì).

表3ARMA(5,2)模型參數(shù)估計(jì)及檢驗(yàn)結(jié)果

圖3　單位根檢驗(yàn)

運(yùn)用EViews采用非線性方法分別對(duì)ARMA(5,2)和ARMA(5,3)模型進(jìn)行參數(shù)估計(jì)與檢驗(yàn)，結(jié)果分別如表3和表4所示，兩表上部分與普通最小二乘估計(jì)結(jié)果一樣.此時(shí)需要根據(jù)模型調(diào)整后的可決系數(shù)以及AIC和SC準(zhǔn)則等判斷模型的整體擬合效果.兩表中最下面三行均為滯后多項(xiàng)式φ(x-1)=0和θ(x-1)=0的倒數(shù)根，當(dāng)?shù)箶?shù)根都位于單位圓之內(nèi)時(shí)，即可判定過程是平穩(wěn)的.由圖3可知，表3和表4中的根都在單位圓內(nèi)，均符合要求.

表4ARMA(5,3)模型參數(shù)估計(jì)及檢驗(yàn)結(jié)果

由表4可知，ARMA(5,3)模型調(diào)整后的R2為0.471741大于表3中ARMA(5,2)模型的0.43824，而AIC和SC值分別為9.564024和9.6399，分別小于表3中的9.622771和9.687808，因此認(rèn)為ARMA(5,3)模型更合適，其對(duì)應(yīng)的模型表達(dá)式為：

Xt=1.28Xt-1-0.36Xt-2+0.07X5+εt+0.52εt-1+0.37εt-2

2.4模型的檢驗(yàn)

參數(shù)估計(jì)后，還需要進(jìn)一步檢驗(yàn)殘差序列et是否為白噪聲，即當(dāng)滯后期k≥1時(shí)，et的樣本自相關(guān)系數(shù)是否大致為0.當(dāng)殘差序列不是白噪聲時(shí)，表明殘差序列中還有少量有用信息沒有被提取出來，此時(shí)需要對(duì)模型進(jìn)行改進(jìn).模型的檢驗(yàn)一般側(cè)重于檢驗(yàn)殘差序列的隨機(jī)性，通過對(duì)殘差序列進(jìn)行χ2檢驗(yàn)，可以大致判斷其是否是純隨機(jī)序列.首先，給出殘差序列的自相關(guān)函數(shù)：

則接受原假設(shè)，認(rèn)為模型通過檢驗(yàn)，即殘差序列{et}之間相互獨(dú)立，；否在檢驗(yàn)不通過[8].

圖4　ARMA(5,3)殘差序列的自相關(guān)—偏自相關(guān)分析圖

2.5模型的預(yù)測(cè)

接下來運(yùn)用已經(jīng)建立好的ARMA(5,3)模型，通過EViews7對(duì)南京市2014年1月1日至12月31日PM2.5濃度進(jìn)行擬合，得到圖5.

圖5　南京市2014年1月1日至12月31日PM2.5濃度擬合預(yù)測(cè)效果

如圖5，紅色、綠色、藍(lán)色曲線分別表示序列觀察值、擬合序列值以及序列的殘差.可以看出，擬合值與觀測(cè)值十分接近，模型的擬合預(yù)測(cè)效果良好.

下面采用線性最小方差，利用ARMA(5,3)模型對(duì)PM2.5序列的未來發(fā)展進(jìn)行預(yù)測(cè)，給出2015年1月1日至5日五天的預(yù)測(cè)值，如表5所示.

表5　2015年1月1日至5日南京市PM2.5預(yù)測(cè)值與真實(shí)值比較

從表中可看出，置信區(qū)間長(zhǎng)度遠(yuǎn)遠(yuǎn)大于預(yù)測(cè)值與真實(shí)值的誤差范圍，所以本文建立模型的預(yù)測(cè)效果較好.由此也進(jìn)一步驗(yàn)證了本文所建立模型的準(zhǔn)確性，對(duì)目前空氣中PM2.5預(yù)測(cè)具有較好的實(shí)際意義.

3結(jié)論

本文利用時(shí)間序列分析的Box-Jenkins建模思想，對(duì)南京市2014年1月1日至12月31日PM2.5濃度這一時(shí)間序列進(jìn)行模型的建立和實(shí)證分析，得到符合南京市PM2.5濃度變化的ARMA(5,3)模型，了解了南京市PM2.5變化的基本特征.

首先，針對(duì)原始數(shù)列進(jìn)行平穩(wěn)性檢驗(yàn)得到原始數(shù)據(jù)基本平穩(wěn)；其次，運(yùn)用原始序列的自相關(guān)圖和偏自相關(guān)圖進(jìn)行模型識(shí)別，確定可建立的模型為ARMA(5,2)和ARMA(5,3)；然后，對(duì)兩個(gè)模型進(jìn)行參數(shù)估計(jì)與檢驗(yàn)，確定合適的模型為ARMA(5,3)，并根據(jù)參數(shù)估計(jì)結(jié)果給出模型的表達(dá)式；隨后，對(duì)已識(shí)別的模型進(jìn)行檢驗(yàn)，充分驗(yàn)證了所設(shè)計(jì)模型的合理性；最后，運(yùn)用已經(jīng)確立好的模型，對(duì)南京市2014年1月1日至12月31日PM2.5濃度進(jìn)行擬合，通過與實(shí)際值進(jìn)行比較，再一次驗(yàn)證了模型的準(zhǔn)確性，并對(duì)南京市2015年1月1日至5日的PM2.5濃度做了短期預(yù)測(cè)，發(fā)現(xiàn)模型具有較好的預(yù)測(cè)效果.

綜上所述，ARMA模型較好地解決了大氣中PM2.5的時(shí)間分布問題，借助EViews軟件能夠十分便捷地將ARMA模型應(yīng)用于大氣中PM2.5濃度的分析與預(yù)測(cè)，為國(guó)家和地方政府進(jìn)行空氣質(zhì)量監(jiān)測(cè)、預(yù)報(bào)、控制以及制定相應(yīng)政策、法規(guī)和管理辦法提供決策指導(dǎo).

參考文獻(xiàn)

[1]國(guó)家發(fā)展改革委環(huán)資司.節(jié)能減排形勢(shì)嚴(yán)峻產(chǎn)業(yè)發(fā)展?jié)摿薮蟆?013年上半年節(jié)能減排形勢(shì)分析[J].中國(guó)經(jīng)貿(mào)導(dǎo)刊.2013:14-15.

[2]楊凌宵.濟(jì)南市大氣PM2.5污染特征、來源解析及其對(duì)能見度的影響[D].山東大學(xué),2008.

[3]劉輝,賀克斌,巧水亮,等.2008年奧運(yùn)前后北京城、郊PM2.5及其水溶性離了變化特征[J].環(huán)境科學(xué)學(xué)報(bào).2011,31(1):177-185.

[4]張振華.PM2.5濃度時(shí)空變化特性、影響因素及來源解析研究[D].浙江大學(xué),2014.

[5]盧鵬,何杰.PM2.5的時(shí)間分布與演變擴(kuò)散研究[J].西南民族大學(xué)學(xué)報(bào)2014,40(1):66:70.

[6]中國(guó)新聞網(wǎng).去年南京霧霾242天“防霾”成南京兩會(huì)熱議詞[2015-11-2].http://www.chinanews.com/sh/2014/01-12/5725911.shtml.

[7]南方日?qǐng)?bào).2014年南京市環(huán)境狀況公報(bào).[2015-11-2].http://njrb.njdaily.cn/njrb/html/2015-06/05/content_158622.htm.

[8]易丹輝.數(shù)據(jù)分析與EViews應(yīng)用[M].北京:中國(guó)統(tǒng)計(jì)出版社,2002.

[9]馮盼,曹顯兵.基于ARMA模型的股價(jià)分析與預(yù)測(cè)的實(shí)證研究[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)2011,41(22):85-89

[10]中國(guó)空氣質(zhì)量在線監(jiān)測(cè)平臺(tái)[2015-10-31]. http://www.aqistudy.cn/.

[責(zé)任編輯：房永磊]

PM2.5 Concentration Analysis and Prediction in Nanjing Based on ARMA Model

LIU Ling, SONG Ma-lin*

(Anhui University of Finance and Economics Institute of Statistics and Applied Mathematics, Bengbu 233030,China)

Abstract:In connection with PM2.5 concentration in different time, use PM2.5 concentration of everyday in Nanjing in 2014, apply Box-Jenkins time series analysis modeling ideas, combined EViews7, get ARMA (5,3) model conform with the changes ofPM2.5 concentration in Nanjing and predict the concentration of PM2.5 from January 1, 2015 to 5th, by contrast with the actual data and found that the model has better fitting, is suitable for the analysis and prediction of PM2.5 concentration in the atmosphere.

Key words:PM2.5; time series; eviews7; ARMA model; prediction

[中圖分類號(hào)]X511

[文獻(xiàn)標(biāo)識(shí)碼]A

[文章編號(hào)]1004-7077(2016)02-0054-09

[作者簡(jiǎn)介]宋馬林(1972-)，男，安徽蚌埠人，安徽財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院教授，碩士生導(dǎo)師，主要從事資源環(huán)境統(tǒng)計(jì)、數(shù)量經(jīng)濟(jì)的研究.

[基金項(xiàng)目]安徽省級(jí)創(chuàng)新創(chuàng)業(yè)項(xiàng)目(項(xiàng)目編號(hào)：AH201410378516).

[收稿日期]2016-02-01