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

课程名称:时间序列分析题目:降水量预测院系：理学院专业班级：数学与应用数学10-1学号：2010303787学生姓名：戴永红指导教师：＿＿潘洁＿2013年12月13日1.问题提出能不能通过以前的降水序列为样本预测出2002的降水量？2.选题以国家黄河水利委员会建站的山西省河曲水文站1952年至2002年51年的资料为例，以1952年至2001年50年的降水序列作为样本，建立线性时间序列模型并预测2002年的降水状态与降水量，并与2002年的实际数据比较说明本模型的具体应用及预测效果。资料数据见表1。表1山西省河曲水文站55年降水量时间序列时段降水量(mm)时段降水量(mm)时段降水量(mm)1952195319541955195619571958195919601961196219631964196519661967196819691970261.6486.4631.5259.0568.0398.2479.6697.6397.7640.4247.1387.7694.2211.4322.6656.6325.3603.8424.81971197219731974197519761977197819791980198119821983198419851986198719881989383.3238.8423.0237.1330.7445.9518.9492.6490.3257.0400.6347.5368.3411.5356.2381.2318.0473.0373.31990199119921993199419951996199719981999200020012002369.0348.3469.2228.1338.8546.1358.9237.1423.3257.4234.4389.6487.33．原理3.1模型表示均值为0，具有有理谱密度的平稳时间序列的线性随机模型的三种形式，描述如下：1、()ARp自回归模型：1122tttptpt由2p个参数刻画；2、()MAq滑动平均模型：1122ttttqtq由2q个参数刻画；3、(,)ARMApq混和模型：11221122tttptptttqtq(,)ARMApq混和模型由3pq个参数刻画；3.2自相关函数k和偏相关函数kk1、自相关函数k刻画了任意两个时刻之间的关系，0/kk2、偏相关函数kk刻画了平稳序列任意一个长1k的片段在中间值11,ttk固定的条件下，两端t，tk的线性联系密切程度。3、线性模型k、kk的性质表2三种线性模型下相关函数性质模型函数()ARp()MAq(,)ARMApqk拖尾kq截尾拖尾kkkp截尾拖尾拖尾3.3模型识别通常平稳时间序列tZ，0,1t仅进行有限n次测量(50)n，得到一个样本函数，且利用平稳序列各态历经性：11njjZZn做变换，ttZ，1,tn，将1,,nZZ样本换算成为样本1,,n，然后再确定平稳时间序列{,0,1}tt的随机线性模型。3.3.1样本自相关函数平稳序列21012,,,,,，()0tE，对于样本，定义自协方差函数：112211ˆnkkknknkjkjjnn，0ˆˆˆ/kk。同时为了保证ˆkk，ˆkk一般取50,/4nkn。常取/10kn。3.3.2确定模型类别和阶数在实际应用中，我们常用有一个样本算出的ˆkk，ˆkkkk判别k，kk是拖尾还是截尾的。随机线性模型的三种形式的判别分别如下：1、若k拖尾，kk截尾在kp处，则线性模型为()ARp模型。k拖尾可以用的点图判断，只要样本自相关函数的绝对值愈变愈小；当kp时，平均20个样本偏相关函数中至多有一个使ˆ2/kkn，则认为kk截尾在kp处。2、若kk截尾，k在kp处截尾，那么线性模型为()MAq滑动平均模型。kk拖尾可以根据样本偏相关函数的点图判断，只要ˆkk愈变愈小。当kq时，若平均20个样本自相关函数中至多有一个使ˆ2/kn。3、若样本自相关函数和样本偏相关函数都是拖尾的，则线性模型可以看成混和模型。3.4模型参数估计1、()ARp模型参数估计：()ARp模型有2p个参数：212,,,,,pp。利用Yule-Walker方程，利用Toeplitz矩阵求逆和作矩阵乘法的方法算样本偏相关函数kk。()ARp模型的参数值不必作专门的计算，只要在样本偏相关函数计算的记录中取出样本参数值即可。此时12,,,p,都已经确定了，经过推理我们可以得到：201pjjj。2、()MAq滑动平均模型参数估计：22221221+1ˆˆˆˆ(1),0ˆˆˆˆˆˆˆ(),1qkkkqkqkkq可得1q个方程，求212ˆˆˆˆ,,q，即解这个非线性方程组。3、(,)ARMApq混和模型参数估计对于满足一个条件：1111......ttptpttptqaaa采用先计算12ˆˆˆ,,,p，在计算212ˆˆˆˆ,,q的方法，具体如下：1）可利用Toeplitz矩阵和作矩阵乘法的方法求出12ˆˆˆ,,,p。2）令'11...tttptp混和模型化为：'11...tttptqaaa这是关于't的()MAq模型，用't的样本协方差函数估计212ˆˆˆˆ,,q的值。4.步骤采用MATLAB处理数据。1、对一个时间序列做n次测量得到一个样本函数12,,nZZZ。实验采用表1中的降水量数据，50n。图1山西省河曲水文站55年降水量时间序列2、数据预先处理：做变换ttZZ，其中501150jjZZ图2将时间序列变为期望为0的平稳时间序列3、计算样本自协方差函数k，样本自方差函数k。0ˆˆˆ/kk，其中0,1,2,3,4,5k，112211ˆnkkknknkjkjjnn。由图-3数据可得：随着k的增大，k越来越小，具有拖尾性。图3计算样本自相关函数接下来计算偏相关函数kk(1k)。利用Yule-Walker方程，利用Toeplitz矩阵求逆和作矩阵乘法的方法算样本偏相关函数kk。2/500.283，由图-4得到的数据可得，2kp时，只有一个偏相关函数大于0.283。所以确定阶数为：2p。图4计算偏相关函数5、由上综述：确定模型为(2)AR模型。下面进行(2)AR模型参数的估计。111ˆˆ0.1695，222ˆˆ0.0190，由图-3的，0ˆ1.6320e+004，由公式201pjjj得：2ˆ1.5855e+004图5噪声方差的计算由上可知模型为：120.16950.0190tttt，又知11402.82njjZZn，12402.820.1695(402.82)0.0190(402.82)ttttZZZ，2ˆ1.5855e+004。最后确定(2)AR模型为：120.16950.0190478.75ttttZZZ，2ˆ1.5855e+0046、通过确定的模型估计2002年的降水量一步估计公式：1ˆˆˆ(1)(1)0.16950.0190478.75kkkZZkZZ。其中，2001年的降水量为234.4mm，2001年的降水量为289.6mm。20020.1695*234.40.0190*389.6478.75431.62Zmm一步预报误差为2ˆ279.66mm，而2002年实际降水量为487.3mm。为了提高预报准确度，可以提供更多样本点，进行预报估计。5．部分程序代码及注释rainfall=[261.6……389.6];b=length(rainfall);z=sum(rainfall)/b;………………………………计算均值w=rainfall-z;………………………………由tZ构造t序列sumw=zeros(1,6);sumw1=0;forj=1:50sumw1=sumw1+w(j)^2;..……………………………..计算0endfork=0:5fori=1:(b-k)sumw(k+1)=sumw(k+1)+w(i)*w(i+k);…………….......计算kendendr=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