时间序列的matlab实例编程.doc

%时间序列ex1 %x-72个源数据 %各个图像需要单独的操作完成(手工进行) x=[9007 8106 8928 9137 10017 10826 11317 10744 9713 9938 9161 8927 7750 6981 8038 8422 8714 9512 10120 9823 8743 9192 8710 8680 8162 7306 8124 7870 9387 9556 10093 9620 8285 8433 8160 8034 7717 7461 7776 7925 8634 8945 10078 9179 8037 8488 7874 8647 7792 6957 7726 8106 8890 9299 10625 9302 8314 8850 8265 8796 7836 6892 7791 8129 9115 9434 10484 9827 9110 9070 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; end R1=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 off S2=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; end R2=x-T2-S2; %随机项R2及绘图 % plot(S2); % hold on % plot(R2,r); % hold off %forecast fcT=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); end end PI=(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-自协