时间序列的matlab实例编程

%时间序列ex1%x-72个源数据%各个图像需要单独的操作完成(手工进行)x=90078106 8928 9137 10017 10826 11317 10744 9713 9938 9161 8927 77506981 8038 8422 8714 9512 10120 9823 87439192 8710 8680 81627306 8124 7870 9387 9556 10093 9620 82858433 8160 8034 77177461 7776 7925 8634 8945 10078 9179 80378488 7874 8647 77926957 7726 8106 8890 9299 10625 9302 83148850 8265 8796 78366892 7791 8129 9115 9434 10484 9827 91109070 8633 9240;T1=zeros(1,72);t=9651.75 8723.75 8585.8+1/30 8396.75 8576.8+1/30 8796.75; %t-各年年平均死亡人数for i=1:6 T1(12*i-11):(12*i)=t(i); %趋势项T1的赋值end%数据图和分段趋势1% plot(x);% hold on% plot(T1,r);% hold off%method1-分段趋势法S1=zeros(1,72);for i=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);% hold on% plot(R1,r);% hold off%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);% hold on% plot(T2,r);% hold offS2=zeros(1,72);for i=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);% hold on% plot(R2,r);% hold off%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,Designed by luli.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)模型的自协方差函数列for i=1:3 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(max(3-i,i-1) r(i)*a;end%变换后的MA(2)模型的参数求解R=zeros(1,k+1);R(1,1:3)=ry;A=0 1;0 0;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);for i=1:k for j=1:k gma(i,j)=R(abs(i-j)+1); endendPI=(oma/gma)*oma;sig2=ry(1)-C*PI*C;b=(r2-A*PI*C)/sig2;function gammak=arma2_5(L,Max)%功能:计算ARMA(2,2)的自协方差函数列%输入:L-变量的计算下标;Max-计算自协方差的最大次数(默认10000)%输出:gamak-自协方差函数%2011-4-17，Designed by lulia1=0.0894; %初始参数的赋值a2=-0.6265;% b0=1;b1=-0.3334;b2=0.8158;sig2=4.0119;psy=eye(1,Max); %wold系数的初值psy(2)=b1+a1;psy(3)=b2+a2+a1*psy(2);for j=4:Max %wold系数列的计算 psy(j)=a1*psy(j-1)+a2*psy(j-2);endgammak=sig2*sum(psy(1:Max-L).*psy(L+1:Max); %自协方差函数的计算function v=MA2(r1,k)%功能:根据自协方差函数求MA(2)模型%输入:r1-已知自协方差函数序列;k-矩阵的阶数%输出:v-模型的参数组合%2011-4-17，Designed by lulir=zeros(1,k+1); %已知参数赋值r(1,1:3)=r1;A=0 1;0 0;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);for i=1:k %gma矩阵的赋值 for j=1:k gma(i,j)=r(abs(i-j)+1); endendPI=(oma/gma)*oma;sig2=r1(1)-C*PI*C;b=(r2-A*PI*C)/sig2;v=k b 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,Designed by luli.format short gH=zeros(5,length(N); %各参数的初始化for j=1:length(N) %不同观测数的循环模拟 Avemu=zeros(1,length(N); Avee=zeros(1,length(N); Stdmu=zeros(1,length(N); Stde=zeros(1,length(N); for k=1:M %多次模拟计算 Y=zeros(1,m+N(j); terr=zeros(1,m+N(j); mu=zeros(1,M); e=zeros(1,M); for l=3:m+N(j) %序列值的计算 terr(l)=randn; Y(l)=2*row*cos(theta)*Y(l-1)-row2*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 %输出类型的判断 if M=1 Avemu(j)=mu; %输出参数的计算 Avee(j)=e; Stdmu(j)=0; Stde(j)=0; else Avemu(j)=sum(mu)/M; Avee(j)=sum(e)/M; Stdmu(j)=sqrt(sum(mu-Avemu(j).2)/(M-1); Stde(j)=sqrt(sum(e-Avee(j).2)/(M-1); end H(:,j)=N(j);Avemu(j);Avee(j);Stdmu(j); Stde(j); %参数值的输出矩阵赋值Endfunction H=Simulink(m,N,M,row,theta)H=zeros(5,length(N);for j=1:length(N) for k=1:M Y=zeros(1,m+N(j); terr=zeros(1,m+N(j); X=zeros(1,N(j); err=zeros(1,N(j); for l=3:m+N(j) terr(l)=randn; Y(l)=2*row*cos(theta)*Y(l-1)-row2*Y(l-2)+terr(l); end X=Y(1,m+1:m+N); err=terr(1,m+1:m+N); mu(k)=mean(X); e(k)=mean(err)