arma模型程序.docx
《arma模型程序.docx》由会员分享，可在线阅读，更多相关《arma模型程序.docx（9页珍藏版）》请在冰豆网上搜索。
arma模型程序
ARMA模型程序（总6页）
clear,clc
closeall
set（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）;%画图并且取得二维图形中最大值和最小值的索引号
gridon
set（gca,'layer','top'）%将栅格放在最上层，防止遮盖
axisequal%使每一格的刻度相同，同时也使长度一样
xlabel（'-\pi\leq\theta\leq\pi','fontsize',20）
ylabel（'tan（sin（\theta））-sin（tan（\theta））','fontsize',20）
title（'\it{-\pi\leq\theta\leq\pi的三角函数图}','fontsize',20）
text（0,0,'\bullet\it原点','fontsize',25）%文本的精确定位
%标注曲线的最大值和最小值
x=get（h,'xdata'）;%获得二维曲线的数据
y=get（h,'ydata'）;
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/1611:
07:
42
%AR,MA模型
clear,clc
closeall
a=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-timeIDPOLYmodel:
%A（q）=1-q^-1-q^-2
b15=bhat{1}（end）
bhat2=predict（cs1,[b;b15;-1]）
宋哲2011/8/1611:
27:
00
%用逻辑斯特预测
clear,clc
a=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,clear
a=load（''）;
amin=min（min（a））,amax=max（max（a））
x0=0:
5:
4000;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-179:
32:
44
clc,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;
fori=1:
m-1
forj=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）;
ifall（[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））;
s=s+s1;
end
ifall（[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
end
end
s
宋哲2011-08-179:
32:
57
%图3的Matlab程序如下
clc,clear
loadpdata
a=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）
holdon
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';
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz1）;
iftd>=995&td<=1000
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz2）;
iftd>=495&td<=500
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz3）;
iftd>=495&td<=500
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz4）;
iftd>=495&td<=500
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
宋哲2011-08-1716:
11:
05
clc,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;
fori=1:
m-1
forj=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）;
ifany（[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））;
s=s+s1;
end
ifany（[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
end
end
s
宋哲2011-08-1716:
41:
37
%求4个特殊点的Matlab程序
clc,clear
a=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）;%求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,j4]=find（a==zz4）;
P4=[（i4-1）*5;（j4-1）*5]
Pz1=[P1;zz1+10];Pz2=[P2;zz2+3];Pz3=[P3;zz3+3];Pz4=[P4;zz4+3];
savepdataP1P2P3P4Pz1Pz2Pz3Pz4
宋哲2011-08-1716:
41:
48
%图3的Matlab程序如下
clc,clear
loadpdata
a=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）
holdon
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';
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz1）;
iftd>=995&td<=1000
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz2）;
iftd>=495&td<=500
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz3）;
iftd>=495&td<=500
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
fori=1:
801
forj=1:
601
txy=[（i-1）*5;（j-1）*5;b（i,j）];
td=norm（txy-Pz4）;
iftd>=495&td<=500
plot（（i-1）*5,（j-1）*5,'o'）
end
end
end
升级会员 