模型的识别与预测VIP

模型的识别与预测实验内容依照某AR模型生成一段数据（1000），同时用另一MA模型生成一段数据（200），合成一段1200长度的数据1）依赖于这1200个数据的前800个数据，识别这段数据背后的AR模型。2）在1）的基础上对新数据进行预测，并通过后续的400个数据进行判别（数据模型是否匹配）或者模型的修正（修正只需要提供思路和方法）。理论基础1.时间序列模型介绍时间序列是随时间改变而随机地变化的序列。时间序列分析的目的是找出它的变化规律，即线性模型，主要有三种：AR模型（自回归模型）、MA模型（滑动平均模型）和ARMA模型（自回归滑动平均模型或混合模型）。设{Xt}为零均值的实平稳时间序列，阶数为p的AR模型定义为其其中{}成为自回归系数，白噪声序列{}成为新信息序列；阶数为q的MA模型定义为其中{}称为滑动平均系数；P阶自回归q阶ARMA模型定义为记为ARMA（p，q）。模型的识别根据教材对平稳时间序列的特性分析，对初步识别平稳时间序列的类型提供了依据，如表1所示：表1各时间序列模型的特性类别模型AR（p）MA（q）ARMA（p，q）自相关函数拖尾截尾k=q处拖尾偏相关函数截尾k=p处拖尾拖尾模型阶数的确定样本自相关函数和样本偏相关函数设有零均值平稳时间序列{}的一段样本观测值，样本协方差函数估计式为同理样本自相关函数定义为MA模型阶数的确定设{}是正态的零均值平稳MA（q）序列，而对于充分大的N，可以认为的分布近似于正态分布，从而，的截尾性判断如下：首先计算（取），因为q值未知，故令q值从小到大，分别检验满足或的比例是否占总个数M的68.3%或95.5%。第一个满足上述条件的q就是的截尾处，及MA（q）模型的阶数。模型参数的估计当选定模型及确定阶数后，进一步的问题就是估计出模型的位置参数。参数估计的方法有矩法，最小二乘法及极大似然估计等。在工程计算中，只要用样本自协方差函数或样本自相关函数中的一部分数值。这里不详细介绍。模型的检验由样本观测序列{}，经过模型的识别，阶数的确定和参数估计，可以初步建立{}的模型。这样建立的模型一般还需要进行统计检验，只有经检验确认模型基本上能反映{}的统计特性时，用它进行预测才能获得良好的效果。模型的所谓自相关函数检验法，其基本思想是，如果模型是正确的，则模型的估计值与实际观测值所产生的残差序列应是随机干扰产生的误差，即{}应是白噪声序列。否则，模型不正确。模型的预报根据时间序列观测数据，建立了一个与实际问题相适应的模型后，就可以利用过去与现在的观测值，对该序列未来时刻的取值进行估计，即预报。关于预报效果优劣的标准，下面采用的是在平稳性最小方差意义下的预报。具体步骤从AR模型中提取1000个样本点以及MA模型中提取200个样本点，分别存储于x2和x1中；取x2的前800个数据作为模型识别的样本观测值，存储于X2中，计算样本自相关函数和样本偏相关函数，并通过它们的特性识别模型；根据自相关或者偏相关确定模型阶数，用最小二乘法求模型系数；提取1000个数据中后200个数样本据检验所得模型是否正确，用卡方假设检验，其自由度为M，显著性水平为0.05；预测模型的新数据，并与原始数据相比较；通过后续的400个数据进行判别（数据模型是否匹配）结论程序基本满足实验要求，能通过一组样本数据通过计算自相关函数和偏相关函数来识别其背后的模型，并确定模型的阶数以及模型的系数，得到一个完整的时间序列。在此基础上，可以通过该数据模型预报一组新的数据，与原始数据相比，能有较好的预报性，还能判别一组数据与这个模型是否匹配。我关于模型修正的思路：为将这后续400个数据与再从原始的800个数据提取出的400个数据合并，重新进行一次模型的建立，理论上能得到较好的修正模型。分析总结尽管程序已基本完成实验内容，但还是有一些不尽人意的地方，还存在如下问题：已给AR模型为2阶，程序中未考虑到更高阶的情况，即使识别到的是更高阶模型，也作二阶来处理，会存在误差，这点有待完善；在确定模型阶数时，偏相关函数的计算不够准确也可能导致阶数的判定出错，但在后续加了检验程序可在一定程度上解决该问题；对于一阶情况也未作处理，样本数据被判别为别的模型时也未处理，这点可以后续完善。附录clearall;clc;%从MA模型X(t)=a(t)-1.3a(t-1)+0.4a(t-2)中提取200个样本点K=200;a1=randn(1,K);x1=zeros(1,K);x1(1)=1;x1(2)=1;fori=3:Kx1(i)=a1(i)-1.3*a1(i-1)+0.4*a1(i-2);endX1=x1;%从AR模型X(t)+0.8X(t-1)=a(t)中提取1000个样本点N=1000;a2=randn(1,N);x2=zeros(1,N);x2(1)=0.1;x2(2)=0.1;x2(3)=1;fori=3:Nx2(i)=1.5*x2(i-1)-0.5*x2(i-2)+a2(i);end%取x2的前800个数据作为模型识别的样本观测值L=800;X2=x2(1,1:L);%计算样本自相关函数和样本偏相关函数r=zeros(1,L);%样本协方差矩阵fork=0:L-1%求样本协方差rk=0;fori=1:L-krk=rk+X2(i)*X2(i+k);endr(k+1)=(1/L)*rk;endp=zeros(1,L);fori=1:L%求样本自相关函数p(i)=r(i)/r(1);endM=L/10;pp=zeros(M,M);%样本偏相关系数矩阵pp(1,1)=p(2);p(1)=[];fork=1:M-1%求样本偏相关函数s1=0;s2=0;forj=1:ks1=s1+p(k+1-j)*pp(k,j);s2=s2+p(j)*pp(k,j);endpp(k+1,k+1)=(p(k+1)-s1)/(1-s2);forj=1:kpp(k+1,j)=pp(k,j)-pp(k+1,k+1)*pp(k,k+1-j);endendpp1=diag(pp);%提取偏相关系数subplot(2,1,1)plot(p);title('自相关图');subplot(2,1,2)plot(pp1);title('偏相关图');%识别模型，返回模型的阶数%取M=L/10个自相关和偏相关系数count=0;forf1=1:Mifabs(p(f1))<=1/sqrt(L)count=count+1;endifcount/M>=0.683break;endendcount_1=0;forf2=1:Mifabs(pp1(f2))<=1/sqrt(L)count_1=count_1+1;endif(count_1)/M>=0.683breakendendiff1==M&&f2<Mdisp('建立AR模型');%确定模型阶数fori=1:Mifabs(pp1(i))<=1/sqrt(L)jieshu=i-1;breakendendfprintf('其阶数为:%d\n',jieshu);%最小二乘法求模型系数forn=1:jieshuY=X2';Y(1:n)=[];m=L-n;X=[];fori=1:mforj=1:nX(i,j)=X2(n+i-j);endendxishu=inv(X'*X)*(X'*Y);end%用后200个数据检验所得模型是否正确X3=x2(1,L+1:end);%提取后200个样本点T=N-L;a3=zeros(1,T);%由观测样本点计算得到的白噪声点ra=ones(1,T+1);%白噪声的协方差a3(1)=X3(1);a3(2)=X3(2)-xishu(1)*X3(1);ifjieshu>1fori=1:T-2%用模型估计值与样本观测值产生白噪声序列a3(i+2)=X3(i+2)-xishu(1)*X3(i+1)-xishu(2)*X3(i);end%求白噪声协方差M=T/10;fori=0:Ms3=0;forj=1:T-is3=s3+a3(j)*a3(j+i);endra(i+1)=s3/T;end%求白噪声自相关rra=ones(1,M);fori=1:Mrra(i)=ra(i+1)/ra(1);end%卡方假设检验，自由度为M，显著性水平为0.05Qm=0;fori=1:MQm=Qm+(rra(i))^2;endQm=T*Qm;kafang=chi2inv(0.95,M);%求得卡方值ifQm<kafangdisp(['统计量','Qm=',num2str(Qm),'小于给定水平的卡方值，故所估计模型合适']);elsedisp(['统计量','Qm=',num2str(Qm),'大于给定水平的卡方值，故所估计模型不合适,请继续执行...']);endendendiff1==M&&f2==Mdisp('建立ARMA模型');endiff1<M&&f2<Mdisp('初始数据处理出错，重新执行');endiff1<M&&f2==Mdisp('建立MA模型');fori=1:Mif(abs(p(i)))<=1/sqrt(L)jieshu=i-1;breakendendfprintf('其阶数为:%d\n',jieshu);end%预测模型的新数据new_X=ones(1,200);new_X(1)=X2(end);new_X(2)=xishu(2)*X2(end)+xishu(1)*new_X(1);fori=1:198new_X(i+2)=xishu(1)*new_X(i+1)+xishu(2)*new_X(i);endfiguresubplot(2,1,1)plot(1:N,x2,'r.')title('原始数据')subplot(2,1,2)plot(1:200,new_X,'b*')holdplot(1:200,X3(1:200),'r.')legend('预测点','后200点')title('预测数据与模型生成后200样本点比较')%通过后续的400个数据进行判别（数据模型是否匹配）X4(1:200)=X1;X4(201:400)=X3;%提取后续400个数据U=400;a4=zeros(1,U);%由观测样本点计算得到的白噪声点ra1=ones(1,U+1);%白噪声的协方差a4(1)=X4(1);a4(2)=X4(2)-xishu(1)*X4(1);ifjieshu>1fori=1:U-2%用模型估计值与样本观测值产生白噪声序列a4(i+2)=X4(i+2)-xishu(1)*X4(i+1)-xishu(2)*X4(i);end%求白噪声协方差M=U/10;fori=0:Ms4=0;forj=1:U-is4=s4+a4(j)*a4(j+i);endra1(i+1)=s4/U;end%求白噪声自相关rra1=ones(1,M);fori=1