时间序列的matlab实例编程

%时间序列ex1%x-72个源数据%各个图像需要单独的操作完成(手工进行)x=[900781068928913710017108261131710744971399389161892777506981803884228714951210120982387439192871086808162730681247870938795561009396208285843381608034771774617776792586348945100789179803784887874864777926957772681068890929910625930283148850826587967836689277918129911594341048498279110907086339240];T1=zeros(1,72);t=[+1/30+1/30];%t-各年年平均死亡人数fori=1:6T1((12*i-11):(12*i))=t(i);%趋势项T1的赋值end%数据图和分段趋势1%plot(x);%holdon%plot(T1,'r');%holdoff%method1-分段趋势法S1=zeros(1,72);fori=1:12%季节项S1的赋值S1(i:12:(i+12*5))=sum(x(i+12.*[0:5])-T1(i+12.*[0:5]))/6;endR1=x-T1-S1;%随机项R1及绘图%plot(S1);%holdon%plot(R1,'r');%holdoff%method2-回归直线法Y=ones(2,72);Y(2,:)=1:72;A=inv(Y*Y')*Y*x';%回归系数T2=A(1)+A(2).*[1:72];%趋势项T2的赋值%数据图和直线趋势2%plot(x);%holdon%plot(T2,'r');%holdoffS2=zeros(1,72);fori=1:12%季节项S2的赋值S2(i:12:(i+12*5))=sum(x(i+12.*[0:5])-T2(i+12.*[0:5]))/6;endR2=x-T2-S2;%随机项R2及绘图%plot(S2);%holdon%plot(R2,'r');%holdoff%forecastfcT=A(1)+A(2).*[73:84];%趋势项预测fcS=S2(1:12);%季节项预测fcx=fcT+fcS;%死亡人数预测function[a,b,sig2]=arma2_4(r,k)%功能：根据自协方差函数列求ARMA(2,2)模型%输入:r-已知自协方差函数列，k-矩阵的阶数%输出：a-模型的数值项系数;b-噪声项系数;sig2-噪声项的方差%2011-4-17,Designedbyluli.atemp=[r(3)r(2);r(4)r(3)]\[r(4);r(5)];%系数a的计算a=[-1;atemp(1);atemp(2)];ry=zeros(1,3);%变换后的MA(2)模型的自协方差函数列fori=1：3r(max(3-i,i-1))ry(i)=a'*[r(i)r(i+1)r(i+2);r(max(3-i,i-1))r(i)r(i+1);r(4-i)r(i)]*a;end%变换后的MA(2)模型的参数求解R=zeros(1,k+1);R(1,1：3)=ry;A=[01;00];C=[1;0];oma=zeros(2,k);oma(1,1)=ry(2);oma(1,2)=ry(3);oma(2,1)=ry(3);r2=[ry(2);ry(3)];gma=zeros(k,k);fori=1：kforj=1：kgma(i,j)=R(abs(i-j)+1);endendPI=(oma/gma)*oma';sig2=ry(1)-C'*PI*C;b=(r2-A*PI*C)/sig2;functiongammak=arma2_5(L,Max)%功能：计算ARMA(2,2)的自协方差函数列%俞入:L-变量的计算下标；Max-计算自协方差的最大次数(默认10000)%输出：gamak-自协方差函数%2011-4-17，Designedbylulia1=；%初始参数的赋值a2=；%b0=1；b1=；b2=；sig2=；psy=eye(1,Max)；%wold系数的初值psy(2)=b1+a1；psy(3)=b2+a2+a1*psy(2)；forj=4:Max%wold系数列的计算psy(j)=a1*psy(j-1)+a2*psy(j-2)；endgammak=sig2*sum(psy(1:Max-L).*psy(L+1:Max))；%自协方差函数的计算functionv=MA2(r1,k)%功能:根据自协方差函数求MA(2)模型%俞入:r1-已知自协方差函数序列；k-矩阵的阶数%输出:v-模型的参数组合%2011-4-17，Designedbylulir=zeros(1,k+1);%已知参数赋值r(1,1:3)=r1;A=[01;00];C=[1;0];oma=zeros(2,k);oma(1,1)=r1(2);oma(1,2)=r1(3);oma(2,1)=r1(3);r2=[r1(2);r1(3)];gma=zeros(k,k);fori=1:k%gma矩阵的赋值forj=1:kgma(i,j)=r(abs(i-j)+1);endendPI=(oma/gma)*oma';sig2=r1(1)-C'*PI*C;b=(r2-A*PI*C)/sig2;v=[kb'sig2];%输出向量function[H]=AR2Simulink(m,N,M,row,theta)%功能:AR(2)的模拟计算.%输入:m-X序列的延后取值数;N-观测数据数;M-模拟次数;row-参数值1;theta-参数值2.%俞出：H-H的第一行是样本数向量N;H的第二行是估计的样本均值Avemu;H的第三行是估计的噪%声项均值Avee;H的第四行是估计的样本方差Stdmu;H的第五行是估计的噪声项方差Stde.%2011-5-1,Designedbyluli.formatshortgH=zeros(5,length(N));%forj=1:length(N)%Avemu=zeros(1,length(N));Avee=zeros(1,length(N));Stdmu=zeros(1,length(N));Stde=zeros(1,length(N));各参数的初始化不同观测数的循环模拟fork=1:M%Y=zeros(1,m+N(j));terr=zeros(1,m+N(j));mu=zeros(1,M);e=zeros(1,M);forl=3:m+N(j)%terr(l)=randn;Y(l)=2*row*cos(theta)*Y(l-1)-rowA2*Y(l-2)+terr(l);end%各参数的赋值多次模拟计算序列值的计算X=Y(1,m+1:m+N(j));err=terr(1,m+1:m+N(j));mu(k)=mean(X);e(k)=mean(err);end%输出类型的判断ifM==1Avemu(j)=mu;%输出参数的计算Avee(j)=e;Stdmu(j)=0;Stde(j)=0;elseAvemu(j)=sum(mu)/M;Avee(j)=sum(e)/M;Stdmu(j)=sqrt(sum((mu-Avemu(j)).A2)/(M-1));Stde(j)=sqrt(sum((e-Avee(j)).A2)/(M-1));参数值的输出矩阵赋值endH(:,j)=[N(j);Avemu(j);Avee(j);Stdmu(j);Stde(j)];%Endfunction[H]=Simulink(m,N,M,row,theta)H=zeros(5,length(N));forj=1:length(N)fork=1:MY=zeros(1,m+N(j));terr=zeros(1,m+N(j));X=zeros(1,N(j));err=zeros(1,N(j));forl=3:m+N(j)terr(l)=randn;Y(l)=2*row*cos(theta)*Y(l-1)-rowA2*Y(l-2)+terr(l);endX=Y(1,m+1:m+N);err=terr(1,m+1:m+N);mu(k)=mean(X);e(k)=mean(err);endifM==1Avemu=mu;Avee=e;Stdmu=0;Stde=0;elseAvemu=sum(mu)/M;Avee=sum(e)/M;Stdmu=sqrt(sum((mu-Avemu)42)/(M-1));Stde=sqrt(sum((e-Avee)42)/(M-1));endH(:,j)=[N(j);Avemu;Avee;Stdmu;Stde];end