基于ARIMA模型的北京市近十年降水量分析

柳向華++楊碩++劉子康++程一帆++楊毅飛

摘 要：降水量是反映一個地區氣候狀況的重要指標，對林木生長有重要影響，與人類生活息息相關，因此，研究年降水量的演變趨勢和預測年降水量具有重大意義。文章對北京地區2004～2015年的年降水量進行時間序列分析，運用R軟件建立預測模型，并通過模型對北京地區2016～2020年的年降水量進行了預測分析，從而為北京地區的各個領域制定相應政策措施提供參考資料。

關鍵詞：降水量；時間序列；白噪聲；ARIMA模型；確定性分析

一、 背景

水資源是人類賴以生存的資源，降水是自然界水循環的一個環節。如能有效干預和控制降水，對于防洪抗旱，涵養水土，指導生產，均為有利。而干預和控制降水的前提是能對降水量進行有效擬合及預測。目前有多種方法應用于降水量的擬合預測，如運用神經網絡模型，建立微分方程動態模型，為考察不確定性對擬合預測效果的影響，用數值模擬方法建模，建立降雨的隨機模型進行擬合預測，根據擬合預測期長短不同，選用不同的擬合預測方法。其中，美國統計學家G.E.P.BOX和英國統計學家G.M.Jenkins聯合提出的求和自回歸移動平均（autoregressive integrated moving average，ARIMA）模型，作為經典時間序列時域分析方法的核心內容，在降水量及其成分的擬合預測中也有應用。為準確描述降水特征，本文將運用R軟件建立ARIMA月降水量時間序列模型，對北京地區月降水量時間序列數據進行模型擬合。

二、 模型介紹

（一） 模型的結構

具有如下結構的模型稱為求和自回歸移動平均（autoregressive integrated moving average）模型，簡記為ARIMA（p，d，q）模型：

Φ（B）

SymbolQC@ dxt=Θ（B）εt

E（εt）=0，Var（εt）=σ2ε，E（εtεs）=0，s≠t

E（xsεt）=0，s

SymbolQC@ d=（1-B）d；Φ（B）=1-φ1B-…-φpBp為平穩可逆ARMA（p，q）模型的自回歸系數多項式；Θ（B）=1-θ1B-…-θqBq 為平穩可逆ARMA（p，q）模型的移動平滑系數多項式。

上式可以簡記為：

SymbolQC@ dxt=Θ（B）Φ（B）εt

{εt} 為零均值白噪聲序列。

（二） 模型的性質

1. 平穩性

假定{xt}服從ARIMA（p，d，q）模型，記φ（B）=Φ（B）

SymbolQC@ d，φ（B）稱為廣義自回歸系數多項式。

ARIMA模型的平穩性完全由廣義自回歸系數多項式的根的性質決定。可以驗證，自回歸系數多項式的根序列xt的特征根的倒數，當xt隨t趨于無窮趨向于0時序列平穩。即要求xt的特征根都在單位圓內。

容易判斷，ARIMA（p，d，q）模型的廣義自回歸系數多項式共有p+d個根，其中p個根1λ1，…，1λp在單位圓外，d個根在單位圓上。

本例中我們采用圖檢驗的方式判斷序列的平穩性。

2. 方差齊性

對于平穩的時間序列，有方差齊性的性質。故白噪聲序列也是典型的平穩序列，具有方差齊性。本例中主要對殘差序列進行白噪聲檢驗，要求方差齊性。對于殘差異方差序列，我們可以采取方差齊性變換或者ARCH模型進一步提取序列有用信息。

（三） 模型建立與定階

1. 定階依據

根據AR（p）模型的自相關系數拖尾性，偏自相關系數截尾性；MA（q）模型自相關系數截尾性，偏自相關系數圖拖尾性；ARMA（p，q）模型自相關系數和偏自相關系數均拖尾的性質進行模型的定階。

2. 確定適當模型

在分析實際實際的過程中，我們往往不會有原序列就平穩的一些數據，根據Creamer定理，任意一個時間序列{xt}都可以分解成兩部分的疊加，其中一部分是由多項式決定的確定性趨勢成分；另一部分是平穩的零均值誤差成分。

（1） 對于確定性趨勢的提取，在本例中我們采用了兩種方式，用ARIMA模型擬合或者用確定性分析提取。

ARIMA模型中，由于Cramer定理的理論保證，我們對非平穩序列進行差分，d階差分一定可以將充分性信息充分提取。

確定性分析中，直接由序列時序圖判斷用幾次多項式擬合原序列。

（2） 對于零均值誤差部分，它的存在即驗證對于數據中的確定信息我們是否提取充分，要求殘差序列是白噪聲序列。對此，我們分別用了LB檢驗統計量和DW檢驗統計量。

3. 季節模型

