柳向华++杨硕++刘子康++程一帆++杨毅飞

摘 要：降水量是反映一个地区气候状况的重要指标，对林木生长有重要影响，与人类生活息息相关，因此，研究年降水量的演变趋势和预测年降水量具有重大意义。文章对北京地区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