随机水文学第六章

第六章自回归滑动平均模型主要内容：第一节滑动平均模型第二节自回归滑动平均模型第三节模型的识别和检验第一节滑动平均模型一、一般滑动平均模型MA(q)滑动平均模型的一般形式MA(q)其中心化形式为：qtqqtqtqttqx,22,11,式中：θq,1,θq,2,…,θq,q——各阶权重又称权重系数是待求的参数εt——白噪声系列，在水文上仍假定为P-Ⅲ分布，均值为零，其参数仅有两个即标准差σε与偏态系数，没有Cvε二、参数估计2,222,,1,3,2,2221,,2,3,1,2,1,2122,21,2)()()1(qqxqqqqqqqqxqqqqqqqqqxqqqxrrr式中：样本序列的和r1,r2,…,rq（各阶自相关系数）均可由样本计算出。由一般形式经推导可得到如下方程组：三、的计算Cs一般形式经变形可得：qtqqtqtqttx,22,11,其中：x1’,x2’,…,xn’为已知中心化径流系列θq,1,θq,2,…,θq,q，已由上面解出'11x11,22'qx12,21,33'qqx11,22,11,'nqnqnqnnx………则三、的计算（2）利用这样计算得到εt系列，系列ε1，ε2，…，εn按照矩法估计：Cs313)3/(nCsnii四、一阶滑动平均模型以中心化变量表示的MA(1)模型形式如下：11,1tttxCs可由上述方法计算得到五、二阶滑动平均模型以中心化变量表示的MA(2)模型形式如下：22,211,2ttttx1))21(214121()21(214121221222221222222rrrrrrrr，θ)1(1,23,222,21θθθrr2,2222θr其中当r20时取+，当r20时取-六、例：同前，江西塞塘站有15年（1964-1978年）年径流资料，试建立MA（1）模型，并生成10年符合P—Ⅲ型分布年径流系列。例计算步骤11,1'tttxxxxtt'1．MA（1）模型的一般形式2．计算样本统计参数同前94.789x19.1xCs03.245x例计算步骤3．用公式计算参数347.01r6659424112122r404.041122111,1rr例计算步骤4．由实测样本经中心化后分别计算εt值04.61'11x855.113'11,122x340.157'31,144x642.129'11,1nnnx844.91'21,133x……例计算步骤5．计算Csε313)3/(nCsnii6．模型的确定1404.0'tttx1404.094.789tttx例计算步骤7．由模型模拟10项新系列11,1'tttx11,1tttx(1)先生成0~1均匀分布伪随机数10个作为P值(2)由上述计算Csε查P—Ⅲ离均系数Φ值表，得到相应的Φ值(3)代入上面公式可生成xt+1同理生成xt+2,…,xt+10计算结束tt第二节自回归滑动平均模型一、一般自回归滑动平均模型一般自回归滑动平均模型ARMA(p,q)，若以中心化变量xt’来表示，则：qtqqtqtqtptpptptptxxxx,22,11,,22,11,''''xxxtt'当时，即q=0时为自回归模型AR(p)0,2,1,qqqq当时，即p=0时为自回归模型MA(q)0,2,1,pppp二、ARMA(1,1)模型的参数估计方法11,111,1''ttttxx1,11221,11,122121,11,11,11,1212))(1(rrrrxxx经推算可得到如下方程：niiniiixxxxxxr122122)())((1)(122nxxniiniiniiixxxxxxr121111)())((方程中二、ARMA(1,1)模型的参数估计方法则前式可由上述三个用实测样本估计的值算出Φ1,1，θ1,1，σε21,121221,11,11,1211,1221,1121,1))(1(xxxxrrrr可采用图解法或迭代法求解θ1,1，σε2三、Csε的计算将一般形式变形为：11,111,1''ttttxx仿照MA(q)的计算εt系列的方法，由实测的系列x1’,x2’,…,xn’（中心化后的）及上面已计算出的Φ1,1,θ1,1，利用上式计算出ε1，ε2，，εn系列，再利用Csε的矩法公式：33)3/(nCsi计算Csε第三节模型的识别和检验对给定的水文序列，应选什么样的模型，是AR(p)，MA(q)还是ARMA(p,q)？若选出一种模型，例如AR(p)模型，又如何确定阶数P。这样的问题在时间序列分析中叫做模型类型的选择及其形式的确定，可以统称为模型的识别。通过识别确定模型，在由实测系列估计出模型中的参数，我们就初步得到了表征其随机变化的模型。对这种模型，尚需作进一步检验，以验证是否符合模型的一些基本假定。一、模型初步认别的原理和方法在模型初步认别中，主要的依据是样本序列的自相关函数和偏相关函数。主要特点归纳如下：(1)AR(p)序列自相关函数ρk，随滞时k的增大逐步变小，自相关图呈拖尾状。它的偏相关函数Φk,k呈截尾状，单调或波动衰减趋向于零。在k=p处出现一个截止点，即在时，当kp时Φk,k=0。pk0,kk(1)AR(p)序列（如下图所示）-1-0.500.51012345678910KρkAR(1)00.40.80123456φk,kKAR(2)-0.200.20.40123456φk,kK(2)MA(q)序列自相关函数ρk呈截尾状，k=q时出现一个截尾点，即在时；当kq时ρk=0。相反，它的偏相关函数随阶数的增加而逐渐变小，呈拖尾状。单调或波动衰减趋向于零。qk0k-0.5-0.2500.250.50.751012345678Kφk,k(3)ARMA(p,q)序列自相关函数和偏相关函数都没有截尾点，均以拖尾状而逐渐变小，趋向于零。-1-0.75-0.5-0.2500.250.50.751024681012kρk-1-0.75-0.5-0.2500.250.50.751024681012kφk,k缺陷：但由样本序列去估计自相关函数ρk和偏相关函数Φk,k抽样误差较大，难以直观判断，必须进行统计推断。例如，即使是出自AR(p)模型的序列，由于抽样的缘故在kp时由样本序列估计出的也不全会为零，而在零上下起伏，而用统计检验方法，推断出它们和零的差异并不显著，则可认为截尾点是k=p，从而推断序列为AR(p)序列。下面分别叙述MA(q)，AR(p)和ARMA(p,q)模型识别的具体方法：1．MA(q)模型由于kq时ρk=0，rk在kq时渐进服从的正态分布（其中为rk的方差D(rk)），根据正态分布的性质：))21(1,0(12qllrnN)21(112qllrn%3.68)21(12112qllkrnrP%0.95)21(96.12112qllkrnrP1．MA(q)模型先假定q0如kq0，rk应明显不为零（k=1,2,…,q0）而当q=q0时，k分别取q0+1，q0+2，…，q0+M（M可取左右）代入上式不等式，如果满足上式不等式，则统计满足个数占总数的比例，如果占总数M的68.3%或95%左右则可判断rk在q0处截尾，并初步认为序列是MA(q0)序列。2．AR(p)模型对AR(p)序列，当kp时Φk,k=0；当kp时将渐近服从正态分布，为的方差D()。kk,ˆkk,ˆ)1,0(nNn1%3.681ˆ,nPkk%0.9596.1ˆ,nPkk同样找到p0，在p0前（k取1,2,…,p0）明显不为零，而p=p0时，k分别取p0+1,p0+2,…,p0+M（M可取左右）代入上式不等式统计满足上式不等式的个数占总数M的比例，如在68.3%或95.0%左右可判断在p0处截尾，并初步认为序列是AR(p0)序列。kk,ˆ3．ARMA(p,q)模型对于这种混合模型，自相关函数和偏相关函数均无截尾的性质，认别较困难，一般识别p,q的方法可以从低阶到高阶逐个尝试，(p,q)分别取(1,1)，(1,2)，(2,1)，(2,2)然后分别进行参数估计，定出模型在检验这个模型是否被接受。若被接受，则模型确定；否则重新调整(p,q)直到接受为止。二、模型识别的AIC准则日本的赤池（Akaike）对于ARMA(p,q)模型中阶数p和q的确定提出了AIC准则：)(2)ln(),(2qpnqpAIC式中：n为序列的长度；ln为自然对数，为残差的方差2q=0为AR(p)模型p=0为MA(q)模型使AIC达到最小值的模型被认为是可以接受的好模型。对AR(p)模型，按AIC准则识别步骤：（1）计算样本序列的自相关函数rknttkntkttkxxxxxxr121)())((（2）按递推计算自回归系数Φp,1,Φp,2,…,Φp,p按AIC准则识别步骤：（3）计算残差方差21RxppppprrrR,22,11,2（4）计算AICPnPAIC2)ln()(2（5）根据不同P计算出相应的AIC(P)使AIC达到最小的P为我们所求例有序列n=45采用AR(p)模型，分别计算得到AIC(k)值，如下：KAICKAICKAIC11.963161.8632111.42321.9271.8211121.333231.899281.7826131.302641.897391.6956141.285651.8895101.5329151.4021试判断AR(p)的阶数pp=14三、模型的检验原假设H0：εt序列相互独立统计量mkkrnQ12)(M——最大滞时，一般m取n/4左右Q服从自由度为（m-p-q）的Χ2分布由样本ε1,ε2,…,εm计算出Q值，取α值由Χ2表查出Χα2，若则接受原假设，认为εt为独立序列的假设成立，反之不成立。例某站有44年径流序列，已识别出AR（1）模型。由样本序列，通过该模型算得εt序列，进而算出它的自相关关系：r1(ε)=-0.13，r2(ε)=0.14，r3(ε)=0.14，r4(ε)=-0.09，r5(ε)=0.04，r6(ε)=0.18，r7(ε)=-0.18，r8(ε)=0.05，r9(ε)=0.02，r10(ε)=-0.09。四、模型实用性的初步分析选定的模型是否反映随机序列真实的统计特性，还必须在模型检验的基础上，作进一步的分析，即实用性分析。若建立的模型是为了随机模拟，这是由模型而得的大量模拟序列，应保持实测（样本）序列的主要统计特性。它们包括均值、方差（或变差系数）、偏态系数和一阶自相关系数。分析方法：1.长序列法由模型模拟出一个很长的模拟序列，例如模拟出长度N为10000，然后仅根据这一个长序列计算各项特征值，并与由实测系列计算的各项特征值进行比较。2.短序列法由模型模拟出许多等长度的短序列，其长度和实测序列长度一样。2.短序列法设有长度为n的m个模拟序列：第一组模拟序列x1,1,x1,2,…,x1,n第二组模拟序列x2,1,x