在實際數據中，數據往往呈現季節效應。最常用的是加法模型和乘法模型。若差分后序列的季節效應部分振幅根據水平的變化不發生變化，則采用加法模型。否則采用乘法模型擬合季節性趨勢。

（四） 模型優化與比較

通過建立不同的ARIMA模型，擬合原序列。根據原序列或差分后的序列顯示出的特征，不斷進行模型的優化。

根據AIC準則，AIC的值越小，說明模型的擬合效果越好。

（五） 模型預測

xt+l的真實值為：

xt+l=（εt+l+Ψ1εt+l-1+…+Ψl-1εt+1）+（Ψlεt+Ψl+1εt-1+…）

由于εt+l，εt+l-1，…，εt+1的不可獲得性，所以xt+l的估計值為：

x^t（l）=Ψ*0εt+Ψ*1εt-1+Ψ*2εt-2+…

我們考慮真實值與預測值之間的均方誤差：

E[xt+l-x^t（l）]2=（1+Ψ21+…+Ψ2l-1）σ2ε+∑∞j=0（Ψl+j-Ψ2j）2σ2ε

在均方誤差最小原則下，l期預測值為：

x^t（l）=Ψlεt+Ψl+1εt-1+Ψl+2εt-2+…

三、 ARIMA模型

（一） 確定觀察值序列

我們針對北京市04～15年間的年降水量做ARIMA模型的擬合，搜集的數據如下：

（二） 平穩性檢驗

得到需要的數據后，首先將數據變成時間序列形式的變量，并做進一步分析。

初步畫出它的時序圖、自相關系數圖和偏自相關系數圖，檢驗序列的平穩性。

選用R語言作為基本的分析工具

a=scan（'降水量.txt'）

a_ts=ts（a，start=c（2004，1），frequency=12）

plot（a_ts，type='l'）

acf（a_ts）

由時序圖和自相關系數圖判斷，該序列呈現出不平穩性。

1. 尋求平穩序列

由于Cramer分解定理的支持，我們對該序列進行差分運算，可將其確定性信息提取出來。對此我們對序列進行一階差分運算，并檢驗其平穩性。

a_ts_diff=diff（a_ts）

plot（a_ts_diff，type='l'）

acf（a_ts_diff）

通過觀察其一階差分的時序圖及自相關系數圖，原序列的一階差分序列趨于平穩。

同時我們看到原序列的一階差分具有季節效應，即降水量隨月份的變化呈現季節特征。

觀察一階差分后序列的時序圖，年降水量隨著年份的增加振幅不斷增大，我們采用季節乘積模型擬合原序列。

2. 模型定階

作出原序列一階差分后的偏自相關系數圖

pacf（a_ts_diff）

在上述的自相關系數圖和偏自相關系數圖中，不考慮季節因素的影響，我們認為原序列滿足ARIMA（11，1，1）

現在考慮季節效應的影響，對一階差分做12步的差分運算，作出其自相關與偏自相關系數圖

a_ts_dif=diff（a_ts，lag=12）

par（mfrow=c（2，1））

acf（a_ts_dif）

pacf（a_ts_dif）

由一階十二步差分的自相關與偏自相關圖，認為季節影響可用ARIMA（0，1，1）12模型擬合。

3. 模型建立與參數估計

根據選定的模型進行參數估計

estimate=arima（a_ts，order=c（11，1，1），seasonal=list（order=c（0，1，1），period=12））

Call：

arima（x = a_ts，order = c（11，1，1），seasonal = list（order = c（0，1，1），period = 12））

Coefficients：

ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8

-0.0340 0.0099 -0.0534 0.0590 0.0798 -0.0373 0.0379 0.1244

s.e. 0.0882 0.0891 0.0862 0.0855 0.0849 0.0907 0.0891 0.0903

ar9 ar10 ar11 ma1 sma1

0.05000.18100.1961-1.0000-0.8640

s.e.0.08910.08650.08770.07660.1454

sigma^2 estimated as 1054： log likelihood = -653.17， aic = 1334.33

從中我們可以得到原序列的擬合模型的表達式為

xt=1+B

1+0.034B-0.0099B2+0.0534B3-

0.0590B4-0.0798B5+0.0373B6-

0.0379B7-0.1244B8-0.05B9-

0.181B10-0.1961B11（1-0.1454B12）

4. 模型檢驗

對模型的殘差序列進行白噪聲檢驗

tsdiag（estimate）

根據殘差序列的自相關系數圖和p值，p值在0.8到0.1之間，我們不能拒絕殘差序列為白噪聲序列。

5. 模型優化

進一步我們可以進行參數的顯著性檢驗，將對時間序列影響不明顯的參數排除，使得模型更精簡。通過偏自相關系數圖我們可以看出其2，3，4，5，12階系數是不明顯的，由此其實我們可以建立疏系數模型ARIMA（（1，6，7，8，9，10，11），1，1）。在本例中，我們沒有進行下一步的參數顯著性檢驗和疏系數模型。

