Harbin Institute of Technology
实验报告
课程名称: 时间序列分析
设计题目: 降水量预测模型
院 系: 电子信息与工程学院
班 级: 电子二班
设 计 者:
学 号:
指导教师: 冀振元
设计时间: 2010/05/07
1. 实验选题
课程设计以国家黄河水利委员会建站的山西省河曲水文站1952年至2002年51年的资料为例,以1952年至2001年50年的降水序列作为样本,建立线性时间序列模型并预测2002年的降水状态与降水量,并与2002年的实际数据比较说明本模型的具体应用及预测效果。资料数据见表1。
表1 山西省河曲水文站55年降水量时间序列
2.实验原理
2.1模型表示
均值为0,具有有理谱密度的平稳时间序列的线性随机模型的三种形式,描述如下:
1、自回归模型:由个参数刻画;
2、滑动平均模型:由个参数刻画;
3、混和模型:
混和模型由个参数刻画;
2.2 自相关函数和偏相关函数
1、自相关函数刻画了任意两个时刻之间的关系,
2、偏相关函数刻画了平稳序列任意一个长的片段在中间值固定的条件下,两端,的线性联系密切程度。
3、线性模型、的性质
表2 三种线性模型下相关函数性质
2.3 模型识别
通常平稳时间序列,仅进行有限次测量,得到一个样本函数,且利用平稳序列各态历经性:做变换,,,将样本换算成为样本,然后再确定平稳时间序列的随机线性模型。
2.3.1 样本自相关函数
平稳序列,,对于样本,定义自协方差函数:
,。同时为了保证,一般取。常取。
2.3.2 确定模型类别和阶数
在实际应用中,我们常用有一个样本算出的,判别,是拖尾还是截尾的。随机线性模型的三种形式的判别分别如下:
1、若拖尾,截尾在处,则线性模型为模型。拖尾可以用的点图判断,只要样本自相关函数的绝对值愈变愈小;当时,平均20个样本偏相关函数中至多有一个使,则认为截尾在处。
2、若截尾,在处截尾,那么线性模型为滑动平均模型。拖尾可以根据样本偏相关函数的点图判断,只要愈变愈小。当时,若平均20个样本自相关函数中至多有一个使。
3、若样本自相关函数和样本偏相关函数都是拖尾的,则线性模型可以看成混和模型。
2.4 模型参数估计
1、模型参数估计:
模型有个参数:。利用Yule-Walker方程,利用Toeplitz矩阵求逆和作矩阵乘法的方法算样本偏相关函数。模型的参数值不必作专门的计算,只要在样本偏相关函数计算的记录中取出样本参数值即可。此时,都已经确定了,经过推理我们可以得到:。
2、滑动平均模型参数估计:
可得个方程,求,即解这个非线性方程组。
3、混和模型参数估计
对于满足一个条件:采用先计算,在计算的方法,具体如下:1)可利用Toeplitz矩阵和作矩阵乘法的方法求出。2)令混和模型化为:这是关于的模型,用的样本协方差函数估计的值。
3. 实验步骤
课程设计采用MATLAB处理数据。
1、对一个时间序列做次测量得到一个样本函数。实验采用表1中的降水量数据,。
图1 山西省河曲水文站55年降水量时间序列
2、数据预先处理:做变换,其中
图2 将时间序列变为期望为0的平稳时间序列
3、计算样本自协方差函数,样本自方差函数。,其中,。由图-3数据可得:随着的增大,越来越小,具有拖尾性。
图3 计算样本自相关函数
接下来计算偏相关函数()。利用Yule-Walker方程,利用Toeplitz矩阵求逆和作矩阵乘法的方法算样本偏相关函数。,由图-4得到的数据可得,时,只有一个偏相关函数大于0.283。所以确定阶数为:。
图4计算偏相关函数
5、由上综述:确定模型为模型。下面进行模型参数的估计。
,,由图-3的,,由公式得:
图5 噪声方差的计算
由上可知模型为:,又知,,。
最后确定模型为:
,
6、通过确定的模型估计2002年的降水量
一步估计公式:。其中,2001年的降水量为234.4mm,2001年的降水量为289.6mm。
mm
一步预报误差为mm,而2002年实际降水量为487.3mm。为了提高预报准确度,可以提供更多样本点,进行预报估计。
4.部分程序代码及注释
rainfall=[261.6 ……389.6];
b=length(rainfall);
z=sum(rainfall)/b; ………………………………计算均值
w=rainfall-z; ………………………………由构造序列
sumw=zeros(1,6);
sumw1=0;
for j=1:50
sumw1=sumw1+w(j)^2; ..……………………………..计算
end
for k=0:5
for i=1:(b-k)
sumw(k+1)=sumw(k+1)+w(i)*w(i+k); …………….......计算
end
end
r=sumw/b;
r0=sumw1/b;
p=r/r0; ……………………….计算自相关函数
kk11=p(2); ………………………计算
a2=[1,p(2);p(2),1]
a22=inv(a2);
kk2=a22*p(1,2:3)'; ………………………计算
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); ………………..计算
kk=zeros(1,5);
kk=[kk11,kk22,kk33,kk44,kk55];
D=r0-kk11*r(2)-kk22*r(3) ………………..计算
本文来源:https://www.2haoxitong.net/k/doc/348bf443a8956bec0975e371.html
文档为doc格式