北京市月平均气温分析以及未来温度预测

北京市月平均气温分析以及未来温度预测

北京交通大学李荣荣、莫海燕、田龙

摘要：本世纪以来所进行的一些科学观测表明，大气中各种温室气体的浓度都在增加。温室气体吸收长波辐射并再反射回地球，从而减少向外层空间的能量净排放，大气层和地球表面将变得热起来，这就是“温室效应”。在过去一个世纪里，全球表面平均温度已经上升了0．3℃到0．6℃，全球海平面上升了10到25厘米。气候变暖是全世界都关注的问题。我们以统计的方法分析北京近40几年来月平均气温的数据，对1951.1-2009.12建立时间序列模型，并结合数据的季节性、趋势性等特点，通过建立乘积性季节模型，对北京市进几十年的温度变化趋势作出统计性的解释，并对北京未来的气温作出预测，与余留12个月的月平均气温相比，进行误差分析。

关键词：时间序列；乘积型季节性模型；趋势分析；气温预测

The Analysis of Monthly Mean Temperature and Forecast of the

Future Temperature of Beijing

Abstract: In this century a number of scientific observations have shown that atmospheric concentrations of greenhouse gases are increasing. Greenhouse gases absorb long wave radiation, and then reflected back to Earth, thereby reducing net emissions of energy to outer space, the atmosphere and the Earth's surface will become hotter, and this is the "greenhouse effect." In the past century, global average surface temperature has increased by 0.3 ℃to 0.6 ℃, the global sea level has risen 10-25 cm. Climate warming is concerned by the world. It is analysis of two scores of years of monthly mean temperature data of Beijing by statistical method. According

to 1951.1-2009.12, build time series models, combined with seasonal data, trends and other characteristics, and creating Seasonal model of product type of a few decades of temperature of Beijing ,get the interpretation of statistical trends, and prediction of future temperature trends in Beijing , with the remaining 12 months compared to the average monthly temperature, finally make error analysis.

Key words: Time Series，Seasonal model of product type，Analysis of trend，Temperature Forecast

目录

目录 (2)

一、问题描述 (3)

二、数据来源以及解释 (3)

2.1 数据来源 (3)

2.2 数据说明 (3)

三、数据分析与模型建立 (3)

3.1 数据初处理 (3)

3.2 数据的季节性以及趋势性 (4)

1.季节性 (4)

2.趋势性 (4)

3.3 平稳化及其检验 (5)

3.3.1 平稳化检验方法 (5)

3.3.2 自相关图检验 (6)

3.4 时间序列模型 (9)

3.4.1 季节调整的概念和作用 (9)

3.4.2 季节调整的模型和方法 (10)

3.4.3 模型的数学理论 (10)

3.4.4 季节性模型 (13)

3.5 乘积模型 (13)

3.5.1 乘积模型构造原理 (13)

3.5.2 乘积模型的定阶和模型拟合 (15)

3.5.3 模型预测 (17)

3.5.4 误差分析 (18)

3.5.5 预测模型的评价 (18)

四、模型讨论 (19)

五、参考文献 (19)

六、附件 (19)

一、问题描述

“全球气候变暖”一直是一个科学界主要课题，然而对于局部环境而言，极端气候的一再出现，很多人去气候变暖提出诸多质疑，以北京为例，2009年得北京就经历的罕见的寒冬，11月份大学飞扬，究竟气候是否在变暖？我们以统计的方法分析北京近40几年来月平均气温的数据，对北京市进几十年的温度变化趋势作出统计性的解释，并对北京未来的气温作出预测。

二、数据来源以及解释

2.1 数据来源