四、 確定性分析方法與ARIMA模型的結合

（一） 對確定性趨勢進行擬合

通過原序列的時序圖，我們大致得到年降水量隨著年份的增長呈線性趨勢，故用線性方式以自回歸的方式擬合確定性趨勢

len=length（a_ts）

t=c（1：len）

Tt1=lm（a_ts～t）

得到結果：

Call：

lm（formula = a_ts ～ t）

Coefficients：

（Intercept） t

36.2128 0.1019endprint

所以趨勢擬合模型為：T=36.2128+0.1019*t

（二） 殘差自相關檢驗

對原始序列進行初步擬合后，我們需要對其進行殘差的自相關檢驗，以確定是否已經將序列中的確定性信息提取充分。

殘差自相關檢驗：dwtest（Tt1）

Durbin-Watson test

data： Tt1

DW = 1.0205，p-value = 9.613e-10

alternative hypothesis： true autocorrelation is greater than 0

根據p值，我們判斷殘差存在自相關性，說明信息提取不夠充分，需要對殘差進行再次擬合

我們仍采取前面ARIMA模型的方式，對于殘差序列進行擬合

res1=Tt1$res

par（mfrow=c（2，1））

acf（res1）

pacf（res1）

觀察自相關和偏自相關系數圖，偏自相關具有截尾性，并且自相關3階后遞減趨于零。對于殘差序列運用模型AR（3）進行擬合。

rfit1=arima（res1，order=c（2，0，0），include.mean=FALSE）

Call：

arima（x = res1，order = c（2，0，0），include.mean = FALSE）

Coefficients：

ar1 ar2

0.5071 -0.0427

s.e.0.0831 0.0830

sigma^2 estimated as 2308： log likelihood=-762.04， aic = 1530.07

（三） 模型比較

此時我們選用1階延遲序列值為變量再次對殘差序列進行擬合

Tt2=lm（a_ts[-1]～a_ts[-len]-1）

可得模型為xt=0.6875x（t-1）

res2=Tt2$res

acf（res2）

pacf（res2）

rfit2=arima（res2，order=c（6，0，0），include.mean=FALSE）

Call：

arima（x = res2，order = c（6，0，0），include.mean = FALSE）

Coefficients：

ar1 ar2 ar3 ar4 ar5 ar6

-0.0455 0.1148 -0.0339 -0.0214 -0.0059 -0.0965

s.e. 0.0829 0.0828 0.0835 0.0830 0.0826 0.0849

sigma^2 estimated as 2557： log likelihood = -763.99， aic = 1541.97

根據AIC準則，在上述三個模型中，我們認為ARIMA（12，1，1）×ARMA（0，1，1）12模型能更好地擬合北京市近十年的降水量序列。

在用確定性分析方法提取趨勢時，Tt1模型能夠更有效擬合北京市近十年降水量的序列分布

（四） 擬合效果圖

對于上述兩種較好的擬合模型做出其擬合效果圖，并與原序列進行比較

我們更進一步的認定了ARIMA（12，1，1）×ARMA（0，1，1）12對于模型的擬合效果是最好的

（五） 最優模型結構

經過上述步驟，我們認為建立的考慮季節因素的ARIMA模型能夠有效擬合北京市近十年降水量的序列分布

其模型的具體結構為

xt=1+B

1+0.034B-0.0099B2+0.0534B3-

0.0590B4-0.0798B5+0.0373B6-

0.0379B7-0.1244B8-0.05B9-

0.181B10-0.1961B11（1-0.1454B12）

五、 基于ARIMA預測

現在利用我們擬合的ARIMA（12，1，1）×ARMA（0，1，1）12模型對未來五年北京市的降水量做預測

fore=predict（estimate，n.ahead=60）

藍色部分即為對序列的預測。

六、 模型評價

經過上述步驟的建模我們可以觀測到，ARIMA（12，1，1）×ARMA（0，1，1）12對于原始序列的擬合仍然存在一些偏差，對此我們可以利用確定性分析與Auto-Regreesive 模型結合進行擬合。更多的還可以利用Holt 三參數模型對于原始序列建模進行擬合。通過我們所建立的模型進行預測可以得到，如上圖所示，北京市未來五年的降水量在0～150ml左右波動，不會出現極端天氣如暴雨或干旱等災害。

注：本論文只是基于模型進行預測，不保證預測的準確性。

參考文獻：

[1] 中國統計年鑒2005～2017.

[2] 王燕.應用時間序列分析（第四版）.中國人民大學出版社.

作者簡介：

柳向華，楊碩，劉子康，程一帆，楊毅飛，北京市，中國礦業大學（北京）理學院2014級。endprint