时间序列分析降水量预测模型

课程名称: 时间序列分析
题目: 降水量预测
院系：理学院
专业班级：数学与应用数学10-1
学号：
学生姓名：戴永红
指导教师：＿＿潘洁＿
2013年 12 月 13日
1.问题提出
能不能通过以前的降水序列为样本预测出2002的降水量？
2.选题
以国家黄河水利委员会建站的山西省河曲水文站1952年至2002年51年的资料为例，以1952年至2001年50年的降水序列作为样本，建立线性时间序列模型并预测2002年的降水状态与降水量，并与2002年的实际数据比较说明本模型的具体应用及预测效果。

资料数据见表1。

表1 山西省河曲水文站55年降水量时间序列
3．原理 3.1模型表示
均值为0，具有有理谱密度的平稳时间序列的线性随机模型的三种形式，描述如下：
1、()AR p 自回归模型：1122t t t p t p t ωφωφωφωα-------=由2p +个参数刻画；
2、()MA q 滑动平均模型：1122t t t t q t q ωαθαθαθα---=----由2q +个参数刻画；
3、(,)ARMA p q 混和模型：
(,)ARMA p q 混和模型由3p q ++个参数刻画；
3.2 自相关函数k ρ和偏相关函数kk φ
1、自相关函数k ρ刻画了任意两个时刻之间的关系，0/k k ργγ=
2、偏相关函数kk φ刻画了平稳序列任意一个长1k +的片段在中间值
11,t t k ωω++-固定的条件下，两端t ω，t k ω+的线性联系密切程度。

3、线性模型k ρ、kk φ的性质
表2 三种线性模型下相关函数性质
3.3 模型识别
通常平稳时间序列t Z ，0,1t =±仅进行有限n 次测量(50)n ≥，得到一个样本
函数，且利用平稳序列各态历经性：1
1n
j j Z Z n μ=≈=∑做变换，t t Z ω=，1,
t n =，将
1,,n Z Z 样本换算成为样本1,,n ωω，然后再确定平稳时间序列{,0,1}t t ω=±的随机
线性模型。

3.3.1 样本自相关函数
平稳序列
21012
,,,,,ωωωωω--， ()0t E ω=，对于样本，定义自协方差函数：
11221
1ˆn k
k k n k n
k j k j j n
n ωωωωωωγωω-++-+=++
+=
=∑，0ˆˆˆ/k k ρ
γγ=。

同时为了保证ˆk k γγ=，ˆk k ρ
ρ=一般取50,/4n k n ><。

常取/10k n =。

3.3.2 确定模型类别和阶数
在实际应用中，我们常用有一个样本算出的ˆk k ρρ=，ˆkk kk
φφ=判别k ρ，kk φ是拖尾还是截尾的。

随机线性模型的三种形式的判别分别如下：
1、若k ρ拖尾，kk φ截尾在k p =处，则线性模型为()AR p 模型。

k ρ拖尾可以用的点图判断，只要样本自相关函数的绝对值愈变愈小；当k p >时，平均20个样本偏
相关函数中至多有一个使ˆ2/kk
φ≥，则认为kk φ截尾在k p =处。

2、若kk φ截尾，
k ρ在k p =处截尾，那么线性模型为()MA q 滑动平均模型。

kk φ拖
尾可以根据样本偏相关函数的点图判断，只要ˆkk
φ愈变愈小。

当k q >时，若平均20个样本自相关函数中至多有一个使ˆ2/k ρ
≥ 3、若样本自相关函数和样本偏相关函数都是拖尾的，则线性模型可以看成混和
模型。

3.4 模型参数估计
1、()AR p 模型参数估计：
()AR p 模型有2p +个参数：2
12,,,,,p p α
φφφσ。

利用Yule-Walker 方程，利用Toeplitz 矩阵求逆和作矩阵乘法的方法算样本偏相关函数kk φ。

()AR p 模型的参数值不必作专门的计算，只要在样本偏相关函数计算的记录中取出样本参数值即可。

此时12,,,p φφφ,都已经确定了，经过推理我们可以得到：2
01p
j j j ασγφγ==-∑。

2、()MA q 滑动平均模型参数估计：
可得1q +个方程，求212ˆˆˆˆ,,q αθθθσ，即解这个非线性方程组。

