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

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

资料数据见表1。

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

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

3.3.1 样本自相关函数平稳序列21012,,,,,ωωωωω--L L ， ()0t E ω=，对于样本，定义自协方差函数：112211ˆn kk k n k nk j k j j nn ωωωωωωγωω-++-+=+++==∑L ，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个样本自相关函数中至多有一个使ˆ2k ρ≥ 3、若样本自相关函数和样本偏相关函数都是拖尾的，则线性模型可以看成混和模型。

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

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

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

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

2、()MA q 滑动平均模型参数估计：22221221+1ˆˆˆˆ(1),0ˆˆˆˆˆˆˆ(),1qk k k q k q k k q αασθθθγσθθθθθ-⎧++++=⎪=⎨-+++≤≤⎪⎩L L 可得1q +个方程，求212ˆˆˆˆ,,qαθθθσL ，即解这个非线性方程组。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

111ˆˆ0.1695φφ==-，222ˆˆ0.0190φφ==-，由图-3的，0ˆ 1.6320e+004γ=，由公式201pj j j ασγφγ==-∑得：2ˆ 1.5855e+004ασ=图5 噪声方差的计算由上可知模型为：120.16950.0190t t t t ωωωα--++=，又知11402.82nj 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=[ ……];b=length(rainfall);z=sum(rainfall)/b; ………………………………计算均值w=rainfall-z; ………………………………由t Z 构造t ω序列 sumw=zeros(1,6); sumw1=0; for j=1:50sumw1=sumw1+w(j)^2; ..……………………………..计算0γ endfor k=0:5 for i=1:(b-k)sumw(k+1)=sumw(k+1)+w(i)*w(i+k); …………….......计算k γ end endr=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);Ma5=[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); ………………..计算55kk=zeros(1,5);kk=[kk11,kk22,kk33,kk44,kk55];σD=r0-kk11*r(2)-kk22*r(3) ………………..计算2α。