arma模型程序.docx

1、arma模型程序AR MA模型程序(总6页)clear,clcclose allset(0,defaultaxeslinestyleorder,:+,-s); %设置默认的线性和标识符set(0,defaultaxescolororder,0,0,1)%设置默认的线条颜色x=-pi:pi/10:pi;y=tan(sin(x)-sin(tan(x);set(gca,xtick,min(x),max(x) %设置x轴的范围set(gca,ylim, min(y) ,max(y),layer,top)%设置y轴的范围h=plot(x,y,markersize,9); %画图并且取得二维图形中最大值和

2、最小值的索引号grid onset(gca,layer,top) %将栅格放在最上层，防止遮盖axis equal %使每一格的刻度相同，同时也使长度一样xlabel(-pileq theta leqpi,fontsize,20)ylabel(tan(sin(theta)-sin(tan(theta),fontsize,20)title(it-pileq thetaleqpi的三角函数图,fontsize,20)text(0,0,bulletit原点,fontsize,25) %文本的精确定位%标注曲线的最大值和最小值x=get(h,xdata); %获得二维曲线的数据y=get(h,ydat

3、a);imin=find(min(y)=y);imax=find(max(y)=y);text(x(imin),y(imin),最小值=,num2str(y(imin),.horizontalalignment,center,verticalalignment,bottom,fontsize,17)text(x(imax),y(imax),最大值=,num2str(y(imax),.horizontalalignment,center,verticalalignment,top,fontsize,17)宋哲 2011/8/16 11:07:42%AR,MA模型clear,clcclose al

4、la=load();a=a(2,4,:);a=a;a=a(:);%将数据税收展成按时间序列排列的% plot(a,-p)b=diff(a);% plot(b,-p)r1=autocorr(b)%计算自相关系数r2=parcorr(b) %计算偏自相关系数figure, subplot(121),autocorr(b)subplot(122), parcorr(b)cs1=ar(b,2) %AR模型，b必须是列向量bhat=predict(cs1,b;20) %Discrete-time IDPOLY model:%A(q) = 1 - q-1 - q-2 b15=bhat1(end)bhat2

5、=predict(cs1,b;b15;-1)宋哲 2011/8/16 11:27:00%用逻辑斯特预测clear,clca=load();a=a(2,4,:);a=a;a=a(:);%将数据税收展成按时间序列排列的x0=a(1);tt=1:14;t0=1;xt=(cs,t)cs(1)/(1+(cs(1)/x0-1)*exp(-cs(2)*(t-t0);cs=lsqcurvefit(xt,rand(2,1),tt(2,end),a(2:end)%以上程序为三维视图和等高线画图程序clc, cleara=load();amin=min(min(a), amax=max(max(a)x0=0:5:4

6、000; y0=0:5:3000;xx0,yy0=meshgrid(x0,y0);mesh(xx0,yy0,a)v=-49,1,12,44,75,107,171,184:80:357,357:100:614,614;figure, c=contour(x0,y0,a,v); clabel(c)宋哲 2011-08-17 9:32:44clc, clear %该程序计算最小面积a=load();a=a;x0=0:5:4000; y0=0:5:3000;pp=csape(x0,y0,a)x=0:4000; y=0:3000;z=fnval(pp,x,y);m,n=size(z);s=0;for i

7、=1:m-1 for j=1:n-1 p1=x(i),y(j),z(i,j); p2=x(i+1),y(j),z(i+1,j); p3=x(i+1),y(j+1),z(i+1,j+1); p4=x(i),y(j+1),z(i,j+1); p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1); p14=norm(p4-p1); p34=norm(p4-p3); if all(z(i,j),z(i+1,j),z(i+1,j+1)=12) z1=(p12+p23+p13)/2; s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13);

8、 s=s+s1; end if all(z(i,j),z(i+1,j+1),z(i,j+1)=12) z2=(p13+p14+p34)/2; s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34); s=s+s2; end endends宋哲 2011-08-17 9:32:57%图3的Matlab程序如下clc, clearload pdataa=load( );x0=0:5:4000; y0=0:5:3000;v=-49,1,12,44,75,107,171,184:80:357,357:100:614,614;c=contour(x0,y0,a,v); clabel(

9、c)hold onplot(P1(1),P1(2),S)plot(P2(1),P2(2),D), plot(P3(1),P3(2),D), plot(P4(1),P4(2),D)b=a; for i=1:801 for j=1:601 txy=(i-1)*5;(j-1)*5;b(i,j); td=norm(txy-Pz1); if td=995 & td=495 & td=495 & td=495 & td=12) z1=(p12+p23+p13)/2; s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13); s=s+s1; end if any(z(i,j),z(i+1

10、,j+1),z(i,j+1)=12) z2=(p13+p14+p34)/2; s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34); s=s+s2; end endends宋哲 2011-08-17 16:41:37%求4个特殊点的Matlab程序clc, cleara=load();a=a; x0=0:5:4000; y0=0:5:3000;i=2500:5:3700/5; j=800:5:2100/5; %求358高地的大致范围的地址z0=a(i,j); %求358高地的大致高程zz1=max(max(z0) %求358高地的最高点高程i1,j1=find(a=zz1

11、); %求358高地在全部数据矩阵a中的地址P1=(i1-1)*5; (j1-1)*5 %求358高地的x坐标和y坐标zz2=max(max(a)i2,j2=find(a=zz2);P2=(i2-1)*5; (j2-1)*5 %求615高地的最高点x坐标和y坐标a0=a(800:5:1100/5,600:5:800/5);zz3=max(max(a0) %求185高地的最高点高程i3,j3=find(a=zz3);P3=(i3-1)*5; (j3-1)*5 %求185高地的最高点x坐标和y坐标z4=a(2100:5:2650/5,170:5:580/5); zz4=max(max(z4)i4,

12、j4=find(a=zz4);P4=(i4-1)*5; (j4-1)*5Pz1=P1;zz1+10; Pz2=P2; zz2+3; Pz3=P3; zz3+3; Pz4=P4; zz4+3;save pdata P1 P2 P3 P4 Pz1 Pz2 Pz3 Pz4宋哲 2011-08-17 16:41:48%图3的Matlab程序如下clc, clearload pdataa=load();x0=0:5:4000; y0=0:5:3000;v=-49,1,12,44,75,107,171,184:80:357,357:100:614,614;c=contour(x0,y0,a,v); clabel(c)hold on plot(P1(1),P1(2),S)plot(P2(1),P2(2),D), plot(P3(1),P3(2),D), plot(P4(1),P4(2),D)b=a; for i=1:801 for j=1:601 txy=(i-1)*5;(j-1)*5;b(i,j); td=norm(txy-Pz1); if td=995 & td=495 & td=495 & td=495 & td=500 plot(i-1)*5,(j-1)*5,o) end endend