3、(,)ARMA p q 混和模型参数估计
对于满足一个条件：1111......t t p t p t t p t q a a a ωφωφωθθ-------=---采用先计算
12ˆˆˆ,,,p φφφ，在计算212
ˆˆˆˆ,,q
αθθθσ的方法，具体如下：1）可利用Toeplitz 矩阵和作矩阵乘法的方法求出12ˆˆˆ,,,p φφφ。

2）令'11...t t t p t p
ωωφωφω--=---混和模型化为：'11...t t t p t q a a a ωθθ--==---这是关于't ω的()MA q 模型，用't ω的样本协方差函数估计
212
ˆˆˆˆ,,q αθθθσ的值。

4. 步骤
采用MATLAB 处理数据。

1、对一个时间序列做n 次测量得到一个样本函数12,,n Z Z Z 。

实验采用表1中
的降水量数据，50n =。

图1 山西省河曲水文站55年降水量时间序列
2、数据预先处理：做变换t t Z Z ω=-，其中50
1
150j j Z Z ==∑
图2 将时间序列变为期望为0的平稳时间序列
3、计算样本自协方差函数k γ，样本自方差函数k ρ。

0ˆˆˆ/k k ρ
γγ=，其中0,1,2,3,4,5k =，11221
1ˆn k
k k n k n
k j k j j n
n ωωωωωωγωω-++-+=++
+=
=∑。

由图-3
数据可得：随着k 的增大，k ρ越来越小，具有拖尾性。

图3 计算样本自相关函数
接下来计算偏相关函数kk φ(1k ≥)。

利用Yule-Walker 方程，利用Toeplitz
矩阵求逆和作矩阵乘法的方法算样本偏相关函数kk φ。

2/0.283=，由图-4得到的数据可得，2k p >=时，只有一个偏相关函数大于0.283。

所以确定阶数为：
2p =。

图4计算偏相关函数
5、由上综述：确定模型为(2)AR 模型。

下面进行(2)AR 模型参数的估计。

111ˆˆ0.1695φφ==-，222
ˆˆ0.0190φφ==-，由图-3的，0ˆ 1.6320e+004γ=，由公式2
01
p
j j j ασγφγ==-∑得：2
ˆ 1.5855e+004ασ=
图5 噪声方差的计算
由上可知模型为：12
0.16950.0190t t t t ωωωα--++=，又知1
1402.82n
j j Z Z n ===∑，
12402.820.1695(402.82)0.0190(402.82)t t t t Z Z Z α---+-+-=，2ˆ 1.5855e+004ασ
=。

最后确定(2)AR 模型为：
120.16950.0190478.75t t t t Z Z Z α--++=+，2ˆ 1.5855e+004ασ
= 6、通过确定的模型估计2002年的降水量
一步估计公式：1
ˆˆˆ(1)(1)0.16950.0190478.75k k k Z Z k Z Z -=+=--+。

其中，2001年的降水量为234.4mm ，2001年的降水量为289.6mm 。

20020.1695*234.40.0190*389.6478.75431.62Z =--+=mm
一步预报误差为79.66=mm ，而2002年实际降水量为487.3mm 。

为了提高预报准确度，可以提供更多样本点，进行预报估计。

5．部分程序代码及注释
rainfall=[261.6 ……389.6];
b=length(rainfall);
z=sum(rainfall)/b; ………………………………计算均值 w=rainfall-z; ………………………………由t Z 构造t ω序列
sumw=zeros(1,6); sumw1=0; for j=1:50
sumw1=sumw1+w(j)^2; ..……………………………..计算0γ end
for k=0:5 for i=1:(b-k)
sumw(k+1)=sumw(k+1)+w(i)*w(i+k); …………….......计算k γ end end
r=sumw/b; r0=sumw1/b;
p=r/r0; ……………………….计算自相关函数k ρ
kk11=p(2); ………………………计算11φ
a2=[1,p(2);p(2),1]
a22=inv(a2);
kk2=a22*p(1,2:3)'; ………………………计算
φ
22
kk22=kk2(2,1);
a5=[1,p(2),p(3),p(4),p(5);p(2),1,p(2),p(3),p(4);p(3),p(2),1,p(2),p(3 );p(4),p(3),p(2),1,p(2);p(5),p(4),p(3),p(2),1];
a55=inv(a5);
kk5=a55*p(1,2:6)';
kk55=kk5(5,1); ………………..计算
φ
55
kk=zeros(1,5);
kk=[kk11,kk22,kk33,kk44,kk55];
D=r0-kk11*r(2)-kk22*r(3) ………………..计算2
σ
α。