本文选用的数据样本为1951年1月到2009年12月的北京月平均气温值，数据来自中国气象科学数据共享服务网(https://www.360docs.net/doc/3f16705398.html,)，中国地面国际交换站气候资料月值数据集，我们这里仅选取地点为北京（54511）。

2.2 数据说明

从数据源中提取1951年1月到2010年12月北京市月平均气温数据，删去其他无效数据，月平均气温数据全部为整形实数，单位为0.1℃，共720条记录。

三、数据分析与模型建立

3.1 数据初处理

记该数据样本为。

记1951年1月为起始月份1，其他依次累计得到即时间序列

{},1,2,720

X n= 。

n

这样得到分析的数据对象。

3.2 数据的季节性以及趋势性

1.季节性

以累计时间序列为x轴，月平均气温为y轴，做出时间序列图：

从时间序列图中可以明显看出数据的季节性，而且数据季节性的特征是显而易见的，气温有一定的增长性趋势，但是气温增长的趋势性并不明显。需对数据做出进一步处理。

2.趋势性

对数据进行行平稳性daniel检验：

Daniel检验方法是建立在Spearman相关系数的基础上，对时间序列的样本

记的秩为,考虑到变量对,t=1,2,……,n的

Spearman秩相关系数，有

做下列检验：

序列：序列

Daniel检验方法：对于显著性水平，由时间序列计算

，t=1,2,……,n的Spearman秩相关系数,若

则拒绝，认为序列非平稳。当时，认为序列有上升趋势，时，认为序列有下降趋势。又当时，接受，可以认为是平稳序列。

对于源数据在显著水平，经由SAS软件得到

由于

故拒绝，认为序列是非平稳的，又>0,故说明近30年北京月平均气温整体有上升趋势。

由以上分析可以得知数据有明显的趋势性和季节性。

3.3 平稳化及其检验

3.3.1 平稳化检验方法

为了使该时间序列具有可预测性，首先我们要消除其季节性和趋势性，使其

变得平稳化并对其检验平稳性，这样才可以进一步的做预测分析。

检验平稳性的一般方法主要有两种：

1、根据时序图和自相关图显示的特征作出的图检法；

2、构造检验统计量进行假设检验。

我们采取利用自相关图检验构造数据的平稳性。

3.3.2 自相关图检验

由于数据有明显的季节性和趋势性，显然不可能得到平稳的自相关图，实际上，其自相关图如下：

显然自相关系数明显大于2倍的标准差，该时间序列不是平稳的序列。由于该时间序列很显然的以12为周期，因此为了消除其季节性，令

得到新的时间序列，其时间序列图如图三所示，

又由于预测的数据有趋势性，因此为了消除趋势性，再做一次差分，即令

得到序列,我们称这样的处理为对原数据的“一阶12步差分”，其时间序列图如图四所示：

显然，由时序图可以直接看出图四的平稳性明显好于图三，因此我们选择用序列，下面得到序列的自相关图，进行进一步的平稳性分析：

显然经过处理后的数据的平稳性明显好于原始数据的平稳性，其相关性系数基本上控制在2倍的标准差范围之内。延迟一步的自相关系数大于2倍的标准差这说明差分后序列具有短期的相关性。还有其延迟12步的自相关系数明显大于2倍的标准差范围，说明差分后的序列仍蕴含非常显著的季节效应。

再做其偏自相关图，如图六：

结论与其自相关图的结论一致，但与自相关图不同的是，偏自相关图显示出非截尾的性质，综合该序列的字相关和偏自相关系数的性质，为拟合模型定阶为MA（1）。

在这种情况下尝试ARMA模型拟合，反复尝试实验结果均不理想，说明简明的ARMA模型并不适合该序列，考虑到该序列既有短期相关性又具有季节效应，短期相关性和季节效应不能简单地、可见性地提取，因此估计该序列的季节效应和短期相关性之间有更为复杂的关联性。

此时，则假定短期相关性和季节效应之间具有乘积关系，尝试使用乘积模型来拟合该序列。

3.4 时间序列模型

3.4.1 季节调整的概念和作用

以月份或季度作为时间观测单位的时间序列通常具有一年一度的周期性变化，这种周期变化是由于季节因素（气候）的影响造成的，称为季节性波动。月度和季度的时间序列的季节性波动是非常显著的，因为季节因素的存在，同一年

中不同月份或季度的数据往往不具有可比性，我们传统上的做法通常是用同比来其增长变化，但它不能及时反映当前变化的走势。因此，在使用月度或季度数据进行分析之前，需要对数据进行“季节调整”，季节调整后的数据消除了季节性的影响，使得不同月份或季度之间的数据具有可比性，可以更及时的反映气温的变化。

所谓“季节调整”，就是一个从时间序列中估计和剔除季节因素影响的过程，目的是反映序列真正的客观规律和趋势。除了季节因素外，一个时间序列通常还受多种因素影响，一般地，我们把所有这些因素分解为趋势因素（T）、循环因素（C）、季节因素（S）和不规则因素（I）。其中，趋势因素反映了气温的长期演变方向，是上升、持平还是下降；循环因素（周期因素）反映了时间序列持续性的周期波动，侧重时间序列是处于周期的上升阶段、下降阶段还是转折阶段，实际工作中趋势与循环因素往往放在一起分析不进行区分；季节因素反映时间序列在不同年份的相同季节（同一季度，同一月份）所呈现出的周期性变化；不规则因素反映的是前三个因素无法解释的误差或随机因素产生的变化。季节调整后的时间序列就是趋势、循环和不规则因素的合成。

3.4.2 季节调整的模型和方法

根据时间序列（Y）的四个构成要素的关系，对时间序列作季节调整的分解模型主要有两种形式：加法模型和乘法模型。

加法模型的表达式是：Y=T+C+S+I。它的特点是，原始序列由四个要素相加构成，各要素都用绝对量表示，与原始序列的计量单位相同，直观性比较好；缺点是不同的变量之间缺乏可比性。加法模型适用于T、S、C相互独立的情形。

乘法模型的表达式是：Y=T×C×S×I。它的特点是，原始序列由四个要素相乘构成，趋势要素是绝对量，其他三个要素是相对量，加强了不同变量间的可比性。乘法模型适合于T、S、C相关的情形。由于时间序列分解的四大要素一般都存在相互影响，因此一般我们都采用乘法模型进行季节调整。

3.4.3 模型的数学理论

1、AR(p)序列

设{t X ，t =0,±1, ±2，…}是零均值平稳序列，满足下列模型：

t X =11-t X ?+22-t X ?+…+p t p X -?+t ε， （3.1）

其中t ε是零均值、方差是2εσ的平稳白噪声，则称t X 是阶数为p 的自回归序列，简记为AR (p )序列，而

?T

),,,(p 21??? =

称为自回归参数向量，其分量j ?，j =1,2，…，p 称为自回归系数。

引进推移算子对描述式（3.1）比较方便，算子B 定义如下： t BX ≡1-t X ,t k X B ≡k t X - 记算子多项式)(B ?=1-B 1?-B 2?-…-p p B ?， 则式（3.1）可以改写为

)(B ?t X =t

ε

2、 MA (q )序列

设{t X ，t =0,±1, ±2，…}是零均值平稳序列，满足下列模型：

t X =t ε-11-t εθ-22-t εθ-…-q t q -εθ， （3.2）

其中t ε是零均值、方差是2εσ的平稳白噪声，则称t X 是阶数为q 的滑动平均序列，简记为MA(q )序列，而 θ=T q ),,,(21θθθ

称为滑动平均参数向量，其分量j θ，j =1,2，…，q 称为滑动平均系数。

对于线性推移算子B ，有

t B ε≡1

-t ε，

t k

B ε≡k

t -ε；

再引进算子多项式

)(B θ=1-B 1θ-22B θ-…-q q B θ， 则式（3.2）可以改写为

t X =)(B θt ε

3、 ARMA （q p ,）序列

设{t X ，t =0,±1, ±2，…}是零均值平稳序列，满足下列模型：

t X -11-t X ?-22-t X ?-…-p t p X -?=t ε-11-t εθ-22-t εθ-…-q t q -εθ， （3.3）

其中t ε是零均值、方差是2εσ的平稳白噪声，则称t X 是阶数为p ，q 的自回归滑动平均序列，简记为ARMA(q p ,)序列。当q =0时，它为AR （p ）序列；当p =0时，它为MA(q )序列。

应用算子多项式)(B ?，)(B θ，式（3.3）可以写为

B ?t X =)(B θt

ε 4 、ARIMA （q d p ,,）模型

ARIMA 模型即为自回归求和滑动平均模型，它的d 阶差分就是属于平稳性质的ARMA （自回归滑动平均混合模型）模型。

设{t X ，t =0,±1, ±2，…}为非平稳序列，但存在正整数d ，使得d ?t X =t

W 为平稳序列，而{t W ，t =0,±1, ±2，…}可以由ARMA （q p ,）模型拟合，则称t X 是ARIMA （q d p ,,）序列。这时，t X 满足：

()B ?d

?

t x = ()B θt ε （3.4）

p

，q 分别为d ?t X =t W 的自回归阶数和滑动平均阶数，d 为差分阶数。

()B ?=1-B 1?-22B ?-…-p p B ?， )(B θ=1-B 1θ-22B θ-…-q q B θ。

3.4.4 季节性模型

在不少实际问题中，时间序列有明显的周期性规律。由季节性因素或其他周期因素引起的周期性变化的时间序列，我们称为季节性时间序列，相应的模型为季节性模型。

对于月度数据，为了消除周期性变化因素，可以进行差分运算：

12?t X =t X -12-t X =)1(s

B -t X

一般地，对周期s 的序列，可进行差分经运算：

s ?t X =)1(s

B -t X

D

s ?t

X

=D s B )1(-t X

对于季节性模型，可以用一类疏系数ARIMA 模型去描述。设s 为一正整数，一个时间序列{t X ，t =0,±1, ±2，…}如果满足下列模型：

)(s B ΦD

s ?t X =)(s

B Θt E ， （3.5）

则称

t

X

是周期为

s

的季节性序列，其中，

?????Θ--Θ-Θ-=ΘΦ--Φ-Φ-=ΦQs Q s s s Ps

P s s s B

B B B B B B B 2212211)(1)(

3.5 乘积模型

3.5.1 乘积模型构造原理

当序列具有短期相关性时，通常可以使用低阶

模型提取。当序列

具有季节效应，季节效应本身还具有相关性时，季节相关性可以使用以周期步长为单位的

模型提取。有的季节性时间序列不仅含有季节性成分,还混

有非季节性成分,若单一用季节性或非季节性ARIMA 模型进行分析,其预测效果往往不理想。季节性ARIMA 与非季节性ARIMA 的混合效应通常表现为相乘的形式。

式（3.5）中t E 一般不必是白噪声，而可设它是另一个ARIMA （q d p ,,）序列。

()B ?d

?t E = ()B θt ε

由式（3.5），

)(s B Φd ?D

s ?t X =)(s B Θd

?

t E

()B ?)(s B Φd ?D

S

?t X =)(s B Θ()B ?d ?t E =)(s B Θ()B θt ε 令t W =d ?D s ?t X 得

()B ?)(s B Φt W =()B θ)(s

B Θt ε （3.6）

式（3.6）称为乘积季节性模型，其阶数常用s Q D P q d p ),,(),,(?表示，式（3.6）中的t W 是由原来序列t X 经差分d ?D s ?得到的，经过这样的运算，可以消除序列

t X 的非平稳性（即趋势性）与季节性（周期性）。式（3.6）是一种疏系数模型，

它反映了序列相邻时刻与相隔为周期s 的时刻之间复杂变化的规律。

式（3.6）又可以记做：

d

?D

S

?t X =()()

()()

s t

s

B B

B B θε

?ΘΦ

其中，

我们简记该模型为

。ARIMA 建模法分为三个阶

段：模型识别，通过不同的p ，q 的参数，比较不同的乘积模型，选取最合适的p ，q 参数，拟合出合适的序列。参数估计和诊断检验，通过这三个处理步骤的反复进行,最终确定一个用于预报或控制的最优模型。这三个阶段都需依赖计算机软件进行大量的计算。

3.5.2 乘积模型的定阶和模型拟合

对于ARMA的定阶可以根据自相关和偏相关图的结尾和拖尾性来定，然后再通过统计量检验拟合的效果。但是一般数据的自相关和偏相关图不会像特殊情况那么有规律，那就根据自相关和偏相关图选择合适的p,q进行试验和相应的统计量检验。因为ARMA模型和乘积型季节性模型不冲突，乘积型季节性模型是对ARMA模型的改进，根据ARMA 模型定性对乘积型季节性模型定阶。

由于原数据有明显的12个月的季节性，且数据有上升趋势，故应该对数据做一下差分运算：

对进行疏系数ARMA模型拟合。代表在自回归部分，

有非零自回归参数；表示在滑

动平均部分；有非零互动平均参数，用选取p，q的各种阶数形式进行，并采取AIC准则寻求最优模型。

由自相关图一阶截尾，十二阶明显大于2倍的标准差，且偏相关图有拖尾的性质，故选择MA（1），MA (12)，AR (1)，AR (2)，AR (12)，ARMA(1,1)，ARMA(1,2)，ARMA(1,12),ARMA(12,1)使用SAS软件T统计量检验参数的显著性。

例如：

将T统计量通过显著性检验的模型，对p，q，P，Q进行初步定阶。

选用对应的p，q，共确定一下五种模型：

①q=（1）（12）

②p=1 q=(1)(12)

③p=(12) q=(1)(12)

④q=(1,2)(12,24)

⑤p=(1) q=(2,12)

⑥p=(1) q=(1,12)

（注：第一个括号内表示低阶模型的p或者q的值，第二个括号内表示季节性的P或者Q。没有括号时表示低阶的p 或者q值。）

下面由PROC ARIMA过程(源程序见附件一)，得到如下结果，及其AIC参数：

①q=（1）（12）AIC=5542.897

②p=1 q=(1)(12) AIC=5521.385

③p=(12) q=(1)(12) AIC=5546.702

④q=(1,2)(12,24) AIC=5523.921

⑤p=(1) q=(2,12) AIC=5539.182

⑥p=(1) q=(1,12) AIC=5519.385

由AIC定阶准则，在以上六种方案中取AIC参数最小值为

⑦p=(1) q=(1,12) AIC=5519.385

模型的参数检验和残差的白噪声检验：

检验结果显示，该模型顺利通过参数检验和残差白噪声检验。

直观地看出该乘积模型对原序列的拟合效果良好。所以最后的选取的模型为：

ARIMA(0,1,(1,12)).

3.5.3 模型预测

由PROC ARIMA过程，预测未来12个月的平均气温，可以得到未来12个月的预报值及其预报标准差如下：

由上表可以得到未来十个月的预测的平均气温，分别为

-28.7205 6.0156 71.7717 153.1023

213.4104 255.2954 273.0346 260.9353

212.9610 142.5534 55.5982 -6.5335

3.5.4 误差分析

2010年1—12个月各月平均气温及误差分析

月份1月2月3月4月5月6月实测值（0.1℃）-48.0000 -10.0000 41.0000 112.0000 217.0000 247.0000 预测值（0.1℃）-28.7205 6.0156 71.7717 153.1023 213.4104 255.2954 绝对误差-19.2795 -16.0156 -30.7717 -41.1023 3.5896 -8.2954 相对误差40.17% 160.16% -75.05% -36.70% 1.65% -3.36%

月份7月8月9月10月11月12月实测值（0.1℃）286.0000 265.0000 213.0000 136.0000 58.0000 -10.0000 预测值（0.1℃）273.0346 260.9353 212.9610 142.5534 55.5982 -6.5335 绝对误差12.9654 4.0647 0.0390 -6.5534 2.4018 -3.4665 相对误差 4.53% 1.53% 0.02% -4.82% 4.14% 34.67%

对于2010年的月平均温度预报值，最大绝对误差出现在4月，最小绝对误差

出现在9月；最大相对误差出现在2月，最小相对误差出现在9月，其中3月份的相

对误差也较大。这是因为每年7月，8月和9月是北方夏天天气，气温较高，一般

在20几摄氏度，所以温度的波动较小，就可使相对误差较小。总的来说，模型拟

合效果不错。

3.5.5 预测模型的评价

由上面可以看出模型拟合和预测结果还是非常不错的，令为预测值，s为

预测值的标准差。几乎所有观测值最终都落在s的0.5摄氏度内。

ARIMA 乘积模型是针对有季节性变动的时间序列提出的建模方法。ARIMA 建

模法的最大特点之一就是其对序列平稳型的要求,在应用模型以前,对序列进行

预处理,以保证序列均值和方差的平稳性。在本例中,我们用一阶季节性差分和一

阶非季节性差分的方法,实现序列的平稳化。差分往往要求大样本的特性,即至少

7个季节周期的序列值。但在更多的实际应用中,我们仍需要不断地改进模型的

估计方法,主要是寻找对序列进行平稳性处理的方法,以提高预测模型的精度。

四、模型讨论

4.1 模型的改进：

本文拟通过建立时间序列模型，去发现时间序列中数据的趋势性，并通过构建合理的数学模型去使得数据平稳化，但是平稳化的过程遇到阻碍，很难使数据序列变得平稳，因此这时候就要想到如何去改进模型，于是进而利用到乘积性季节模型，通过乘积型季节模型的建立，对未来的气温做出一定误差范围里的合理预测。

对乘积型季节模型的定阶进行了改进，不是只通过观测自相关图和偏自相关图来定阶，而是用通过了T检验的ARMA模型进行定阶。通过了T检验的ARMA 模型的p，q，使得乘积型季节模型拟合得更好，使结果更准确。

4.2 模型的缺点：

尽管预测结果还不错，但是预测时的标准差比较差，使得预测值得取值范围比较大，导致有几个月份的相对误差偏大。所以该模型还有待于继续改进。

五、参考文献

1. 范金城,梅长林. 数据分析[M].北京: 科学出版社，2010.

2. 汪晓银,周保平.数学建模与数学实验[M].北京: 科学出版社,2010.

3．王燕.应用时间序列分析.[M].北京: 中国人民大学出版社,2008.

4. 张世英,许启发,周红.金融时间序列分析[M].北京: 清华大学出版社,2008

5. 潘红宇.金融时间序列分析.[M]北京: 对外经济贸易大学出版社,2008

6. 黄燕, 吴平等.SAS统计分析及应用[M].北京: 机械工业出版社,2006.

7. Margaret H.Dunham.数据挖掘教程[M].北京: 清华大学出版社,2005.

六、附件

附件一：

data sjfx;

input y@@;

date=intnx('month','31dec1950'd,_n_);

format date monyy.;

cards; 原数据略

;

proc print;

run;

proc gplot data=sjfx;

plot y*date;

symbol1i=join v=dot c=black;

run;

proc arima;

identify var=y(1,12) nlag=15;

estimate q=(1)(12) noconstant method=ml;

forecast lead=24interval=month id=date out=a;

estimate q=(1)(12) p=1method=ml;

forecast lead=24interval=month id=date out=a1;

estimate q=(1)(12) p=(12) method=ml;

forecast lead=24interval=month id=date out=a2;

estimate q=(1,2)(12,24) noconstant method=ml; forecast lead=24interval=month id=date out=a3;

estimate q=(2)(12) p=(1) noconstant method=ml; forecast lead=24interval=month id=date out=a7;

estimate q=(1)(12) p=(1) noconstant method=ml; forecast lead=24interval=month id=date out=a8;

run;

symbol2 v=plus i=join c=green;

symbol3 v=none i=join c=blue;

run;